MODULE module_microphysics_driver
CONTAINS

SUBROUTINE microphysics_driver( &
                       th, rho, pi_phy, p &
                      ,ht, dz8w, p8w, dt,dx,dy &
                      ,mp_physics, spec_zone &
                      ,specified, channel_switch &
                      ,warm_rain &
                      ,t8w &
                      ,chem_opt, progn &
                      ,cldfra, cldfra_old, exch_h, nsource &
                      ,qlsink, precr, preci, precs, precg &
                      ,xland,snowh,itimestep &
                      ,f_ice_phy,f_rain_phy,f_rimef_phy &
                      ,lowlyr,sr, id &
                      ,ids,ide, jds,jde, kds,kde &
                      ,ims,ime, jms,jme, kms,kme &
                      ,ips,ipe, jps,jpe, kps,kpe &
                      ,i_start,i_end,j_start,j_end,kts,kte &
                      ,num_tiles, naer &

                      ,dlf,dlf2,t_phy,p_hyd,p8w_hyd,tke_pbl,z_at_w,qfx &
                      ,rliq,turbtype3d,smaw3d,wsedl3d,cldfra_old_mp &
                      ,cldfra_mp,cldfra_mp_all,lradius,iradius &
                      ,cldfrai,cldfral,cldfra_conv &
                      ,alt &
                      ,accum_mode,aitken_mode,coarse_mode &
                      ,icwmrsh3d,icwmrdp3d,shfrc3d,cmfmc3d,cmfmc2_3d &
                      ,config_flags,fnm,fnp,rh_old_mp,lcd_old_mp &

                      ,chem &
                      ,qme3d,prain3d,nevapr3d,rate1ord_cw2pr_st3d &
                      ,dgnum4D,dgnumwet4D &

                      ,qv_curr,qc_curr,qr_curr,qi_curr,qs_curr,qg_curr &
                      ,qndrop_curr,qni_curr,qh_curr,qnh_curr &
                      ,qzr_curr,qzi_curr,qzs_curr,qzg_curr,qzh_curr &
                      ,qns_curr,qnr_curr,qng_curr,qnn_curr,qnc_curr &
                      ,qvolg_curr &
                      ,f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop,f_qni &
                      ,f_qns,f_qnr,f_qng,f_qnc,f_qnn,f_qh,f_qnh &
                      , f_qzr,f_qzi,f_qzs,f_qzg,f_qzh &
                      ,f_qvolg &
                      ,qrcuten, qscuten, qicuten, mu &
                      ,qt_curr,f_qt &
                      ,mp_restart_state,tbpvs_state,tbpvs0_state &
                      ,hail,ice2 &

                      ,w ,z &
                      ,rainnc, rainncv &
                      ,snownc, snowncv &
                      ,hailnc, hailncv &
                      ,graupelnc, graupelncv &

                      ,rainprod, evapprod &
                      ,qv_b4mp, qc_b4mp, qi_b4mp, qs_b4mp &

                      ,refl_10cm &

                      ,ri_curr &
                      ,diagflag, do_radar_ref &
                      ,re_cloud, re_ice, re_snow &
                      ,has_reqc, has_reqi, has_reqs &
                                                   )
   USE module_state_description, ONLY : &
                     KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME &
                    ,WSM6SCHEME, ETAMPNEW, ETAMPOLD, THOMPSON, MORR_TWO_MOMENT &
                    ,GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, NSSL_2MOM, NSSL_2MOMCCN &
                    ,NSSL_1MOM,NSSL_1MOMLFO &
                    ,MILBRANDT2MOM , CAMMGMPSCHEME
   USE module_model_constants
   USE module_wrf_error
   USE module_configure, only: grid_config_rec_type
   USE module_state_description, only: num_chem
   USE modal_aero_data, only: ntot_amode_cam_mam => ntot_amode
   USE module_configure, only: &
       p_na_a01, p_na_a02, p_na_a03, p_na_a04, &
       p_na_a05, p_na_a06, p_na_a07, p_na_a08, &
       p_cl_a01, p_cl_a02, p_cl_a03, p_cl_a04, &
       p_cl_a05, p_cl_a06, p_cl_a07, p_cl_a08, &
       p_oin_a01, p_oin_a02, p_oin_a03, p_oin_a04, &
       p_oin_a05, p_oin_a06, p_oin_a07, p_oin_a08, &
       p_num_a01, p_num_a02, p_num_a03, p_num_a04, &
       p_num_a05, p_num_a06, p_num_a07, p_num_a08
   USE module_mp_kessler
   USE module_mp_lin
   USE module_mp_sbu_ylin
   USE module_mp_wsm3
   USE module_mp_wsm5
   USE module_mp_wsm6
   USE module_mp_etanew
   USE module_mp_etaold
   USE module_mp_thompson
   USE module_mp_gsfcgce
   USE module_mp_morr_two_moment
   USE module_mp_wdm5
   USE module_mp_wdm6
   USE module_mp_milbrandt2mom
   USE module_mp_cammgmp_driver, ONLY: CAMMGMP
   USE module_mp_nssl_2mom
   USE module_mp_HWRF
   USE module_mixactivate, only: prescribe_aerosol_mixactivate
   USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
   USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
   IMPLICIT NONE
   TYPE(grid_config_rec_type), INTENT(IN ) , OPTIONAL :: config_flags
   INTEGER, INTENT(IN ) :: mp_physics
   LOGICAL, INTENT(IN ) :: specified
   INTEGER, OPTIONAL, INTENT(IN ) :: chem_opt, progn
   INTEGER, OPTIONAL, INTENT(IN ) :: hail, ice2
   INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
   INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
   INTEGER, OPTIONAL, INTENT(IN ) :: ips,ipe, jps,jpe, kps,kpe
   INTEGER, INTENT(IN ) :: kts,kte
   INTEGER, INTENT(IN ) :: itimestep,num_tiles,spec_zone
   INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
     & i_start,i_end,j_start,j_end
   LOGICAL, INTENT(IN ) :: warm_rain
   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
         INTENT(INOUT) :: th
   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
         INTENT(IN ) :: &
                                                                 rho, &
                                                                dz8w, &
                                                                 p8w, &
                                                              pi_phy, &
                                                                   p
   REAL,INTENT(IN), OPTIONAL ::accum_mode,aitken_mode,coarse_mode
   REAL , DIMENSION( kms:kme ) , &
        INTENT(IN ) , OPTIONAL :: fnm, &
                                                                fnp
   REAL, DIMENSION( ims:ime, jms:jme ), &
        INTENT(IN), OPTIONAL :: &
                                                                 qfx, &
                                                                 rliq
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
      INTENT(IN), OPTIONAL :: &
                                                                 dlf, &
                                                                dlf2, &
                                                               t_phy, &
                                                               p_hyd, &
                                                             p8w_hyd, &
                                                              z_at_w, &
                                                             tke_pbl, &
                                                          turbtype3d, &
                                                              smaw3d, &
                                                                 alt, &
                                                           icwmrsh3d, &
                                                           icwmrdp3d, &
                                                             shfrc3d, &
                                                             cmfmc3d, &
                                                           cmfmc2_3d
 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme,ntot_amode_cam_mam ), &
        INTENT(IN) :: &
                                                             dgnum4D, &
                                                          dgnumwet4D
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
      INTENT(INOUT) , OPTIONAL :: &
                                                       cldfra_old_mp, &
                                                           rh_old_mp, &
                                                          lcd_old_mp
 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem), &
      INTENT(INOUT) :: &
                                                                 chem
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
      INTENT(INOUT) , OPTIONAL:: &
                                                            wsedl3d, &
                                                          cldfra_mp, &
                                                      cldfra_mp_all, &
                                                            cldfrai, &
                                                            cldfral, &
                                                            lradius, &
                                                            iradius, &
                                                            cldfra_conv
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
      INTENT(INOUT), OPTIONAL :: &
                                                              qme3d, &
                                                            prain3d, &
                                                           nevapr3d, &
                                                rate1ord_cw2pr_st3d
   REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
                                     F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
   REAL, OPTIONAL, INTENT(OUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
         qlsink, &
         precr, &
         preci, &
         precs, &
         precg
   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: XLAND
   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN), OPTIONAL :: SNOWH
   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: SR
   REAL, INTENT(IN ) :: dt,dx,dy
   INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: LOWLYR
   REAL, OPTIONAL, DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(OUT) :: refl_10cm
   LOGICAL, OPTIONAL, INTENT(IN ) :: channel_switch
   REAL, OPTIONAL, INTENT(INOUT ) :: naer
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
         OPTIONAL, &
         INTENT(INOUT ) :: &
                  w, z, t8w &
                 ,cldfra, cldfra_old, exch_h &
                 ,qv_curr,qc_curr,qr_curr,qi_curr,qs_curr,qg_curr &
                 ,qt_curr,qndrop_curr,qni_curr,qh_curr,qnh_curr &
                 ,qns_curr,qnr_curr,qng_curr,qnn_curr,qnc_curr &
                 ,qzr_curr,qzi_curr,qzs_curr,qzg_curr,qzh_curr &
                 ,qvolg_curr
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
         OPTIONAL, &
         INTENT(IN) :: qrcuten, qscuten, qicuten
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
         INTENT(INOUT) :: rainprod, evapprod
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
         INTENT(INOUT) :: qv_b4mp, qc_b4mp, qi_b4mp, qs_b4mp
   REAL, DIMENSION( ims:ime, jms:jme ), &
         OPTIONAL, &
         INTENT(IN) :: mu
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
         OPTIONAL, &
         INTENT(INOUT) :: ri_curr
   REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), &
         OPTIONAL, &
         INTENT(OUT ) :: &
                  nsource
   REAL, DIMENSION( ims:ime , jms:jme ), &
         INTENT(INOUT), &
         OPTIONAL :: &
                                                           RAINNC &
                                                         ,RAINNCV &
                                                          ,SNOWNC &
                                                         ,SNOWNCV &
                                                       ,GRAUPELNC &
                                                      ,GRAUPELNCV &
                                                          ,HAILNC &
                                                         ,HAILNCV
   INTEGER,OPTIONAL,INTENT(IN ) :: id
   REAL , DIMENSION( ims:ime , jms:jme ) , OPTIONAL , &
         INTENT(IN) :: ht
   REAL, DIMENSION (:), OPTIONAL, INTENT(INOUT) :: mp_restart_state &
                                         ,tbpvs_state,tbpvs0_state
   LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop,f_qni,f_qt &
                       ,f_qns,f_qnr,f_qng,f_qnn,f_qnc,f_qh,f_qnh,f_qzr &
                       ,f_qzi,f_qzs,f_qzg,f_qzh,f_qvolg
   LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
   INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
   REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: &
                 re_cloud, re_ice, re_snow
   INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs
   INTEGER :: i,j,k,its,ite,jts,jte,ij,sz,n
   LOGICAL :: channel
   REAL :: z0, z1, z2, w1, w2
   channel = .FALSE.
   IF ( PRESENT ( channel_switch ) ) channel = channel_switch
   if (mp_physics .eq. 0) return
   IF( specified ) THEN
     sz = spec_zone
   ELSE
     sz = 0
   ENDIF
   !$OMP PARALLEL DO &
   !$OMP PRIVATE ( ij, its, ite, jts, jte, i,j,k,n )
   DO ij = 1 , num_tiles
       IF (channel) THEN
         its = max(i_start(ij),ids)
         ite = min(i_end(ij),ide-1)
       ELSE
         its = max(i_start(ij),ids+sz)
         ite = min(i_end(ij),ide-1-sz)
       ENDIF
       jts = max(j_start(ij),jds+sz)
       jte = min(j_end(ij),jde-1-sz)
       IF( PRESENT(qlsink) ) qlsink(its:ite,kts:kte,jts:jte) = 0.
       IF( PRESENT(precr ) ) precr(its:ite,kts:kte,jts:jte) = 0.
       IF( PRESENT(preci ) ) preci(its:ite,kts:kte,jts:jte) = 0.
       IF( PRESENT(precs ) ) precs(its:ite,kts:kte,jts:jte) = 0.
       IF( PRESENT(precg ) ) precg(its:ite,kts:kte,jts:jte) = 0.
       IF( PRESENT(chem_opt) .AND. PRESENT(progn) ) THEN
       IF( chem_opt==0 .AND. progn==1 .AND. (mp_physics==LINSCHEME .OR. mp_physics==MORR_TWO_MOMENT)) THEN
          IF( PRESENT( QNDROP_CURR ) ) THEN
             CALL wrf_debug ( 100 , 'microphysics_driver: calling prescribe_aerosol_mixactivate' )
             call prescribe_aerosol_mixactivate ( &
                  id, itimestep, dt, naer, &
                  rho, th, pi_phy, w, cldfra, cldfra_old, &
                  z, dz8w, p8w, t8w, exch_h, &
                  qv_curr, qc_curr, qi_curr, qndrop_curr, &
                  nsource, &
                  ids,ide, jds,jde, kds,kde, &
                  ims,ime, jms,jme, kms,kme, &
                  its,ite, jts,jte, kts,kte, &
                  F_QC=f_qc, F_QI=f_qi )
          END IF
       ELSE IF( progn==1 .AND. mp_physics/=LINSCHEME .AND. mp_physics/=MORR_TWO_MOMENT) THEN
             call wrf_error_fatal3("<stdin>",559,&
             "SETTINGS ERROR: Prognostic cloud droplet number can only be used with the mp_physics=LINSCHEME or MORRISON.")
       END IF
       END IF
     micro_select: SELECT CASE(mp_physics)
        CASE (KESSLERSCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling kessler' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT( QC_CURR ) .AND. &
                                           PRESENT( QR_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                                           PRESENT( Z )) THEN
               CALL kessler( &
                  T=th &
                 ,QV=qv_curr &
                 ,QC=qc_curr &
                 ,QR=qr_curr &
                 ,RHO=rho, PII=pi_phy,DT_IN=dt, Z=z, XLV=xlv, CP=cp &
                 ,EP2=ep_2,SVP1=svp1,SVP2=svp2 &
                 ,SVP3=svp3,SVPT0=svpt0,RHOWATER=rhowater &
                 ,DZ8W=dz8w &
                 ,RAINNC=rainnc,RAINNCV=rainncv &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",587,&
'arguments not present for calling kessler' )
             ENDIF
        CASE (THOMPSON)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling thompson' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. &
                  PRESENT( QNR_CURR) .AND. PRESENT ( QNI_CURR) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
                 qv_b4mp(its:ite,kts:kte,jts:jte) = qv_curr(its:ite,kts:kte,jts:jte)
                 qc_b4mp(its:ite,kts:kte,jts:jte) = qc_curr(its:ite,kts:kte,jts:jte)
                 qi_b4mp(its:ite,kts:kte,jts:jte) = qi_curr(its:ite,kts:kte,jts:jte)
                 qs_b4mp(its:ite,kts:kte,jts:jte) = qs_curr(its:ite,kts:kte,jts:jte)
             CALL mp_gt_driver( &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QG=qg_curr, &
                     NI=qni_curr, &
                     NR=qnr_curr, &
                     TH=th, &
                     PII=pi_phy, &
                     P=p, &
                     DZ=dz8w, &
                     DT_IN=dt, &
                     ITIMESTEP=itimestep, &
                     RAINNC=RAINNC, &
                     RAINNCV=RAINNCV, &
                     SNOWNC=SNOWNC, &
                     SNOWNCV=SNOWNCV, &
                     GRAUPELNC=GRAUPELNC, &
                     GRAUPELNCV=GRAUPELNCV, &
                     SR=SR, &
                     RAINPROD=rainprod, &
                     EVAPPROD=evapprod, &
                     REFL_10CM=refl_10cm, &
                     diagflag=diagflag, &
                     do_radar_ref=do_radar_ref, &
                     re_cloud=re_cloud, &
                     re_ice=re_ice, &
                     re_snow=re_snow, &
                     has_reqc=has_reqc, &
                     has_reqi=has_reqi, &
                     has_reqs=has_reqs, &
                 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, &
                 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, &
                 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte)
             ELSE
                CALL wrf_error_fatal3("<stdin>",640,&
'arguments not present for calling thompson_et_al' )
             ENDIF
    CASE (MORR_TWO_MOMENT)
         CALL wrf_debug(100, 'microphysics_driver: calling morrison two moment')
         IF (PRESENT (QV_CURR) .AND. PRESENT (QC_CURR) .AND. &
             PRESENT (QR_CURR) .AND. PRESENT (QI_CURR) .AND. &
         PRESENT (QS_CURR) .AND. PRESENT (QG_CURR) .AND. &
         PRESENT (QR_CURR) .AND. PRESENT (QI_CURR) .AND. &
         PRESENT (QNS_CURR) .AND. PRESENT (QNI_CURR).AND. &
         PRESENT (QNR_CURR) .AND. PRESENT (QNG_CURR).AND. &
         PRESENT (MU) .AND. PRESENT (QSCUTEN).AND. &
         PRESENT (QRCUTEN) .AND. PRESENT (QICUTEN).AND. &
         PRESENT (RAINNC ) .AND. PRESENT (RAINNCV) .AND. &
         PRESENT (Z ) .AND.PRESENT ( W ) ) THEN
         CALL mp_morr_two_moment( &
                     ITIMESTEP=itimestep, &
                     TH=th, &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QG=qg_curr, &
                     NI=qni_curr, &
                     NS=qns_curr, &
                     NR=qnr_curr, &
                     NG=qng_curr, &
                     RHO=rho, &
                     PII=pi_phy, &
                     P=p, &
                     DT_IN=dt, &
                     DZ=dz8w, &
                     HT=ht, &
                     W=w &
                    ,RAINNC=RAINNC &
                    ,RAINNCV=RAINNCV &
                    ,SNOWNC=SNOWNC &
                    ,SNOWNCV=SNOWNCV &
                    ,GRAUPELNC=GRAUPELNC &
                    ,GRAUPELNCV=GRAUPELNCV &
                    ,SR=SR &
                    ,REFL_10CM=refl_10cm &
                    ,diagflag=diagflag &
                    ,do_radar_ref=do_radar_ref &
                    ,qrcuten=qrcuten &
                    ,qscuten=qscuten &
                    ,qicuten=qicuten &
                    ,mu=mu &
                    ,F_QNDROP=f_qndrop &
                 ,QNDROP=qndrop_curr &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                 ,QLSINK=qlsink &
                 ,PRECR=precr,PRECI=preci,PRECS=precs,PRECG=precg &
                                                                    )
        ELSE
           Call wrf_error_fatal3("<stdin>",700,&
'arguments not present for calling morrison two moment')
        ENDIF
    CASE (MILBRANDT2MOM)
         CALL wrf_debug(100, 'microphysics_driver: calling milbrandt2mom')
         IF (PRESENT (QV_CURR) .AND. &
             PRESENT (QC_CURR) .AND. PRESENT (QNC_CURR) .AND. &
             PRESENT (QR_CURR) .AND. PRESENT (QNR_CURR) .AND. &
             PRESENT (QI_CURR) .AND. PRESENT (QNI_CURR) .AND. &
             PRESENT (QS_CURR) .AND. PRESENT (QNS_CURR) .AND. &
             PRESENT (QG_CURR) .AND. PRESENT (QNG_CURR) .AND. &
             PRESENT (QH_CURR) .AND. PRESENT (QNH_CURR) .AND. &
             PRESENT (RAINNC ) .AND. PRESENT (RAINNCV) .AND. &
             PRESENT (SNOWNC ) .AND. PRESENT (SNOWNCV) .AND. &
             PRESENT (HAILNC ) .AND. PRESENT (HAILNCV) .AND. &
             PRESENT (GRAUPELNC).AND.PRESENT (GRAUPELNCV).AND. &
             PRESENT (Z ) .AND. PRESENT ( W ) ) THEN
         CALL mp_milbrandt2mom_driver( &
                     ITIMESTEP=itimestep, &
                     p8w=p8w, &
                     TH=th, &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QG=qg_curr, &
                     QH=qh_curr, &
                     NC=qnc_curr, &
                     NR=qnr_curr, &
                     NI=qni_curr, &
                     NS=qns_curr, &
                     NG=qng_curr, &
                     NH=qnh_curr, &
                     PII=pi_phy, &
                     P=p, &
                     DT_IN=dt, &
                     DZ=dz8w, &
                     W=w, &
                     RAINNC = RAINNC, &
                     RAINNCV = RAINNCV, &
                     SNOWNC = SNOWNC, &
                     SNOWNCV = SNOWNCV, &
                     HAILNC = HAILNC, &
                     HAILNCV = HAILNCV, &
                     GRPLNC = GRAUPELNC, &
                     GRPLNCV = GRAUPELNCV, &
                     SR=SR, &
                     Zet = refl_10cm, &
                  chem = chem, &
                  IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, &
                  IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, &
                  ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
        ELSE
           Call wrf_error_fatal3("<stdin>",760,&
'arguments not present for calling milbrandt2mom')
        ENDIF
    CASE (NSSL_1MOM)
         CALL wrf_debug(100, 'microphysics_driver: calling nssl1mom')
         IF (PRESENT (QV_CURR) .AND. &
             PRESENT (QC_CURR) .AND. &
             PRESENT (QR_CURR) .AND. &
             PRESENT (QI_CURR) .AND. &
             PRESENT (QS_CURR) .AND. &
             PRESENT (QG_CURR) .AND. &
             PRESENT (QH_CURR) .AND. &
             PRESENT (RAINNC ) .AND. PRESENT (RAINNCV) .AND. &
             PRESENT (SNOWNC ) .AND. PRESENT (SNOWNCV) .AND. &
             PRESENT (HAILNC ) .AND. PRESENT (HAILNCV) .AND. &
             PRESENT (GRAUPELNC).AND.PRESENT (GRAUPELNCV).AND. &
             PRESENT (Z ) .AND. PRESENT ( W ) .AND. &
             PRESENT (QVOLG_CURR) ) THEN
         CALL nssl_2mom_driver( &
                     ITIMESTEP=itimestep, &
                     TH=th, &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QH=qg_curr, &
                     QHL=qh_curr, &
                     CCW=qnc_curr, &
                     CRW=qnr_curr, &
                     CCI=qni_curr, &
                     CSW=qns_curr, &
                     CHW=qng_curr, &
                     CHL=qnh_curr, &
                     VHW=qvolg_curr, &
                     PII=pi_phy, &
                     P=p, &
                     W=w, &
                     DZ=dz8w, &
                     DTP=dt, &
                     DN=rho, &
                     RAINNC = RAINNC, &
                     RAINNCV = RAINNCV, &
                     SNOWNC = SNOWNC, &
                     SNOWNCV = SNOWNCV, &
                     HAILNC = HAILNC, &
                     HAILNCV = HAILNCV, &
                     GRPLNC = GRAUPELNC, &
                     GRPLNCV = GRAUPELNCV, &
                     SR=SR, &
                     dbz = refl_10cm, &
                     diagflag = diagflag, &
                  IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, &
                  IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, &
                  ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
        ELSE
           Call wrf_error_fatal3("<stdin>",869,&
'arguments not present for calling nssl_2mom')
        ENDIF
    CASE (NSSL_1MOMLFO)
         CALL wrf_debug(100, 'microphysics_driver: calling nssl1mom')
         IF (PRESENT (QV_CURR) .AND. &
             PRESENT (QC_CURR) .AND. &
             PRESENT (QR_CURR) .AND. &
             PRESENT (QI_CURR) .AND. &
             PRESENT (QS_CURR) .AND. &
             PRESENT (QG_CURR) .AND. &
             PRESENT (RAINNC ) .AND. PRESENT (RAINNCV) .AND. &
             PRESENT (SNOWNC ) .AND. PRESENT (SNOWNCV) .AND. &
             PRESENT (GRAUPELNC).AND.PRESENT (GRAUPELNCV).AND. &
             PRESENT (Z ) .AND. PRESENT ( W ) ) THEN
         CALL nssl_2mom_driver( &
                     ITIMESTEP=itimestep, &
                     TH=th, &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QH=qg_curr, &
                     PII=pi_phy, &
                     P=p, &
                     W=w, &
                     DZ=dz8w, &
                     DTP=dt, &
                     DN=rho, &
                     RAINNC = RAINNC, &
                     RAINNCV = RAINNCV, &
                     SNOWNC = SNOWNC, &
                     SNOWNCV = SNOWNCV, &
                     GRPLNC = GRAUPELNC, &
                     GRPLNCV = GRAUPELNCV, &
                     SR=SR, &
                     dbz = refl_10cm, &
                     diagflag = diagflag, &
                  IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, &
                  IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, &
                  ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
        ELSE
           Call wrf_error_fatal3("<stdin>",917,&
'arguments not present for calling nssl_2mom')
        ENDIF
    CASE (NSSL_2MOM)
         CALL wrf_debug(100, 'microphysics_driver: calling nssl2mom')
         IF (PRESENT (QV_CURR) .AND. &
             PRESENT (QC_CURR) .AND. PRESENT (QNC_CURR) .AND. &
             PRESENT (QR_CURR) .AND. PRESENT (QNR_CURR) .AND. &
             PRESENT (QI_CURR) .AND. PRESENT (QNI_CURR) .AND. &
             PRESENT (QS_CURR) .AND. PRESENT (QNS_CURR) .AND. &
             PRESENT (QG_CURR) .AND. PRESENT (QNG_CURR) .AND. &
             PRESENT (QH_CURR) .AND. PRESENT (QNH_CURR) .AND. &
             PRESENT (RAINNC ) .AND. PRESENT (RAINNCV) .AND. &
             PRESENT (SNOWNC ) .AND. PRESENT (SNOWNCV) .AND. &
             PRESENT (HAILNC ) .AND. PRESENT (HAILNCV) .AND. &
             PRESENT (GRAUPELNC).AND.PRESENT (GRAUPELNCV).AND. &
             PRESENT (Z ) .AND. PRESENT ( W ) .AND. &
             PRESENT (QVOLG_CURR) ) THEN
         CALL nssl_2mom_driver( &
                     ITIMESTEP=itimestep, &
                     TH=th, &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QH=qg_curr, &
                     QHL=qh_curr, &
                     CCW=qnc_curr, &
                     CRW=qnr_curr, &
                     CCI=qni_curr, &
                     CSW=qns_curr, &
                     CHW=qng_curr, &
                     CHL=qnh_curr, &
                     VHW=qvolg_curr, &
                     PII=pi_phy, &
                     P=p, &
                     W=w, &
                     DZ=dz8w, &
                     DTP=dt, &
                     DN=rho, &
                     RAINNC = RAINNC, &
                     RAINNCV = RAINNCV, &
                     SNOWNC = SNOWNC, &
                     SNOWNCV = SNOWNCV, &
                     HAILNC = HAILNC, &
                     HAILNCV = HAILNCV, &
                     GRPLNC = GRAUPELNC, &
                     GRPLNCV = GRAUPELNCV, &
                     SR=SR, &
                     dbz = refl_10cm, &
                     diagflag = diagflag, &
                  IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, &
                  IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, &
                  ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
        ELSE
           Call wrf_error_fatal3("<stdin>",977,&
'arguments not present for calling nssl_2mom')
        ENDIF
    CASE (NSSL_2MOMCCN)
         CALL wrf_debug(100, 'microphysics_driver: calling nssl_2momccn')
         IF (PRESENT (QV_CURR) .AND. &
             PRESENT (QC_CURR) .AND. PRESENT (QNC_CURR) .AND. &
             PRESENT (QR_CURR) .AND. PRESENT (QNR_CURR) .AND. &
             PRESENT (QI_CURR) .AND. PRESENT (QNI_CURR) .AND. &
             PRESENT (QS_CURR) .AND. PRESENT (QNS_CURR) .AND. &
             PRESENT (QG_CURR) .AND. PRESENT (QNG_CURR) .AND. &
             PRESENT (QH_CURR) .AND. PRESENT (QNH_CURR) .AND. &
             PRESENT (RAINNC ) .AND. PRESENT (RAINNCV) .AND. &
             PRESENT (SNOWNC ) .AND. PRESENT (SNOWNCV) .AND. &
             PRESENT (HAILNC ) .AND. PRESENT (HAILNCV) .AND. &
             PRESENT (GRAUPELNC).AND.PRESENT (GRAUPELNCV).AND. &
             PRESENT (Z ) .AND. PRESENT ( W ) .AND. &
             PRESENT (QVOLG_CURR) .AND. PRESENT( QNN_CURR ) ) THEN
         CALL nssl_2mom_driver( &
                     ITIMESTEP=itimestep, &
                     TH=th, &
                     QV=qv_curr, &
                     QC=qc_curr, &
                     QR=qr_curr, &
                     QI=qi_curr, &
                     QS=qs_curr, &
                     QH=qg_curr, &
                     QHL=qh_curr, &
                     CCW=qnc_curr, &
                     CRW=qnr_curr, &
                     CCI=qni_curr, &
                     CSW=qns_curr, &
                     CHW=qng_curr, &
                     CHL=qnh_curr, &
                     VHW=qvolg_curr, &
                     cn=qnn_curr, &
                     PII=pi_phy, &
                     P=p, &
                     W=w, &
                     DZ=dz8w, &
                     DTP=dt, &
                     DN=rho, &
                     RAINNC = RAINNC, &
                     RAINNCV = RAINNCV, &
                     SNOWNC = SNOWNC, &
                     SNOWNCV = SNOWNCV, &
                     HAILNC = HAILNC, &
                     HAILNCV = HAILNCV, &
                     GRPLNC = GRAUPELNC, &
                     GRPLNCV = GRAUPELNCV, &
                     SR=SR, &
                     dbz = refl_10cm, &
                     diagflag = diagflag, &
                  IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, &
                  IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, &
                  ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
        ELSE
           Call wrf_error_fatal3("<stdin>",1038,&
'arguments not present for calling nssl_2momccn')
        ENDIF
        CASE (GSFCGCESCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling GSFCGCE' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                  PRESENT( HAIL ) .AND. PRESENT ( ICE2 ) .AND. &
                  PRESENT( Z ) .AND. PRESENT ( W ) ) THEN
               CALL gsfcgce( &
                  TH=th &
                 ,QV=qv_curr &
                 ,QL=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,RHO=rho, PII=pi_phy, P=p, DT_IN=dt, Z=z &
                 ,HT=ht, DZ8W=dz8w, GRAV=G &
                 ,RHOWATER=rhowater, RHOSNOW=rhosnow &
                 ,ITIMESTEP=itimestep &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                 ,RAINNC=rainnc, RAINNCV=rainncv &
                 ,SNOWNC=snownc, SNOWNCV=snowncv ,SR=sr &
                 ,GRAUPELNC=graupelnc ,GRAUPELNCV=graupelncv &
                 ,REFL_10CM=refl_10cm &
                 ,diagflag=diagflag &
                 ,do_radar_ref=do_radar_ref &
                 ,F_QG=f_qg &
                 ,QG=qg_curr &
                 ,IHAIL=hail, ICE2=ice2 &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1082,&
'arguments not present for calling GSFCGCE' )
             ENDIF
        CASE (LINSCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling lin_et_al' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                  PRESENT( Z ) ) THEN
               CALL lin_et_al( &
                  TH=th &
                 ,QV=qv_curr &
                 ,QL=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,QLSINK=qlsink &
                 ,RHO=rho, PII=pi_phy, P=p, DT_IN=dt, Z=z &
                 ,HT=ht, DZ8W=dz8w, GRAV=G, CP=cp &
                 ,RAIR=r_d, RVAPOR=R_v &
                 ,XLS=xls, XLV=xlv, XLF=xlf &
                 ,RHOWATER=rhowater, RHOSNOW=rhosnow &
                 ,EP2=ep_2,SVP1=svp1,SVP2=svp2 &
                 ,SVP3=svp3,SVPT0=svpt0 &
                 ,RAINNC=rainnc, RAINNCV=rainncv &
                 ,SNOWNC=snownc, SNOWNCV=snowncv &
                 ,GRAUPELNC=graupelnc, GRAUPELNCV=graupelncv, SR=sr &
                 ,REFL_10CM=refl_10cm &
                 ,diagflag=diagflag &
                 ,do_radar_ref=do_radar_ref &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                 ,PRECR=precr,PRECI=preci,PRECS=precs,PRECG=precg &
                 ,F_QG=f_qg, F_QNDROP=f_qndrop &
                 ,QG=qg_curr &
                 ,QNDROP=qndrop_curr &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1123,&
'arguments not present for calling lin_et_al' )
             ENDIF
       CASE (SBU_YLINSCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling sbu_ylin' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. &
                  PRESENT( RI_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                  PRESENT( Z ) ) THEN
               CALL sbu_ylin( &
                  TH=th &
                 ,QV=qv_curr &
                 ,QL=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,RI3D=ri_curr &
                 ,RHO=rho, PII=pi_phy, P=p, DT_IN=dt, Z=z &
                 ,HT=ht, DZ8W=dz8w &
                 ,RAINNC=rainnc, RAINNCV=rainncv &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                     )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1165,&
'arguments not present for calling sbu_ylin' )
             ENDIF
        CASE (WSM3SCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling wsm3' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                  PRESENT( W ) ) THEN
             CALL wsm3( &
                  TH=th &
                 ,Q=qv_curr &
                 ,QCI=qc_curr &
                 ,QRS=qr_curr &
                 ,W=w,DEN=rho,PII=pi_phy,P=p,DELZ=dz8w &
                 ,DELT=dt,G=g,CPD=cp,CPV=cpv &
                 ,RD=r_d,RV=r_v,T0C=svpt0 &
                 ,EP1=ep_1, EP2=ep_2, QMIN=epsilon &
                 ,XLS=xls, XLV0=xlv, XLF0=xlf &
                 ,DEN0=rhoair0, DENR=rhowater &
                 ,CLIQ=cliq,CICE=cice,PSAT=psat &
                 ,RAIN=rainnc ,RAINNCV=rainncv &
                 ,SNOW=snownc ,SNOWNCV=snowncv &
                 ,SR=sr &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1196,&
'arguments not present for calling wsm3' )
             ENDIF
        CASE (WSM5SCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling wsm5' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
             CALL wsm5( &
                  TH=th &
                 ,Q=qv_curr &
                 ,QC=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,DEN=rho,PII=pi_phy,P=p,DELZ=dz8w &
                 ,DELT=dt,G=g,CPD=cp,CPV=cpv &
                 ,RD=r_d,RV=r_v,T0C=svpt0 &
                 ,EP1=ep_1, EP2=ep_2, QMIN=epsilon &
                 ,XLS=xls, XLV0=xlv, XLF0=xlf &
                 ,DEN0=rhoair0, DENR=rhowater &
                 ,CLIQ=cliq,CICE=cice,PSAT=psat &
                 ,RAIN=rainnc ,RAINNCV=rainncv &
                 ,SNOW=snownc ,SNOWNCV=snowncv &
                 ,SR=sr &
                 ,REFL_10CM=refl_10cm &
                 ,diagflag=diagflag &
                 ,do_radar_ref=do_radar_ref &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1231,&
'arguments not present for calling wsm5' )
             ENDIF
        CASE (WSM6SCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling wsm6' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
             CALL wsm6( &
                  TH=th &
                 ,Q=qv_curr &
                 ,QC=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,QG=qg_curr &
                 ,DEN=rho,PII=pi_phy,P=p,DELZ=dz8w &
                 ,DELT=dt,G=g,CPD=cp,CPV=cpv &
                 ,RD=r_d,RV=r_v,T0C=svpt0 &
                 ,EP1=ep_1, EP2=ep_2, QMIN=epsilon &
                 ,XLS=xls, XLV0=xlv, XLF0=xlf &
                 ,DEN0=rhoair0, DENR=rhowater &
                 ,CLIQ=cliq,CICE=cice,PSAT=psat &
                 ,RAIN=rainnc ,RAINNCV=rainncv &
                 ,SNOW=snownc ,SNOWNCV=snowncv &
                 ,SR=sr &
                 ,REFL_10CM=refl_10cm &
                 ,diagflag=diagflag &
                 ,do_radar_ref=do_radar_ref &
                 ,GRAUPEL=graupelnc ,GRAUPELNCV=graupelncv &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1268,&
'arguments not present for calling wsm6' )
             ENDIF
        CASE (WDM5SCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling wdm5' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. PRESENT( QNN_CURR ) .AND. &
                  PRESENT ( QNC_CURR ) .AND. PRESENT( QNR_CURR ).AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
             CALL wdm5( &
                  TH=th &
                 ,Q=qv_curr &
                 ,QC=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,NN=qnn_curr &
                 ,NC=qnc_curr &
                 ,NR=qnr_curr &
                 ,DEN=rho,PII=pi_phy,P=p,DELZ=dz8w &
                 ,DELT=dt,G=g,CPD=cp,CPV=cpv,CCN0=n_ccn0 &
                 ,RD=r_d,RV=r_v,T0C=svpt0 &
                 ,EP1=ep_1, EP2=ep_2, QMIN=epsilon &
                 ,XLS=xls, XLV0=xlv, XLF0=xlf &
                 ,DEN0=rhoair0, DENR=rhowater &
                 ,CLIQ=cliq,CICE=cice,PSAT=psat &
                 ,RAIN=rainnc ,RAINNCV=rainncv &
                 ,SNOW=snownc ,SNOWNCV=snowncv &
                 ,SR=sr &
                 ,REFL_10CM=refl_10cm &
                 ,diagflag=diagflag &
                 ,do_radar_ref=do_radar_ref &
                 ,ITIMESTEP=itimestep &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1308,&
'arguments not present for calling wdm5')
             ENDIF
       CASE (WDM6SCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling wdm6' )
             IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
                  PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
                  PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. &
                  PRESENT( QNN_CURR ) .AND. PRESENT ( QNC_CURR ) .AND. &
                  PRESENT( QNR_CURR ).AND. &
                 PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
             CALL wdm6( &
                  TH=th &
                 ,Q=qv_curr &
                 ,QC=qc_curr &
                 ,QR=qr_curr &
                 ,QI=qi_curr &
                 ,QS=qs_curr &
                 ,QG=qg_curr &
                 ,NN=qnn_curr &
                 ,NC=qnc_curr &
                 ,NR=qnr_curr &
                 ,DEN=rho,PII=pi_phy,P=p,DELZ=dz8w &
                 ,DELT=dt,G=g,CPD=cp,CPV=cpv,CCN0=n_ccn0 &
                 ,RD=r_d,RV=r_v,T0C=svpt0 &
                 ,EP1=ep_1, EP2=ep_2, QMIN=epsilon &
                 ,XLS=xls, XLV0=xlv, XLF0=xlf &
                 ,DEN0=rhoair0, DENR=rhowater &
                 ,CLIQ=cliq,CICE=cice,PSAT=psat &
                 ,RAIN=rainnc ,RAINNCV=rainncv &
                 ,SNOW=snownc ,SNOWNCV=snowncv &
                 ,SR=sr &
                 ,REFL_10CM=refl_10cm &
                 ,diagflag=diagflag &
                 ,do_radar_ref=do_radar_ref &
                 ,GRAUPEL=graupelnc ,GRAUPELNCV=graupelncv &
                 ,ITIMESTEP=itimestep &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                                                                    )
             ELSE
               CALL wrf_error_fatal3("<stdin>",1351,&
'arguments not present for calling wdm6')
             ENDIF
        CASE (ETAMPNEW)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling etampnew')
             IF ( PRESENT( qv_curr ) .AND. PRESENT( qt_curr ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                  PRESENT( mp_restart_state ) .AND. &
                  PRESENT( tbpvs_state ) .AND. &
                  PRESENT( tbpvs0_state ) ) THEN
               CALL ETAMP_NEW( &
                  ITIMESTEP=itimestep,DT=dt,DX=dx,DY=dy &
                 ,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,PI_PHY=pi_phy,TH_PHY=th &
                 ,QV=qv_curr &
                 ,QC=qc_curr &
                 ,QS=qs_curr &
                 ,QR=qr_curr &
                 ,QT=qt_curr &
                 ,LOWLYR=LOWLYR,SR=SR &
                 ,F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY &
                 ,F_RIMEF_PHY=F_RIMEF_PHY &
                 ,RAINNC=rainnc,RAINNCV=rainncv &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                 ,MP_RESTART_STATE=mp_restart_state &
                 ,TBPVS_STATE=tbpvs_state,TBPVS0_STATE=tbpvs0_state &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1381,&
'arguments not present for calling etampnew' )
             ENDIF
          CASE (CAMMGMPSCHEME)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling CAMMGMPSCHEME')
             IF ( PRESENT( z ) .AND. PRESENT( ht ) .AND. &
                  PRESENT( qs_curr ) .AND. &
                  PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
                  PRESENT( qi_curr ) .AND. PRESENT( f_qc ) .AND. &
                  PRESENT( qr_curr ) .AND. PRESENT( qndrop_curr ) .AND. &
                  PRESENT( f_qi ) .AND. PRESENT( qnc_curr ) .AND. &
                  PRESENT( RAINNCV ) .AND. PRESENT( SNOWNCV ) .AND. &
                  PRESENT( qns_curr ) .AND. PRESENT( qnr_curr ) .AND. &
                  PRESENT( chem ) .AND. PRESENT(dgnum4D ) .AND. &
                  PRESENT( dgnumwet4D ) .AND. &
                  PRESENT( qni_curr ) .AND. PRESENT( RAINNC ) ) THEN
                qv_b4mp(its:ite,kts:kte,jts:jte) = qv_curr(its:ite,kts:kte,jts:jte)
                qc_b4mp(its:ite,kts:kte,jts:jte) = qc_curr(its:ite,kts:kte,jts:jte)
                qi_b4mp(its:ite,kts:kte,jts:jte) = qi_curr(its:ite,kts:kte,jts:jte)
                qs_b4mp(its:ite,kts:kte,jts:jte) = qs_curr(its:ite,kts:kte,jts:jte)
                CALL CAMMGMP(ITIMESTEP=itimestep,DT=dt,P8W=p8w_hyd,P_HYD=p_hyd &
                     ,T_PHY=t_phy,PI_PHY=pi_phy,Z_AT_W=z_at_w,QFX=qfx &
                     ,TKE_PBL=tke_pbl,TURBTYPE3D=turbtype3d,SMAW3D=smaw3d &
                     ,DLF3D=dlf,DLF2_3D=dlf2,RLIQ2D=rliq,Z_SEA_LEVEL=z &
                     ,KVH3D=exch_h,HT=ht,ALT=alt,ACCUM_MODE=accum_mode &
                     ,AITKEN_MODE=aitken_mode,COARSE_MODE=coarse_mode &
                     ,ICWMRSH3D=icwmrsh3d,ICWMRDP3D=icwmrdp3d,SHFRC3D=shfrc3d &
                     ,CMFMC3D=cmfmc3d,CMFMC2_3D=cmfmc2_3d &
                     ,CONFIG_FLAGS=config_flags,F_ICE_PHY=f_ice_phy &
                     ,F_RAIN_PHY=f_rain_phy &
                     ,DGNUM4D=dgnum4D,DGNUMWET4D=dgnumwet4D &
                     ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                     ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                     ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                     ,TH=th,CLDFRA_OLD_MP=cldfra_old_mp,CLDFRA_MP=cldfra_mp &
                     ,CLDFRA_MP_ALL=cldfra_mp_all,lradius=lradius,iradius=iradius &
                     ,CLDFRAI=cldfrai,CLDFRAL=cldfral &
                     ,CLDFRA_CONV=cldfra_conv,WSEDL3D=wsedl3d &
                     ,RAINNC=rainnc,RAINNCV=rainncv,SNOWNC=snownc,SNOWNCV=snowncv &
                     ,SR=sr,QV_CURR=qv_curr,QC_CURR=qc_curr,QI_CURR=qi_curr &
                     ,QS_CURR=qs_curr,QR_CURR=qr_curr,NC3D=qnc_curr &
                     ,NI3D=qni_curr,NS3D=qns_curr,NR3D=qnr_curr,QNDROP=qndrop_curr&
                     ,RH_OLD_MP=rh_old_mp,LCD_OLD_MP=lcd_old_mp &
                     ,CHEM=chem &
                     ,QME3D=qme3d,PRAIN3D=prain3d,NEVAPR3D=nevapr3d &
                     ,RATE1ORD_CW2PR_ST3D=rate1ord_cw2pr_st3d &
                     ,XLAND=XLAND,SNOWH=SNOWH &
                     )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1432,&
'arguments not present for calling CAMMGMP SCHEME' )
             ENDIF
        CASE (ETAMPOLD)
             CALL wrf_debug ( 100 , 'microphysics_driver: calling etampold')
             IF ( PRESENT( qv_curr ) .AND. PRESENT( qt_curr ) .AND. &
                  PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) .AND. &
                  PRESENT( mp_restart_state ) .AND. &
                  PRESENT( tbpvs_state ) .AND. &
                  PRESENT( tbpvs0_state ) ) THEN
               CALL ETAMP_OLD( &
                  ITIMESTEP=itimestep,DT=dt,DX=dx,DY=dy &
                 ,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,PI_PHY=pi_phy,TH_PHY=th &
                 ,QV=qv_curr &
                 ,QC=qc_curr &
                 ,QS=qs_curr &
                 ,QR=qr_curr &
                 ,QT=qt_curr &
                 ,LOWLYR=LOWLYR,SR=SR &
                 ,F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY &
                 ,F_RIMEF_PHY=F_RIMEF_PHY &
                 ,RAINNC=rainnc,RAINNCV=rainncv &
                 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
                 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
                 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
                 ,MP_RESTART_STATE=mp_restart_state &
                 ,TBPVS_STATE=tbpvs_state,TBPVS0_STATE=tbpvs0_state &
                                                                    )
             ELSE
                CALL wrf_error_fatal3("<stdin>",1463,&
'arguments not present for calling etampold' )
             ENDIF
      CASE DEFAULT
         WRITE( wrf_err_message , * ) 'The microphysics option does not exist: mp_physics = ', mp_physics
         CALL wrf_error_fatal3("<stdin>",1470,&
wrf_err_message )
      END SELECT micro_select
   ENDDO
   !$OMP END PARALLEL DO
   CALL wrf_debug ( 200 , 'microphysics_driver: returning from' )
   RETURN
   END SUBROUTINE microphysics_driver
END MODULE module_microphysics_driver


text:
c*****************************************************************************
c     Subroutine TIECHV
C
C     CHECK FOR 'TIEPOINTS' OPTION ERRORS AND
C     OPTIONALLY CORRECT DISTORTION
C
C     ARGUMENTS:
C       IND   - INDICATOR ON RETURN (0 IF OK)
C       ITYPE - 0 = NO DISTORTION CORRECTION
C               1 = CORRECT FOR DISTORTION (USE GEOMA ROUTINE GEMDST)
C       PAR   - PARAMETER BUFFER
C               IN ORDER OF LINES,SAMPLES,LAT,LONG,......
C       NPOINT- NUMBER OF POINTS
C       CONV  - GEOMA parameters to be used by TRITRA
C               for geometric correction of image space
C               tiepoints.
C               The GEOMA parameters are in the format:
C               conv(1)='NAH '
C               conv(2)='    '
C               conv(3)=value (value of 'NAH',number of horizontal areas)
C               conv(4)='NAV '
C               conv(5)='    '
C               conv(6)=value (value of 'NAV',number of vertical areas)
C               conv(7)='TIEP'
C               conv(8)='    '
C               conv(9)=   beginning of GEOMA tiepoint parameters
C
C     HISTORY:
C      13 JUL 93  T. L. TRUONG  PORTED TO UNIX
C      16 JAN 78    ...JJL...   INITIAL RELEASE
c*****************************************************************************
C
      SUBROUTINE TIECHV(IND,ITYPE,PAR,NPOINT,CONV)
      CHARACTER*132 MSG
      REAL*4 PAR(*)
      INTEGER*4 CONV(*)
      IND=0
      J=NPOINT*4
      IF(NPOINT.LT.3) GO TO 140
      IF(NPOINT.GT.20) GO TO 130
      WRITE (MSG,9900) NPOINT
9900  FORMAT (' ',I2,' TIEPOINTS SUPPLIED BY USER')
      CALL XVMESSAGE(MSG(2:30),' ')
131   DO 101 I=1,J,4
      IF(ABS(PAR(I+2)).LT.90.005) GO TO 120
      IND=-1
      CALL XVMESSAGE(' TIEPTS LATITUDE ERROR',' ')
      RETURN
120   IF(ABS(PAR(I+3)).LT.360.005) GO TO 110
      IND=-1
      CALL XVMESSAGE(' TIEPTS LONGITUDE ERROR',' ')
      RETURN
110   IF(ITYPE.EQ.0) GO TO 101
C***************************************
C     CORRECT LINE SAMPLE FOR DISTORTION
C***************************************
      CALL TRITRA(IND,CONV(9),CONV(3)+1,CONV(6)+1,PAR(I),PAR(I+1),1)
101   CONTINUE
      RETURN
140   CALL XVMESSAGE(' ***TOO FEW TIEPOINTS',' ')
      IND=-1
      RETURN
130   CALL XVMESSAGE(' ***TOO MANY TIEPOINTS SUPPLIED--FIRST 20 '//
     +			'USED, REST IGNORED',' ')
      NPOINT=20
      J=NPOINT*4
      GO TO 131
      END

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
| vim: ft=reva
|
| process.f
|
| This program demonstrates usage of the os/process library.
| Ensure that reva\bin is in your library search path!

needs os/process
needs os/fs
with~ ~process

: reva-path appname zcount ;

create mystr 100 allot
quote %
  ." hello from the started process!" cr
  bye
% mystr lplace
: the-string mystr lcount ;

: start-it ( -- process )
." starting process " reva-path type cr
  reva-path create-process0
  dup RSS_PROCESS_ALL start-process drop

." the stdin  fd of the child process: " dup get-to-fd   . cr
." the stdout fd of the child process: " dup get-from-fd . cr
;

: check-alive ( process -- process )
  dup is-alive
  ." status of the process: " . cr
;

: send-it ( process -- process )
 dup the-string send-to-process ." bytes sent: " . cr
 dup close-to drop
;

create mytmp 1024 allot
: read-it ( process -- process )
  dup mytmp 1024 receive-from-process
  ." bytes read: " dup . cr
  dup 0if drop else
    mytmp swap type cr
    read-it
  ;then
;

: stop-it ( process -- )
  RSS_PROCESS_WAIT release-process drop
;

: go
  start-it
  send-it
  read-it
  stop-it
;

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
*----------------------------------------------------------------------*
      subroutine spinsym_traf(mode,ndim,mat,map,nsing,sing,trip,half)
*----------------------------------------------------------------------*
*     does a unitary (half-)transformation using a spin-symmetrization
*     map as input for the unitary matrix U.
*     the result is split into a "singlet" and a "triplet" block;
*     in case of the back-transformation (mode=2), these two blocks
*     are rejoined.
*
*     matthias, feb 2010
*----------------------------------------------------------------------*
      implicit none

      include 'stdunit.h'

      integer, parameter ::
     &     ntest = 00
      integer, intent(in) ::
     &     mode, ndim, nsing, map(ndim)
      real(8), intent(inout), target ::
     &     mat(ndim,ndim), sing(nsing,nsing),
     &     trip(ndim-nsing,ndim-nsing)
      logical, intent(in) ::
     &     half
      real(8) ::
     &     sig, fac
      real(8), pointer ::
     &     scr(:,:)

      integer ::
     &     icol, iline, ntrip, ising, itrip
      fac = 1d0/sqrt(2d0)

      ntrip = ndim - nsing
      if (mode.lt.1.or.mode.gt.2.or.ntrip.lt.0
     &    .or.(half.and.mode.eq.1))
     &     call quit(1,'spinsym_traf','are you kidding me?')

      if (ntest.ge.100) then
        write(lulog,'(x,a)') '-------------------------------'
        write(lulog,'(x,a,i1,a)') 'spinsym_traf at work (mode = ',
     &                            mode,')'
        write(lulog,'(x,a)') '-------------------------------'
      end if

      if (mode.eq.1) then
        if (ntest.ge.100) then
          write(lulog,*) 'input matrix:'
          call wrtmat2(mat,ndim,ndim,ndim,ndim)
        end if
        allocate(scr(ndim,ndim))
        scr = 0d0
        sing = 0d0
        trip = 0d0

        ! first mat*U
        ising = 0
        itrip = 0
        do icol = 1, ndim
          if (abs(map(icol)).gt.ndim) call quit(1,'spinsym_traf',
     &           'map must stay within matrix dimension!')
          if (abs(map(icol)).eq.icol) then
            if (map(icol).gt.0) then
              ising = ising + 1
              scr(1:ndim,ising) = mat(1:ndim,icol)
            else
              itrip = itrip + 1
              scr(1:ndim,nsing+itrip) = mat(1:ndim,icol)
            end if
          else if (abs(map(icol)).gt.icol) then
            sig = dble(sign(1,map(icol)))
            ising = ising + 1
            scr(1:ndim,ising)
     &       = fac*mat(1:ndim,icol)
     &       + sig*fac*mat(1:ndim,abs(map(icol)))
          else
            sig = - dble(sign(1,map(icol)))
            itrip = itrip + 1
            scr(1:ndim,nsing+itrip)
     &       = fac*mat(1:ndim,abs(map(icol)))
     &       + sig*fac*mat(1:ndim,icol)
          end if
        end do
        if (ising.ne.nsing.or.ising+itrip.ne.ndim) call quit(1,
     &      'spinsym_traf','error in singlet triplet splitting')

        ! now U^+*(mat*U) --> (sing|trip)
        ising = 0
        itrip = 0
        do iline = 1, ndim
          if (abs(map(iline)).eq.iline) then
            if (map(iline).gt.0) then
              ising = ising + 1
              sing(ising,1:nsing) = scr(iline,1:nsing)
            else
              itrip = itrip + 1
              trip(itrip,1:ntrip) = scr(iline,nsing+1:ndim)
            end if
          else if (abs(map(iline)).gt.iline) then
            sig = dble(sign(1,map(iline)))
            ising = ising + 1
            sing(ising,1:nsing)
     &       = fac*scr(iline,1:nsing)
     &       + sig*fac*scr(abs(map(iline)),1:nsing)
          else
            sig = - dble(sign(1,map(iline)))
            itrip = itrip + 1
            trip(itrip,1:ntrip)
     &       = fac*scr(abs(map(iline)),nsing+1:ndim)
     &       + sig*fac*scr(iline,nsing+1:ndim)
          end if
        end do
        deallocate(scr)
        if (ntest.ge.100) then
          write(lulog,*) 'output "singlet" block:'
          call wrtmat2(sing,nsing,nsing,nsing,nsing)
          write(lulog,*) 'output "triplet" block:'
          call wrtmat2(trip,ntrip,ntrip,ntrip,ntrip)
        end if

      else if (mode.eq.2) then
        if (ntest.ge.100) then
          write(lulog,*) 'input "singlet" block:'
          call wrtmat2(sing,nsing,nsing,nsing,nsing)
          write(lulog,*) 'input "triplet" block:'
          call wrtmat2(trip,ntrip,ntrip,ntrip,ntrip)
        end if

        if (half) then
          scr => mat
        else
          allocate(scr(ndim,ndim))
        end if

        ! first U*(sing|trip)
        scr = 0d0
        ising = 0
        itrip = 0
        do iline = 1, ndim
          if (abs(map(iline)).gt.ndim) call quit(1,'spinsym_traf',
     &           'map must stay within matrix dimension!')
          if (abs(map(iline)).eq.iline) then
            if (map(iline).gt.0) then
              ising = ising + 1
              scr(iline,1:nsing) = sing(ising,1:nsing)
            else
              itrip = itrip + 1
              scr(iline,nsing+1:ndim) = trip(itrip,1:ntrip)
            end if
          else if (abs(map(iline)).gt.iline) then
            sig = dble(sign(1,map(iline)))
            ising = ising + 1
            scr(iline,1:nsing) = fac*sing(ising,1:nsing)
            scr(abs(map(iline)),1:nsing) = sig*fac*sing(ising,1:nsing)
          else
            sig = - dble(sign(1,map(iline)))
            itrip = itrip + 1
            scr(abs(map(iline)),nsing+1:ndim) = fac*trip(itrip,1:ntrip)
            scr(iline,nsing+1:ndim) = sig*fac*trip(itrip,1:ntrip)
          end if
        end do
        if (ising.ne.nsing.or.ising+itrip.ne.ndim) call quit(1,
     &      'spinsym_traf','error in singlet triplet splitting')

        ! if full trafo is required, do [U*(sing|trip)]*U^+
        if (.not.half) then
          mat = 0d0
          ising = 0
          itrip = 0
          do icol = 1, ndim
            if (abs(map(icol)).eq.icol) then
              if (map(icol).gt.0) then
                ising = ising + 1
                mat(1:ndim,icol) = scr(1:ndim,ising)
              else
                itrip = itrip + 1
                mat(1:ndim,icol) = scr(1:ndim,nsing+itrip)
              end if
            else if (abs(map(icol)).gt.icol) then
              sig = dble(sign(1,map(icol)))
              ising = ising + 1
              mat(1:ndim,icol) = fac*scr(1:ndim,ising)
              mat(1:ndim,abs(map(icol))) = sig*fac*scr(1:ndim,ising)
            else
              sig = - dble(sign(1,map(icol)))
              itrip = itrip + 1
              mat(1:ndim,abs(map(icol))) = mat(1:ndim,abs(map(icol)))
     &                           + fac*scr(1:ndim,nsing+itrip)
              mat(1:ndim,icol) = mat(1:ndim,icol)
     &                           + sig*fac*scr(1:ndim,nsing+itrip)
            end if
          end do
          deallocate(scr)
        end if
c dbg   can be used to set 1
c        mat(1:ndim,1:ndim) = 0d0
c        do iline = 1, ndim
c          mat(iline,iline) = 1d0
c        end do
c dbgend

        if (ntest.ge.100) then
          write(lulog,*) 'output matrix:'
          call wrtmat2(mat,ndim,ndim,ndim,ndim)
        end if
      end if

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
c=====================================================================72
      subroutine write_tec(cfile,ct,cx,cy,cv, nx,ny,sof)
      implicit none
      character*(*) cfile,ct,cx,cy,cv
      integer nx, ny
      real sof(0:nx, 0:ny)
      integer i1, i2, i, j

      write (6,*) "at top of write_tec"

c      cv = "PSD     "
c      cx = "x       "
c      cy = "y       "
c      ct = "Is MINE "
c      cfile = "a.pdt   "

      write (6,*) trim(ct)
      write (6,*) trim(cx)
      write (6,*) trim(cy)
      write (6,*) trim(cv)
      write (6,*) trim(cfile)

      open(22,file=cfile,form="formatted")
900   format(a,a,a,a,a)
      write (22,900) " TITLE = ", '"', trim(ct), '"'
      write (22,900) "     VARIABLES=", '"', trim(cx), '"', ','
      write (22,900)                    '"', trim(cy), '"', ','
      write (22,900)                    '"', trim(cv), '"'

901   format("     ZONE  I=",i4,", J=",i4,", F=BLOCK,")
      write (22,901) nx, ny
      write (22,900) "           T=",'"', "2-D DATA", '"'

902   format(1p7e11.3)
      do j  = 1, ny
      do i1 = 1, nx, 7
        i2 = min(i1+6, nx)
        write (22,902) (sof(i,0),i=i1,i2)
      enddo
      enddo

      do j  = 1, ny
      do i1 = 1, nx, 7
        i2 = min(i1+6, nx)
        write (22,902) (sof(0,j),i=i1,i2)
      enddo
      enddo

      do j  = 1, ny
      do i1 = 1, nx, 7
        i2 = min(i1+6, nx)
        write (22,902) (sof(i,j),i=i1,i2)
      enddo
      enddo
      close(22)

      return
      end
c=====================================================================72

max_stars_repo_name: nannau/astrohubs

max_issues_repo_name: nannau/astrohubs

max_forks_repo_name: nannau/astrohubs


text:
! { dg-do run }
! PR33162 INTRINSIC functions as ACTUAL argument
! Test case adapted from PR by Jerry DeLisle <jvdelisle@gcc.gnu.org>
real function t(x)
  real ::x
  t = x
end function

program p
  implicit none
  intrinsic sin
  procedure(sin):: t
  if (t(1.0) /= 1.0) call abort
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


module DataInputM
  use UtilitiesM
  use DebuggerM

  use NodeM
  use NodePtrM
  use GeometryM
  use SourceM
  use StructuralMaterialM
  use StructuralElementM
  use PressureM
  use Structural3DApplicationM

  use MeshM
  use ModelM

  implicit none

  private
  public :: initFEM3D

  integer(ikind), parameter    :: projectData = 1
  integer(ikind), parameter    :: project = 2
  integer(ikind), parameter    :: functions = 4
  integer(ikind), dimension(8) :: date_time
  integer(ikind), parameter    :: nDof = 3
  integer(ikind)               :: nElem
  integer(ikind)               :: nTetraElem
  integer(ikind)               :: nHexaElem
  integer(ikind)               :: nPoint
  integer(ikind)               :: iPoint
  integer(ikind)               :: nPressure
  integer(ikind)               :: nDirichletX
  integer(ikind)               :: nDirichletY
  integer(ikind)               :: nDirichletZ
  integer(ikind)               :: nMaterial
  integer(ikind)               :: nGauss
  integer(ikind)               :: nConvection
  integer(ikind)               :: isQuadratic
  integer(ikind)               :: nSourceOnPoints
  integer(ikind)               :: nSourceOnVolumes
  integer(ikind)               :: nPointSource
  integer(ikind)               :: nVolumeSource
  character(100)               :: projectName
  character(100)               :: aux
  logical       , parameter    :: verbose = .false.
  logical                      :: isMaterialAsigned = .true.

  interface initFEM3D
     procedure :: initFEM3D
  end interface initFEM3D

contains

  subroutine initFEM3D(structuralAppl)
    implicit none
    type(Structural3DApplicationDT), intent(inout) :: structuralAppl
    print'(A)', 'Initializing Structural3D application'
    call initLog(.true., 'log.dat')
    call debugLog('  Reading project data')
    call readProjectData
    call debugLog('  Reading mesh data')
    call initMesh(structuralAppl)
    call debugLog('  Reading materials properties')
    call initMaterials(structuralAppl)
    call debugLog('  Reading elements')
    call initElements(structuralAppl)
    call debugLog('  Reading point and line Sources')
    call readSources(structuralAppl)
    call debugLog('  Reading Boundary Conditions')
    call readBoundaryConditions(structuralAppl)
    call debugLog('End loading data')
  end subroutine initFEM3D

  subroutine readProjectData
    implicit none
    open(projectData, file = 'projectData.dat')
    read(projectData, '(*(A))') projectName
    close(projectData)
    call debugLog('    Project name: ', trim(projectName))
  end subroutine readProjectData

  subroutine initMesh(structuralAppl)
    implicit none
    type(Structural3DApplicationDT), intent(inout) :: structuralAppl
    integer(ikind) :: i
    integer(ikind) :: nnz
    real(rkind)    :: x, y, z
    open(project, file = trim(projectName)//'.dat')
    do i = 1, 9
       read(project,*)
    end do
    read(project,*)  aux, nElem
    read(project,*)  aux, nPoint
    read(project,*)  aux, isQuadratic
    read(project,*)  aux, nTetraElem
    read(project,*)  aux, nHexaElem
    read(project,*)  aux, nMaterial
    read(project,*)  aux, nGauss
    read(project,*)  aux, nDirichletX
    read(project,*)  aux, nDirichletY
    read(project,*)  aux, nDirichletZ
    read(project,*)  aux, nPressure
    read(project,*)  aux, nSourceOnPoints
    read(project,*)  aux, nSourceOnVolumes
    read(project,*)  aux, nPointSource
    read(project,*)  aux, nVolumeSource

    if(verbose) print'(A,I0)','Number of Elements.............................: ', nElem
    if(verbose) print'(A,I0)','Are Elements Quadratic.........................: ', isQuadratic
    if(verbose) print'(A,I0)','Number of Tetrahedral elements.................: ', nTetraElem
    if(verbose) print'(A,I0)','Number of Hexahedral elements..................: ', nHexaElem
    if(verbose) print'(A,I0)','Number of Nodes................................: ', nPoint
    if(verbose) print'(A,I0)','Number of Dirichlet X conditions...............: ', nDirichletX
    if(verbose) print'(A,I0)','Number of Dirichlet Y conditions...............: ', nDirichletY
    if(verbose) print'(A,I0)','Number of Dirichlet Z conditions...............: ', nDirichletZ
    if(verbose) print'(A,I0)','Number of Pressure conditions..................: ', nPressure
    if(verbose) print'(A,I0)','Number of Loads on points......................: ', nSourceOnPoints
    if(verbose) print'(A,I0)','Number of Loads on volumes.....................: ', nSourceOnVolumes
    if(verbose) print'(A,I0)','Number of points with pointSource..............: ', nPointSource
    if(verbose) print'(A,I0)','Number of volumes with volumeSource............: ', nVolumeSource
    if(verbose) print'(A,I0)','Number of Materials............................: ', nMaterial
    if(verbose) print'(A,I0)','Gauss cuadrature order.........................: ', nGauss

    call debugLog('    Number of Elements.............................: ', nElem)
    call debugLog('    Are Elements Quadratic.........................: ', isQuadratic)
    call debugLog('    Number of Tetrahedral elements.................: ', nTetraElem)
    call debugLog('    Number of Hexahedral elements..................: ', nHexaElem)
    call debugLog('    Number of Nodes................................: ', nPoint)
    call debugLog('    Number of Dirichlet X conditions...............: ', nDirichletX)
    call debugLog('    Number of Dirichlet Y conditions...............: ', nDirichletY)
    call debugLog('    Number of Dirichlet Z conditions...............: ', nDirichletZ)
    call debugLog('    Number of Pressure conditions..................: ', nPressure)
    call debugLog('    Number of Sources on points....................: ', nSourceOnPoints)
    call debugLog('    Number of Sources on surfaces..................: ', nSourceOnVolumes)
    call debugLog('    Number of points with pointSource..............: ', nPointSource)
    call debugLog('    Number of Volumes with volumeSource............: ', nVolumeSource)
    call debugLog('    Number of Materials............................: ', nMaterial)
    call debugLog('    Gauss cuadrature order.........................: ', nGauss)

    if(isQuadratic == 0) then
       nnz = nElem*8*8*3*3
    else if (isQuadratic == 1) then
       nnz = nElem*20*20*3*3
    end if

    structuralAppl = structural3DApplication(                  &
           nNode = nPoint                                      &
         , nElement = nElem                                    &
         , nPressure = nPressure                               &
         , nSource = nSourceOnPoints + nSourceOnVolumes        &
         , nMaterial = nMaterial                               &
         , nGauss = nGauss                                     &
         , nnz = nnz                                           )

    do i = 1, 6
       read(project,*)
    end do
    if(verbose) print'(A)', 'Coordinates:'
    if(verbose) print'(1X,A)', 'NODES      X          Y'
    do i = 1, nPoint
       read(project,*) iPoint, x, y, z
       if(verbose) print'(3X,I0,3X,E10.3,3X,E10.3,3X,E10.3)', iPoint, x, y, z
       structuralAppl%node(iPoint) = node(iPoint, 3, x, y, z)
       call structuralAppl%node(iPoint)%assignDof(1, structuralAppl%model%dof(iPoint*nDof-2))
       call structuralAppl%node(iPoint)%assignDof(2, structuralAppl%model%dof(iPoint*nDof-1))
       call structuralAppl%node(iPoint)%assignDof(3, structuralAppl%model%dof(iPoint*nDof))
       call structuralAppl%model%addNode(iPoint, structuralAppl%node(iPoint))
    end do
  end subroutine initMesh

  subroutine initMaterials(structuralAppl)
    implicit none
    type(Structural3DApplicationDT), intent(inout) :: structuralAppl
    integer(ikind) :: i, iMat
    real(rkind)    :: alpha, E, nu
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(A)', 'Material         alpha        E        nu'
    do i = 1, nMaterial
       read(project,*) iMat, alpha, E, nu
       structuralAppl%material(iMat) = structuralMaterial(E, nu, alpha)
       if(verbose) print'(4X,I0,7X,5(E10.3,3X))', iMat, alpha, E, nu
    end do
  end subroutine initMaterials

  subroutine initElements(structuralAppl)
    type(Structural3DApplicationDT), intent(inout) :: structuralAppl
    type(NodePtrDT), dimension(:), allocatable :: auxNode
    integer(ikind) :: i, j, iElem, iMat, nNode, Conectivities(20)
    character(len=13) :: type
    Conectivities = 0
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(A)', 'Element  |      Type      |  material index  |  nNodes  |  connectivities'
    do i = 1, nElem
       read(project,*) iElem, type, iMat, nNode, (Conectivities(j),j=1,nNode)
       if(verbose) print'(I5,A15,I18,I14,5X,*(I5,X))', iElem, type, iMat, nNode, (Conectivities(j),j=1,nNode)
       allocate(auxNode(nNode))
       do j = 1, nNode
          call auxNode(j)%associate(structuralAppl%node(conectivities(j)))
       end do
       structuralAppl%element(iElem) = structuralElement(iElem, auxNode, structuralAppl%material(iMat))
       call structuralAppl%model%addElement(i, structuralAppl%element(iElem))
       deallocate(auxNode)
    end do
  end subroutine initElements

  subroutine readSources(structuralAppl)
    implicit none
    type(Structural3DApplicationDT), intent(inout) :: structuralAppl
    integer(ikind)                              :: i, countSource, auxInt
    integer(ikind)                              :: iNode, iElem, iSource
    character(150), dimension(3)                :: func
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'nSource'
    if(verbose) print'(A)', 'Source    Function'
    do i = 1, nSourceOnPoints+nSourceOnVolumes
       read(project,*) iSource, func(1), func(2), func(3)
       structuralAppl%source(iSource) = source(3, 3, (/'x', 'y', 'z'/), func)
       if(verbose) print'(I0,5X,30A,30A,30A)', iSource, func(1), func(2), func(3)
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'pointSources'
    if(verbose) print'(A)', 'Node    Load'
    do i = 1, nPointSource
       read(project,*) iNode, iSource
       if(verbose) print'(I0,5X,I0)', iNode, iSource
       call structuralAppl%node(iNode)%assignSource(structuralAppl%source(iSource))
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'volumeSources'
    if(verbose) print'(A)', 'Element   Load'
    do i = 1, nVolumeSource
       read(project,*) iElem, iSource
       if(verbose) print'(I0,5X,I0)', iElem, iSource
       call structuralAppl%element(iElem)%assignSource(structuralAppl%source(iSource))
    end do
  end subroutine readSources

  subroutine readBoundaryConditions(structuralAppl)
    implicit none
    type(Structural3DApplicationDT), intent(inout)  :: structuralAppl
    integer(ikind)                               :: i, j, id, elemID, nPointID
    integer(ikind)                               :: iPoint, conditionCounter
    integer(ikind), dimension(:), allocatable    :: pointID
    real(rkind)                                  :: value
    real(rkind)                                  :: coef, temp
    type(StructuralElementDT)                    :: element
    type(NodePtrDT)  , dimension(:), allocatable :: node
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'Dirichlet X conditions'
    if(verbose) print'(A)', 'Node    Value'
    do i = 1, nDirichletX
       read(Project,*) id, value
       if(verbose) print'(I0,5X,E10.3)', id, value
       call structuralAppl%node(id)%fixDof(1, value)
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'Dirichlet Y conditions'
    if(verbose) print'(A)', 'Node    Value'
    do i = 1, nDirichletY
       read(Project,*) id, value
       if(verbose) print'(I0,5X,E10.3)', id, value
       call structuralAppl%node(id)%fixDof(2, value)
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'Dirichlet Z conditions'
    if(verbose) print'(A)', 'Node    Value'
    do i = 1, nDirichletZ
       read(Project,*) id, value
       if(verbose) print'(I0,5X,E10.3)', id, value
       call structuralAppl%node(id)%fixDof(3, value)
    end do
    do i = 1, 7
       read(project,*)
    end do
    allocate(pointID(20))
    conditionCounter = 0
    if(verbose) print'(/,A)', 'Pressure On Surfaces conditions'
    if(verbose) print'(A)', 'Elem    Nodes     Value'
    do i = 1, nPressure
       conditionCounter = conditionCounter + 1
       read(Project,*) elemID, nPointID, (pointID(j),j=1,nPointID), value
       if(verbose) print*, elemID, (pointID(j),j=1,nPointID), value
       element = structuralAppl%element(elemID)
       allocate(node(nPointID))
       do j = 1, nPointID
          node(j) = element%node(pointID(j))
       end do
       structuralAppl%pressure(i) = &
            pressure(i, pointID, value, node, element%geometry, element%material)
       call structuralAppl%model%addCondition(conditionCounter, structuralAppl%pressure(i))
       deallocate(node)
    end do
    close(project)
  end subroutine readBoundaryConditions

end module DataInputM


text:
      SUBROUTINE CROSSV(VI1,VI2,VO)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION VI1(3),VI2(3),VO(3)
C
      VO(1)=VI1(2)*VI2(3)-VI1(3)*VI2(2)
      VO(2)=VI1(3)*VI2(1)-VI1(1)*VI2(3)
      VO(3)=VI1(1)*VI2(2)-VI1(2)*VI2(1)
C
      RETURN
C
      END

ext: for


text:
!-------------------------------------------------------------------------------
! implementation of a convolutional/deconvolutional layer type (ConvLayer)
!-------------------------------------------------------------------------------

module conv_layer_definitions
use net_helper_procedures
use pool_layer_definitions
implicit none

!===============================================================================
!===============================================================================
! procedures with 4D array input require (rows, columns, channels, batches) form
!===============================================================================
!===============================================================================

!===============================================================================
! types
!===============================================================================

! represents a Convolutional Layer in a neural network
!
! the ConvLayer keeps track of the kernels that pass over (cross-correlate with)
! the input channels, and the biases, weighted inputs, activations, and deltas
! (errors) that are used on the result of passing the kernels over the input
! channels
!
! each kernel is shared between images in a batch, and each produces an output
! channel. kernels have 4D matrix dimensions: (rows, columns, channels, batches)
!
! the implementation here uses 'tied' biases, where there is 1 bias per channel;
! bias value repeated into higher dimension array so we can add array later
! rather than using nested loops
!
! weighted inputs, activations, and deltas match the output channel produced by
! passing a kernel over an input, and are passed between ConvLayers to preserve
! batch information. each have 4D matrices with dimensions:
! (rows, columns, channels, batches)
!
! an input image can have any number of channels (1 for grayscale, 3 for RGB);
! the third dimension for all of the matrices is the channels, so no special
! processing must be done; an input must just be specified with 3 channels
!
! ConvLayers are connected to each other as a doubly-linked list, and each
! ConvLayer can have a PoolLayer that it feeds into before reaching the next
! ConvLayer
type :: ConvLayer
    integer                   :: in_dims(3), k_dims(2), stride(2), &
                                 out_rows, out_cols, out_channels, out_count, &
                                 batch_size
    class(ConvLayer), pointer :: prev_layer, next_layer
    class(PoolLayer), pointer :: next_pool
    character(len=20)         :: pad, activ
    ! dimension order: (rows, columns, channels, batches or kernel items)
    real(kind=8), allocatable :: k(:,:,:,:), & ! kernels
                                 b(:,:,:),   & ! biases
                                 z(:,:,:,:), & ! weighted inputs
                                 a(:,:,:,:), & ! activations
                                 d(:,:,:,:), & ! deltas (errors)
                                 drop(:,:,:,:) ! dropout inputs (0=drop, 1=keep)
    real                      :: drop_rate     ! % of input nodes to drop
contains
    ! procedures that traverse through linked list of ConvLayers
    procedure, pass           :: conv_init, conv_add_pool_layer, &
                                 conv_forw_prop, conv_back_prop, conv_update, &
                                 conv_dropout_rand
end type
contains

!===============================================================================
! constructors / destructors
!===============================================================================

!-------------------------------------------------------------------------------
! constructs a new ConvLayer;
!
! 'same' padding maintains dimensions;
! 'valid' padding reduces dimensions;
! 'full' padding applies opposite operations to 'valid'; upscaling
!-------------------------------------------------------------------------------
! input_dims:  (integer(3)) (height, width, channels) of one input
! kernels:     (integer) kernels in this layer
! kernel_dims: (integer(2)) (height, width) of kernel
! stride:      (integer(2)) size of kernel moves in (y, x) directions
! activ:       (characters) activation function
! padding:     (characters) padding type
! drop_rate:   (real) % of input nodes to dropout
!-------------------------------------------------------------------------------
! returns ::   (ConvLayer pointer) new ConvLayer
!-------------------------------------------------------------------------------
function create_conv_layer(input_dims, kernels, kernel_dims, stride, &
                           activ, padding, drop_rate)
    class(ConvLayer), pointer :: create_conv_layer
    integer, intent(in)       :: input_dims(3), kernels, kernel_dims(2), &
                                 stride(2)
    character(*), intent(in)  :: activ, padding
    real, intent(in)          :: drop_rate
    integer                   :: pad_rows, final_rows, pad_cols, final_cols, &
                                 exp_rows, exp_cols, stride_rows, stride_cols

    if (.not. (padding == 'valid' .or. &
               padding == 'same' .or. &
               padding == 'full')) then
        print *, '---------------------------------------------'
        print *, '(conv_layer_definitions :: create_conv_layer)'
        print *, 'supported: valid, same, full'
        print *, '---------------------------------------------'
        stop -1
    end if

    allocate(create_conv_layer)
    create_conv_layer%in_dims    =  input_dims
    create_conv_layer%k_dims     =  kernel_dims
    create_conv_layer%stride     =  stride
    create_conv_layer%activ      =  activ
    create_conv_layer%pad        =  padding
    create_conv_layer%drop_rate  =  drop_rate
    create_conv_layer%prev_layer => null()
    create_conv_layer%next_layer => null()
    create_conv_layer%next_pool  => null()

    if (padding == 'full') then
        ! for full padding, we first expand the input before kernel
        ! to upscale, then kernel is only applied with strides of 1
        exp_rows    = expand_calc(input_dims(1), stride(1))
        exp_cols    = expand_calc(input_dims(2), stride(2))
        stride_rows = 1
        stride_cols = 1
    else
        ! not full padding; fill in with normal values
        exp_rows    = input_dims(1)
        exp_cols    = input_dims(2)
        stride_rows = stride(1)
        stride_cols = stride(2)
    end if

    ! determine output shape from input pad then applying kernel
    pad_rows   = pad_calc(exp_rows, kernel_dims(1), stride_rows, padding)
    pad_cols   = pad_calc(exp_cols, kernel_dims(2), stride_cols, padding)
    final_rows = res_calc(exp_rows+pad_rows, kernel_dims(1), stride_rows)
    final_cols = res_calc(exp_cols+pad_cols, kernel_dims(2), stride_cols)

    create_conv_layer%out_rows     = final_rows
    create_conv_layer%out_cols     = final_cols
    create_conv_layer%out_channels = kernels

    ! outward facing node count
    create_conv_layer%out_count = create_conv_layer%out_rows * &
                                  create_conv_layer%out_cols * &
                                  create_conv_layer%out_channels
end function

!-------------------------------------------------------------------------------
! deallocate a ConvLayer
!-------------------------------------------------------------------------------
! l:        (ConvLayer pointer)
!-------------------------------------------------------------------------------
! alters :: l is deallocated
!-------------------------------------------------------------------------------
subroutine deallocate_conv_layer(l)
    class(ConvLayer), pointer :: l

    if (associated(l%next_pool)) then
        call deallocate_pool_layer(l%next_pool) ! deallocate trailing pool layer
    end if

    deallocate(l%k, l%b, l%z, l%a, l%d)
    deallocate(l)
end subroutine

!===============================================================================
! ConvLayer procedures
!===============================================================================

!-------------------------------------------------------------------------------
! add a PoolLayer for this ConvLayer to feed into
!-------------------------------------------------------------------------------
! this:        (ConvLayer - implicitly passed)
! kernel_dims: (integer(2)) (height, width) of the pool kernel
! stride:      (integer(2)) size of pool kernel moves in (y, x) directions
! pool:        (characters) pool type
! padding:     (characters) padding type
!-------------------------------------------------------------------------------
! alters ::    this ConvLayer now has a next_pool PoolLayer
!-------------------------------------------------------------------------------
subroutine conv_add_pool_layer(this, kernel_dims, stride, pool, padding)
    class(ConvLayer)         :: this
    integer, intent(in)      :: kernel_dims(2), stride(2)
    character(*), intent(in) :: pool, padding

    ! feed this ConvLayers activations (output) into PoolLayer
    this%next_pool => create_pool_layer([this%out_rows, this%out_cols, &
                                         this%out_channels], &
                                        kernel_dims, stride, pool, padding)
end subroutine

!-------------------------------------------------------------------------------
! helper subroutine to initialize the various matrices of ConvLayers
!-------------------------------------------------------------------------------
! this:       (ConvLayer - implicitly passed)
! batch_size: (integer) examples to process before back prop
!-------------------------------------------------------------------------------
! alters ::   this ConvLayer's matrices are allocated and become usable
!-------------------------------------------------------------------------------
subroutine conv_init(this, batch_size)
    class(ConvLayer)    :: this
    integer, intent(in) :: batch_size

    this%batch_size = batch_size

    allocate(this%k(this%k_dims(1), this%k_dims(2), this%in_dims(3), &
                    this%out_channels), &
             this%b(this%out_rows, this%out_cols, this%out_channels), &

             this%z(this%out_rows, this%out_cols, this%out_channels, &
                    batch_size), &
             this%a(this%out_rows, this%out_cols, this%out_channels, &
                    batch_size), &
             this%d(this%out_rows, this%out_cols, this%out_channels, &
                    batch_size), &

             this%drop(this%in_dims(1), this%in_dims(2), this%in_dims(3), &
                       batch_size))

    ! initialize kernels and biases
    call random_number(this%k)   ! random weights, uniform range: [0,1)
    this%k = (this%k - 0.5) / 10 ! shift and scale to [-0.05, 0.05)
    this%b = 0

    ! init next PoolLayer
    if (associated(this%next_pool)) then
        call this%next_pool%pool_init(batch_size)
    end if

    ! traverse next ConvLayers
    if (associated(this%next_layer)) then
        call this%next_layer%conv_init(batch_size)
    end if
end subroutine

!-------------------------------------------------------------------------------
! helper subroutine to random initialize/overwrite a ConvLayers dropout array
!-------------------------------------------------------------------------------
! this:     (ConvLayer - implicitly passed)
!-------------------------------------------------------------------------------
! alters :: this ConvLayer's drop array has randomized values
!-------------------------------------------------------------------------------
subroutine conv_dropout_rand(this)
    class(ConvLayer)   :: this

    ! initialize dropout layer; 0=drop, 1=keep (we will be multiplying)
    if (this%drop_rate > 0) then
        call random_number(this%drop) ! random weights, uniform range: [0,1)

        ! this%drop_rate are dropped, so 1 - this%drop_rate are kept;
        ! think of keep = 1 - this%drop_rate; add keep, then cast to int to
        ! truncate kept values to 1 and dropped values to 0;
        !
        ! shift to [keep, 1+keep), so after truncating, so
        ! values in top keep proportion are kept at 1
        this%drop = int(this%drop + 1 - this%drop_rate)
    end if
end subroutine

!-------------------------------------------------------------------------------
! helper subroutine to forward propagate through ConvLayers
!-------------------------------------------------------------------------------
! this:     (ConvLayer - implicitly passed)
! input0:   (real(:,:,:,:)) input batch to forward propagate
! is_train: (logical) in training iteration
!-------------------------------------------------------------------------------
! alters :: this ConvLayer's z and a are calculated
!-------------------------------------------------------------------------------
subroutine conv_forw_prop(this, input0, is_train)
    class(ConvLayer)          :: this
    real(kind=8), intent(in)  :: input0(:,:,:,:)
    logical, intent(in)       :: is_train
    real(kind=8), allocatable :: input(:,:,:,:), z_slice(:,:,:)
    integer                   :: i

    if (this%drop_rate > 0 .and. is_train) then
        call this%conv_dropout_rand() ! randomize dropout
        input = input0 * this%drop
    else
        input = input0
    end if

    do i = 1, this%batch_size
        if (this%pad == 'full') then
            ! z(l) = transpose_convolution(a(l-1), k(l)) + b(l)
            call transpose_convolve_3D_kernels(input(:,:,:,i), this%pad, &
                                               this%k, this%stride, z_slice)
        else
            ! z(l) = cross_correlation(a(l-1), k(l)) + b(l)
            call cross_correlate_3D_kernels(input(:,:,:,i), this%pad, this%k, &
                                            this%stride, z_slice)
        end if

        this%z(:,:,:,i) = z_slice + this%b
    end do

    ! output activation cannot be softmax (or any "special case" function),
    ! so we can call all activations in the same way
    ! a(l) = activ(z(l))
    this%a = activfunc(this%z, this%activ)

    ! forward prop activations through pool
    if (associated(this%next_pool)) then
        call this%next_pool%pool_forw_prop(this%a)
    end if

    ! traverse next layers
    if (associated(this%next_layer)) then
        if (associated(this%next_pool)) then
            ! has next_pool; forward prop pooled output
            call this%next_layer%conv_forw_prop(this%next_pool%a, is_train)
        else
            ! forward prop regular output
            call this%next_layer%conv_forw_prop(this%a, is_train)
        end if
    end if
end subroutine

!-------------------------------------------------------------------------------
! helper subroutine to back propagate through ConvLayers
!-------------------------------------------------------------------------------
! this:     (ConvLayer - implicitly passed)
!-------------------------------------------------------------------------------
! alters :: this ConvLayer's d is calculated
!-------------------------------------------------------------------------------
subroutine conv_back_prop(this)
    class(ConvLayer)          :: this
    real(kind=8), allocatable :: d_slice(:,:,:)
    integer                   :: i

    if (.not. associated(this%next_layer)) then
        print *, '-------------------------------------------'
        print *, '(conv_layer_definitions :: conv_back_prop)'
        print *, 'cannot call conv_back_prop on output layer.'
        print *, '-------------------------------------------'
        stop -1
    end if

    do i = 1, this%batch_size
        if (this%next_layer%pad == 'full') then
            call cross_correlate_3D_perms_sum_kernel( &
                                                this%next_layer%d(:,:,:,i), &
                                                this%next_layer%pad, &
                                                this%next_layer%k, &
                                                this%next_layer%stride, d_slice)
        else
            ! derivative d(l+1) wrt a(l) OR pool(l) (if present)
            call transpose_convolve_3D_perms_sum_kernel( &
                                                this%next_layer%d(:,:,:,i), &
                                                this%next_layer%pad, &
                                                this%next_layer%k, &
                                                this%next_layer%stride, d_slice)
        end if

        if (associated(this%next_pool)) then
            ! must undo the pooling to find derivative wrt a's kept by pool
            call this%next_pool%pool_back_prop(d_slice, i, this%d)
        else
            ! no pooling to process
            this%d(:,:,:,i) = d_slice
        end if
    end do

    ! derivative wrt z
    this%d = this%d * activfunc_deriv(this%z, this%activ)

    ! traverse previous layers
    if (associated(this%prev_layer)) then
        call this%prev_layer%conv_back_prop()
    end if
end subroutine

!-------------------------------------------------------------------------------
! helper subroutine to adjust kernels and biases in ConvLayers
!-------------------------------------------------------------------------------
! this:       (ConvLayer - implicitly passed)
! input0:     (real(:,:,:,:)) previous layer activations
! learn_rate: (real) scale factor for change in kernels and biases
! is_train:   (logical) in training iteration
! avg_deltas: (logical) average deltas across batch for update
!-------------------------------------------------------------------------------
! alters ::   this ConvLayer's kernels and biases adjusted to minimize loss
!-------------------------------------------------------------------------------
subroutine conv_update(this, input0, learn_rate, is_train, avg_deltas)
    class(ConvLayer)              :: this
    real(kind=8), intent(in)      :: input0(:,:,:,:)
    real, intent(in)              :: learn_rate
    logical, intent(in)           :: is_train, avg_deltas
    real(kind=8), allocatable     :: input(:,:,:,:), total_k_change(:,:,:,:), &
                                     k_change(:,:,:,:), chan_avgs(:,:,:), &
                                     chan_biases(:)
    real                          :: scale
    integer                       :: i

    if (this%drop_rate > 0 .and. is_train) then
        call this%conv_dropout_rand() ! randomize dropout
        input = input0 * this%drop
    else
        input = input0
    end if

    ! initial call to start counter array (total_k_change)
    if (this%pad == 'full') then
        call cross_correlate_3D_perms_group_base( &
                                            this%d(:,:,:,1), this%pad, &
                                            input(:,:,:,1), this%stride, &
                                            'right', total_k_change)
    else
        call cross_correlate_3D_perms_group_kernel( &
                                            input(:,:,:,1), this%pad, &
                                            this%d(:,:,:,1), this%stride, &
                                            'right', total_k_change)
    end if

    ! sum kernel changes across batch
    do i = 2, this%batch_size
        if (this%pad == 'full') then
            call cross_correlate_3D_perms_group_base( &
                                                this%d(:,:,:,i), this%pad, &
                                                input(:,:,:,i), this%stride, &
                                                'right', k_change)
        else
            call cross_correlate_3D_perms_group_kernel( &
                                                input(:,:,:,i), this%pad, &
                                                this%d(:,:,:,i), this%stride, &
                                                'right', k_change)
        end if

        total_k_change = total_k_change + k_change
    end do

    ! multiply deltas by learn_rate
    scale = learn_rate

    if (avg_deltas) then
        scale = scale / this%batch_size
    end if

    ! weights updated on average change
    this%k = this%k - total_k_change * scale

    ! calculate 1 bias per channel per batch;
    ! start with 4D: (row,col,channel,batch):

    ! sum then average by batch, and scale, 3D: (row,col,channel)
    chan_avgs = sum(this%d, dim=4) * scale

    ! sum then average channels, 1D: (channel), where each val is a bias
    chan_biases = sum(sum(chan_avgs, dim=1), dim=1) / &       ! sum
                  (size(this%d, dim=1) * size(this%d, dim=2)) ! count

    ! one bias per channel -> 2D array; can add bias array instead of loops
    do i = 1, size(this%b, dim=3)
        this%b(:,:,i) = this%b(:,:,i) - chan_biases(i)
    end do

    ! traverse next layers
    if (associated(this%next_layer)) then
        if (associated(this%next_pool)) then
            ! prop pooled output
            call this%next_layer%conv_update(this%next_pool%a, learn_rate, &
                                             is_train, avg_deltas)
        else
            call this%next_layer%conv_update(this%a, learn_rate, &
                                             is_train, avg_deltas)
        end if
    end if
end subroutine
end module

max_stars_repo_name: welchmatt/Neural-Network-Framework

max_issues_repo_name: welchmatt/Neural-Network-Framework

max_forks_repo_name: welchmatt/Neural-Network-Framework