diff --git a/Registry/Registry.EM b/Registry/Registry.EM index abfc5030bd..1ed503b6ee 100644 --- a/Registry/Registry.EM +++ b/Registry/Registry.EM @@ -20,6 +20,7 @@ state real SST_INPUT ij misc 1 - rh "SST_IN # Chem Scalars rconfig integer chem_opt namelist,physics max_domains 0 rh "chem_opt" "" "" +rconfig integer use_ml_activate namelist,physics max_domains 1 rh "use_ml_activate" "Use ML activation (1) or traditional activation (0)" "" state real - ikjftb chem 1 - - - # # Tracer Scalars: If you want full resolved and non-resolved dispersion, compile WRF-Chem diff --git a/Registry/registry.chem b/Registry/registry.chem index 6cd996156b..de3dee320d 100644 --- a/Registry/registry.chem +++ b/Registry/registry.chem @@ -24,12 +24,132 @@ state real ph_aer02 ikj misc 1 - - "p state real ph_aer03 ikj misc 1 - - "ph_aer03" "H+ of aerosol bin 3" "mol/kg" state real ph_aer04 ikj misc 1 - - "ph_aer04" "H+ of aerosol bin 4" "mol/kg" -state real ccn1 ikj misc 1 - r "ccn1" "CCN concentration at S=0.02%" "#/cm3" -state real ccn2 ikj misc 1 - r "ccn2" "CCN concentration at S=0.05%" "#/cm3" -state real ccn3 ikj misc 1 - r "ccn3" "CCN concentration at S=0.1%" "#/cm3" -state real ccn4 ikj misc 1 - r "ccn4" "CCN concentration at S=0.2%" "#/cm3" -state real ccn5 ikj misc 1 - r "ccn5" "CCN concentration at S=0.5%" "#/cm3" -state real ccn6 ikj misc 1 - r "ccn6" "CCN concentration at S=1.0%" "#/cm3" +state real ccn1 ikj misc 1 - rh "ccn1" "CCN concentration at S=0.02%" "#/cm3" +state real ccn2 ikj misc 1 - rh "ccn2" "CCN concentration at S=0.05%" "#/cm3" +state real ccn3 ikj misc 1 - rh "ccn3" "CCN concentration at S=0.1%" "#/cm3" +state real ccn4 ikj misc 1 - rh "ccn4" "CCN concentration at S=0.2%" "#/cm3" +state real ccn5 ikj misc 1 - rh "ccn5" "CCN concentration at S=0.5%" "#/cm3" +state real ccn6 ikj misc 1 - rh "ccn6" "CCN concentration at S=1.0%" "#/cm3" + +state real fn11 ijk misc 1 - rh "fn11" "aerosol activation fraction" " " +state real fn12 ijk misc 1 - rh "fn12" "aerosol activation fraction" " " +state real fn13 ijk misc 1 - rh "fn13" "aerosol activation fraction" " " +state real fn14 ijk misc 1 - rh "fn14" "aerosol activation fraction" " " +state real fn21 ijk misc 1 - rh "fn21" "aerosol activation fraction" " " +state real fn22 ijk misc 1 - rh "fn22" "aerosol activation fraction" " " +state real fn23 ijk misc 1 - rh "fn23" "aerosol activation fraction" " " +state real fn24 ijk misc 1 - rh "fn24" "aerosol activation fraction" " " +state real fn31 ijk misc 1 - rh "fn31" "aerosol activation fraction" " " +state real fn32 ijk misc 1 - rh "fn32" "aerosol activation fraction" " " +state real fn33 ijk misc 1 - rh "fn33" "aerosol activation fraction" " " +state real fn34 ijk misc 1 - rh "fn34" "aerosol activation fraction" " " +state real fn41 ijk misc 1 - rh "fn41" "aerosol activation fraction" " " +state real fn42 ijk misc 1 - rh "fn42" "aerosol activation fraction" " " +state real fn43 ijk misc 1 - rh "fn43" "aerosol activation fraction" " " +state real fn44 ijk misc 1 - rh "fn44" "aerosol activation fraction" " " + +state real na11 ijk misc 1 - rh "na11" "aerosol number concentration" " " +state real na12 ijk misc 1 - rh "na12" "aerosol number concentration" " " +state real na13 ijk misc 1 - rh "na13" "aerosol number concentration" " " +state real na14 ijk misc 1 - rh "na14" "aerosol number concentration" " " +state real na21 ijk misc 1 - rh "na21" "aerosol number concentration" " " +state real na22 ijk misc 1 - rh "na22" "aerosol number concentration" " " +state real na23 ijk misc 1 - rh "na23" "aerosol number concentration" " " +state real na24 ijk misc 1 - rh "na24" "aerosol number concentration" " " +state real na31 ijk misc 1 - rh "na31" "aerosol number concentration" " " +state real na32 ijk misc 1 - rh "na32" "aerosol number concentration" " " +state real na33 ijk misc 1 - rh "na33" "aerosol number concentration" " " +state real na34 ijk misc 1 - rh "na34" "aerosol number concentration" " " +state real na41 ijk misc 1 - rh "na41" "aerosol number concentration" " " +state real na42 ijk misc 1 - rh "na42" "aerosol number concentration" " " +state real na43 ijk misc 1 - rh "na43" "aerosol number concentration" " " +state real na44 ijk misc 1 - rh "na44" "aerosol number concentration" " " + +state real sg11 ijk misc 1 - rh "sg11" "geometric standard deviation of aerosol size distribution" " " +state real sg12 ijk misc 1 - rh "sg12" "geometric standard deviation of aerosol size distribution" " " +state real sg13 ijk misc 1 - rh "sg13" "geometric standard deviation of aerosol size distribution" " " +state real sg14 ijk misc 1 - rh "sg14" "geometric standard deviation of aerosol size distribution" " " +state real sg21 ijk misc 1 - rh "sg21" "geometric standard deviation of aerosol size distribution" " " +state real sg22 ijk misc 1 - rh "sg22" "geometric standard deviation of aerosol size distribution" " " +state real sg23 ijk misc 1 - rh "sg23" "geometric standard deviation of aerosol size distribution" " " +state real sg24 ijk misc 1 - rh "sg24" "geometric standard deviation of aerosol size distribution" " " +state real sg31 ijk misc 1 - rh "sg31" "geometric standard deviation of aerosol size distribution" " " +state real sg32 ijk misc 1 - rh "sg32" "geometric standard deviation of aerosol size distribution" " " +state real sg33 ijk misc 1 - rh "sg33" "geometric standard deviation of aerosol size distribution" " " +state real sg34 ijk misc 1 - rh "sg34" "geometric standard deviation of aerosol size distribution" " " +state real sg41 ijk misc 1 - rh "sg41" "geometric standard deviation of aerosol size distribution" " " +state real sg42 ijk misc 1 - rh "sg42" "geometric standard deviation of aerosol size distribution" " " +state real sg43 ijk misc 1 - rh "sg43" "geometric standard deviation of aerosol size distribution" " " +state real sg44 ijk misc 1 - rh "sg44" "geometric standard deviation of aerosol size distribution" " " + +state real vo11 ijk misc 1 - rh "vo11" "total aerosol volume concentration (m3/m3)" " " +state real vo12 ijk misc 1 - rh "vo12" "total aerosol volume concentration (m3/m3)" " " +state real vo13 ijk misc 1 - rh "vo13" "total aerosol volume concentration (m3/m3)" " " +state real vo14 ijk misc 1 - rh "vo14" "total aerosol volume concentration (m3/m3)" " " +state real vo21 ijk misc 1 - rh "vo21" "total aerosol volume concentration (m3/m3)" " " +state real vo22 ijk misc 1 - rh "vo22" "total aerosol volume concentration (m3/m3)" " " +state real vo23 ijk misc 1 - rh "vo23" "total aerosol volume concentration (m3/m3)" " " +state real vo24 ijk misc 1 - rh "vo24" "total aerosol volume concentration (m3/m3)" " " +state real vo31 ijk misc 1 - rh "vo31" "total aerosol volume concentration (m3/m3)" " " +state real vo32 ijk misc 1 - rh "vo32" "total aerosol volume concentration (m3/m3)" " " +state real vo33 ijk misc 1 - rh "vo33" "total aerosol volume concentration (m3/m3)" " " +state real vo34 ijk misc 1 - rh "vo34" "total aerosol volume concentration (m3/m3)" " " +state real vo41 ijk misc 1 - rh "vo41" "total aerosol volume concentration (m3/m3)" " " +state real vo42 ijk misc 1 - rh "vo42" "total aerosol volume concentration (m3/m3)" " " +state real vo43 ijk misc 1 - rh "vo43" "total aerosol volume concentration (m3/m3)" " " +state real vo44 ijk misc 1 - rh "vo44" "total aerosol volume concentration (m3/m3)" " " + +state real dl11 ijk misc 1 - rh "dl11" "minimum size of section (cm)" " " +state real dl12 ijk misc 1 - rh "dl12" "minimum size of section (cm)" " " +state real dl13 ijk misc 1 - rh "dl13" "minimum size of section (cm)" " " +state real dl14 ijk misc 1 - rh "dl14" "minimum size of section (cm)" " " +state real dl21 ijk misc 1 - rh "dl21" "minimum size of section (cm)" " " +state real dl22 ijk misc 1 - rh "dl22" "minimum size of section (cm)" " " +state real dl23 ijk misc 1 - rh "dl23" "minimum size of section (cm)" " " +state real dl24 ijk misc 1 - rh "dl24" "minimum size of section (cm)" " " +state real dl31 ijk misc 1 - rh "dl31" "minimum size of section (cm)" " " +state real dl32 ijk misc 1 - rh "dl32" "minimum size of section (cm)" " " +state real dl33 ijk misc 1 - rh "dl33" "minimum size of section (cm)" " " +state real dl34 ijk misc 1 - rh "dl34" "minimum size of section (cm)" " " +state real dl41 ijk misc 1 - rh "dl41" "minimum size of section (cm)" " " +state real dl42 ijk misc 1 - rh "dl42" "minimum size of section (cm)" " " +state real dl43 ijk misc 1 - rh "dl43" "minimum size of section (cm)" " " +state real dl44 ijk misc 1 - rh "dl44" "minimum size of section (cm)" " " + +state real dh11 ijk misc 1 - rh "dh11" "maximum size of section (cm)" " " +state real dh12 ijk misc 1 - rh "dh12" "maximum size of section (cm)" " " +state real dh13 ijk misc 1 - rh "dh13" "maximum size of section (cm)" " " +state real dh14 ijk misc 1 - rh "dh14" "maximum size of section (cm)" " " +state real dh21 ijk misc 1 - rh "dh21" "maximum size of section (cm)" " " +state real dh22 ijk misc 1 - rh "dh22" "maximum size of section (cm)" " " +state real dh23 ijk misc 1 - rh "dh23" "maximum size of section (cm)" " " +state real dh24 ijk misc 1 - rh "dh24" "maximum size of section (cm)" " " +state real dh31 ijk misc 1 - rh "dh31" "maximum size of section (cm)" " " +state real dh32 ijk misc 1 - rh "dh32" "maximum size of section (cm)" " " +state real dh33 ijk misc 1 - rh "dh33" "maximum size of section (cm)" " " +state real dh34 ijk misc 1 - rh "dh34" "maximum size of section (cm)" " " +state real dh41 ijk misc 1 - rh "dh41" "maximum size of section (cm)" " " +state real dh42 ijk misc 1 - rh "dh42" "maximum size of section (cm)" " " +state real dh43 ijk misc 1 - rh "dh43" "maximum size of section (cm)" " " +state real dh44 ijk misc 1 - rh "dh44" "maximum size of section (cm)" " " + +state real hg11 ijk misc 1 - rh "hg11" "bulk hygroscopicity of aerosol mode" " " +state real hg12 ijk misc 1 - rh "hg12" "bulk hygroscopicity of aerosol mode" " " +state real hg13 ijk misc 1 - rh "hg13" "bulk hygroscopicity of aerosol mode" " " +state real hg14 ijk misc 1 - rh "hg14" "bulk hygroscopicity of aerosol mode" " " +state real hg21 ijk misc 1 - rh "hg21" "bulk hygroscopicity of aerosol mode" " " +state real hg22 ijk misc 1 - rh "hg22" "bulk hygroscopicity of aerosol mode" " " +state real hg23 ijk misc 1 - rh "hg23" "bulk hygroscopicity of aerosol mode" " " +state real hg24 ijk misc 1 - rh "hg24" "bulk hygroscopicity of aerosol mode" " " +state real hg31 ijk misc 1 - rh "hg31" "bulk hygroscopicity of aerosol mode" " " +state real hg32 ijk misc 1 - rh "hg32" "bulk hygroscopicity of aerosol mode" " " +state real hg33 ijk misc 1 - rh "hg33" "bulk hygroscopicity of aerosol mode" " " +state real hg34 ijk misc 1 - rh "hg34" "bulk hygroscopicity of aerosol mode" " " +state real hg41 ijk misc 1 - rh "hg41" "bulk hygroscopicity of aerosol mode" " " +state real hg42 ijk misc 1 - rh "hg42" "bulk hygroscopicity of aerosol mode" " " +state real hg43 ijk misc 1 - rh "hg43" "bulk hygroscopicity of aerosol mode" " " +state real hg44 ijk misc 1 - rh "hg44" "bulk hygroscopicity of aerosol mode" " " + # cloud water fractional removal rate needed for wet scavenging state real qlsink ikj misc 1 - rdu "qlsink" "CLOUD WATER SINK" "/S" diff --git a/chem/chem_driver.F b/chem/chem_driver.F index 8650b9444a..0453f28c1d 100755 --- a/chem/chem_driver.F +++ b/chem/chem_driver.F @@ -1021,7 +1021,7 @@ end SUBROUTINE sum_pm_driver end if #endif endif - + call wrf_debug(15,'Liran start vertical mixing with dry deposition') ! ! do vertical mixing with dry deposition ! @@ -1048,8 +1048,22 @@ end SUBROUTINE sum_pm_driver call dry_dep_driver(grid%id,curr_secs,ktau,grid%dt,config_flags, & grid%gmt,ijulian,t_phy,moist,scalar,p8w,t8w,vvel, & rri,p_phy,chem,tracer,rho,dz8w,rh,grid%exch_h,grid%hfx,grid%dx, & - grid%cldfra, grid%cldfra_old,grid%raincv_b,seasin,dustin, & - grid%ccn1, grid%ccn2, grid%ccn3, grid%ccn4, grid%ccn5, grid%ccn6, & + grid%cldfra,grid%cldfra_old,grid%raincv_b,seasin,dustin, & + grid%fn11,grid%fn12,grid%fn13,grid%fn14,grid%fn21,grid%fn22,grid%fn23,grid%fn24,& + grid%fn31,grid%fn32,grid%fn33,grid%fn34,grid%fn41,grid%fn42,grid%fn43,grid%fn44,& + grid%na11,grid%na12,grid%na13,grid%na14,grid%na21,grid%na22,grid%na23,grid%na24,& + grid%na31,grid%na32,grid%na33,grid%na34,grid%na41,grid%na42,grid%na43,grid%na44,& + grid%vo11,grid%vo12,grid%vo13,grid%vo14,grid%vo21,grid%vo22,grid%vo23,grid%vo24,& + grid%vo31,grid%vo32,grid%vo33,grid%vo34,grid%vo41,grid%vo42,grid%vo43,grid%vo44,& + grid%dl11,grid%dl12,grid%dl13,grid%dl14,grid%dl21,grid%dl22,grid%dl23,grid%dl24,& + grid%dl31,grid%dl32,grid%dl33,grid%dl34,grid%dl41,grid%dl42,grid%dl43,grid%dl44,& + grid%dh11,grid%dh12,grid%dh13,grid%dh14,grid%dh21,grid%dh22,grid%dh23,grid%dh24,& + grid%dh31,grid%dh32,grid%dh33,grid%dh34,grid%dh41,grid%dh42,grid%dh43,grid%dh44,& + grid%sg11,grid%sg12,grid%sg13,grid%sg14,grid%sg21,grid%sg22,grid%sg23,grid%sg24,& + grid%sg31,grid%sg32,grid%sg33,grid%sg34,grid%sg41,grid%sg42,grid%sg43,grid%sg44,& + grid%hg11,grid%hg12,grid%hg13,grid%hg14,grid%hg21,grid%hg22,grid%hg23,grid%hg24,& + grid%hg31,grid%hg32,grid%hg33,grid%hg34,grid%hg41,grid%hg42,grid%hg43,grid%hg44,& + grid%ccn1, grid%ccn2, grid%ccn3, grid%ccn4, grid%ccn5, grid%ccn6, & grid%qndropsource,grid%ivgtyp,grid%tsk,grid%gsw,grid%vegfra,pbl_h, & grid%rmol,grid%ust,grid%znt,grid%xlat,grid%xlong, & zmid,z_at_w,grid%xland,grid%ash_fall, & diff --git a/chem/dry_dep_driver.F b/chem/dry_dep_driver.F index b5c2f66501..47dc333c96 100755 --- a/chem/dry_dep_driver.F +++ b/chem/dry_dep_driver.F @@ -14,7 +14,21 @@ MODULE module_dry_dep_driver subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & gmt,julday,t_phy,moist,scalar,p8w,t8w,w,alt, & p_phy,chem,tracer,rho_phy,dz8w,rh,exch_h,hfx,dx, & - cldfra, cldfra_old,raincv,seasin,dustin, & + cldfra,cldfra_old,raincv,seasin,dustin,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ivgtyp,tsk,gsw,vegfra,pbl,rmol,ust,znt,xlat,xlong,z,z_at_w,& xland,ash_fall,h2oaj,h2oai,nu3,ac3,cor3,asulf,ahno3, & @@ -144,7 +158,20 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & REAL, INTENT(OUT), dimension(ims:ime,kms:kme,jms:jme) :: nsource, & ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat - + REAL, INTENT(OUT), dimension(ims:ime,kms:kme,jms:jme) :: fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 REAL, INTENT(IN ) :: & dtstep,gmt,dx @@ -1068,18 +1095,31 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & end where qsrflx(:,:,:) = 0.0 - + CALL wrf_debug(15,'Liran test aerosol') mixactivate_select: SELECT CASE(config_flags%chem_opt) - CASE (RADM2SORG_AQ, RADM2SORG_AQCHEM, RACMSORG_AQ, RACMSORG_AQCHEM_KPP, RACM_ESRLSORG_AQCHEM_KPP, CBMZSORG_AQ, & CB05_SORG_AQ_KPP) - CALL wrf_debug(15,'call mixactivate for sorgam aerosol') + CALL wrf_debug(15,'call mixactivate for LR sorgam aerosol') call sorgam_mixactivate ( & id, ktau, dtstep, config_flags, idrydep_onoff, & dryrho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p8w, t8w, exch_h, & moist(ims,kms,jms,P_QV), moist(ims,kms,jms,P_QC), moist(ims,kms,jms,P_QI), & - scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem, & + scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -1091,7 +1131,21 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & dryrho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p8w, t8w, exch_h, & moist(ims,kms,jms,P_QV), moist(ims,kms,jms,P_QC),moist(ims,kms,jms,P_QI), & - scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem, & + scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -1104,6 +1158,20 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & ddvel, z, dz8w, p8w, t8w, exch_h, & moist(ims,kms,jms,P_QV), moist(ims,kms,jms,P_QC), moist(ims,kms,jms,P_QI), & scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem, & + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -1118,7 +1186,21 @@ subroutine dry_dep_driver(id,curr_secs,ktau,dtstep,config_flags, & dryrho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p8w, t8w, exch_h, & moist(ims,kms,jms,P_QV), moist(ims,kms,jms,P_QC), moist(ims,kms,jms,P_QI), & - scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem, & + scalar(ims,kms,jms,P_QNDROP), f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & qsrflx, & ids,ide, jds,jde, kds,kde, & diff --git a/chem/module_mixactivate_wrappers.F b/chem/module_mixactivate_wrappers.F index 71d5292b6e..6375ea6809 100644 --- a/chem/module_mixactivate_wrappers.F +++ b/chem/module_mixactivate_wrappers.F @@ -23,7 +23,21 @@ subroutine mosaic_mixactivate ( & id, ktau, dtstep, config_flags, idrydep_onoff, & rho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p_at_w, t_at_w, exch_h, & - qv, qc, qi, qndrop3d, f_qc, f_qi, chem, & + qv, qc, qi, qndrop3d, f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & qsrflx, & ids,ide, jds,jde, kds,kde, & @@ -75,7 +89,20 @@ subroutine mosaic_mixactivate ( & chem real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource,& ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat - + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 type(grid_config_rec_type), intent(in) :: config_flags real, intent(out) :: qsrflx(ims:ime, jms:jme, num_chem) ! wet deposition flux of aerosol ! local vars @@ -119,9 +146,23 @@ subroutine mosaic_mixactivate ( & its,ite, jts,jte, kts,kte, & rho_phy, z, dz8w, p_at_w, t_at_w, exch_h, & cldfra, cldfra_old, qsrflx, & + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & id, ktau, dtstep, & - f_qc, f_qi ) + f_qc, f_qi, config_flags ) end subroutine mosaic_mixactivate @@ -207,7 +248,21 @@ subroutine sorgam_mixactivate ( & id, ktau, dtstep, config_flags, idrydep_onoff, & rho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p_at_w, t_at_w, exch_h, & - qv, qc, qi, qndrop3d, f_qc, f_qi, chem, & + qv, qc, qi, qndrop3d, f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -258,7 +313,20 @@ subroutine sorgam_mixactivate ( & chem real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, & ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat - + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 type(grid_config_rec_type), intent(in) :: config_flags ! local vars @@ -301,10 +369,24 @@ subroutine sorgam_mixactivate ( & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & rho_phy, z, dz8w, p_at_w, t_at_w, exch_h, & - cldfra, cldfra_old, qsrflx, & + cldfra, cldfra_old, qsrflx,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & id, ktau, dtstep, & - f_qc, f_qi ) + f_qc, f_qi, config_flags ) end subroutine sorgam_mixactivate @@ -313,7 +395,21 @@ subroutine soa_vbs_mixactivate ( & id, ktau, dtstep, config_flags, idrydep_onoff, & rho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p_at_w, t_at_w, exch_h, & - qv, qc, qi, qndrop3d, f_qc, f_qi, chem, & + qv, qc, qi, qndrop3d, f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -364,7 +460,20 @@ subroutine soa_vbs_mixactivate ( & chem real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, & ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat - + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 type(grid_config_rec_type), intent(in) :: config_flags ! local vars real qsrflx(ims:ime, jms:jme, num_chem) ! dry deposition flux of aerosol real sumhygro,sumvol @@ -405,10 +514,24 @@ subroutine soa_vbs_mixactivate ( & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & rho_phy, z, dz8w, p_at_w, t_at_w, exch_h, & - cldfra, cldfra_old, qsrflx, & + cldfra, cldfra_old, qsrflx,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & id, ktau, dtstep, & - f_qc, f_qi ) + f_qc, f_qi, config_flags ) end subroutine soa_vbs_mixactivate @@ -416,7 +539,21 @@ subroutine sorgam_vbs_mixactivate ( & id, ktau, dtstep, config_flags, idrydep_onoff, & rho_phy, t_phy, w, cldfra, cldfra_old, & ddvel, z, dz8w, p_at_w, t_at_w, exch_h, & - qv, qc, qi, qndrop3d, f_qc, f_qi, chem, & + qv, qc, qi, qndrop3d, f_qc, f_qi, chem,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -467,7 +604,20 @@ subroutine sorgam_vbs_mixactivate ( & chem real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, & ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat - + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 type(grid_config_rec_type), intent(in) :: config_flags ! local vars @@ -510,10 +660,24 @@ subroutine sorgam_vbs_mixactivate ( & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & rho_phy, z, dz8w, p_at_w, t_at_w, exch_h, & - cldfra, cldfra_old, qsrflx, & + cldfra, cldfra_old, qsrflx,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & id, ktau, dtstep, & - f_qc, f_qi ) + f_qc, f_qi, config_flags ) end subroutine sorgam_vbs_mixactivate diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index 53514d346e..9fc09a32b4 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -928,7 +928,7 @@ SUBROUTINE microphysics_driver( & 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 ) + F_QC=f_qc, F_QI=f_qi, CONFIG_FLAGS=config_flags ) END IF ELSEIF ( (chem_opt==0 .OR. chem_opt==401) .AND. progn==1 .AND. & (mp_physics==NSSL_2MOM .and. config_flags%nssl_2moment_on==1)) THEN diff --git a/phys/module_mixactivate.F b/phys/module_mixactivate.F index a44c0cc279..ac1c311f99 100644 --- a/phys/module_mixactivate.F +++ b/phys/module_mixactivate.F @@ -11,7 +11,7 @@ MODULE module_mixactivate PRIVATE -PUBLIC prescribe_aerosol_mixactivate, mixactivate, activate !BSINGH - added 'activate' for WRFCuP scheme +PUBLIC prescribe_aerosol_mixactivate, mixactivate, activate, activate_ml CONTAINS @@ -20,109 +20,125 @@ MODULE module_mixactivate ! 06-nov-2005 rce - grid_id & ktau added to arg list ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3) subroutine prescribe_aerosol_mixactivate ( & - grid_id, ktau, dtstep, naer, & - ccn_conc, chem_opt, & ! RAS - rho_phy, th_phy, pi_phy, w, cldfra, cldfra_old, & - z, dz8w, p_at_w, t_at_w, exch_h, & + grid_id, ktau, dtstep, naer, & + ccn_conc, chem_opt, & ! RAS + rho_phy, th_phy, pi_phy, w, cldfra, cldfra_old, & + z, dz8w, p_at_w, t_at_w, exch_h, & qv, qc, qi, qndrop3d, & nsource, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, & - f_qc, f_qi, cn ) + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte, & + f_qc, f_qi, cn, config_flags ) -! USE module_configure + USE module_configure, only: grid_config_rec_type ! wrapper to call mixactivate for mosaic description of aerosol - implicit none + implicit none ! subr arguments - integer, intent(in) :: & - grid_id, ktau, & - ids, ide, jds, jde, kds, kde, & - ims, ime, jms, jme, kms, kme, & - its, ite, jts, jte, kts, kte - - real, intent(in) :: dtstep - real, intent(inout) :: naer ! aerosol number (/kg) + integer, intent(in) :: & + grid_id, ktau, & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + its, ite, jts, jte, kts, kte + + real, intent(in) :: dtstep + real, intent(inout) :: naer ! aerosol number (/kg) real, intent(in) :: ccn_conc ! CCN conc set within namelist integer, optional, intent(in) :: chem_opt - real, intent(in), & - dimension( ims:ime, kms:kme, jms:jme ) :: & - rho_phy, th_phy, pi_phy, w, & - z, dz8w, p_at_w, t_at_w, exch_h + real, intent(in), & + dimension( ims:ime, kms:kme, jms:jme ) :: & + rho_phy, th_phy, pi_phy, w, & + z, dz8w, p_at_w, t_at_w, exch_h - real, intent(inout), & - dimension( ims:ime, kms:kme, jms:jme ) :: cldfra, cldfra_old + real, intent(inout), & + dimension( ims:ime, kms:kme, jms:jme ) :: cldfra, cldfra_old - real, intent(in), & - dimension( ims:ime, kms:kme, jms:jme ) :: & - qv, qc, qi + real, intent(in), & + dimension( ims:ime, kms:kme, jms:jme ) :: & + qv, qc, qi - real, intent(inout), optional, & - dimension( ims:ime, kms:kme, jms:jme ) :: & - cn ! single-size CCN concentration + real, intent(inout), optional, & + dimension( ims:ime, kms:kme, jms:jme ) :: & + cn ! single-size CCN concentration - real, intent(inout), & - dimension( ims:ime, kms:kme, jms:jme ) :: & - qndrop3d + real, intent(inout), & + dimension( ims:ime, kms:kme, jms:jme ) :: & + qndrop3d - real, intent(out), & - dimension( ims:ime, kms:kme, jms:jme) :: nsource + real, intent(out), & + dimension( ims:ime, kms:kme, jms:jme) :: nsource LOGICAL, OPTIONAL :: f_qc, f_qi + type(grid_config_rec_type), intent(in) :: config_flags + ! local vars - integer maxd_aphase, maxd_atype, maxd_asize, maxd_acomp, max_chem - parameter (maxd_aphase=2,maxd_atype=1,maxd_asize=1,maxd_acomp=1, max_chem=10) - real ddvel(its:ite, jts:jte, max_chem) ! dry deposition velosity - real qsrflx(ims:ime, jms:jme, max_chem) ! dry deposition flux of aerosol - real chem(ims:ime, kms:kme, jms:jme, max_chem) ! chem array - integer i,j,k,l,m,n,p - real hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk - integer ntype_aer, nsize_aer(maxd_atype),ncomp_aer(maxd_atype), nphase_aer - integer massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), & - waterptr_aer( maxd_asize, maxd_atype ), & - numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), & - ai_phase, cw_phase + integer maxd_aphase, maxd_atype, maxd_asize, maxd_acomp, max_chem + parameter (maxd_aphase=2,maxd_atype=1,maxd_asize=1,maxd_acomp=1, max_chem=10) + real ddvel(its:ite, jts:jte, max_chem) ! dry deposition velosity + real qsrflx(ims:ime, jms:jme, max_chem) ! dry deposition flux of aerosol + real chem(ims:ime, kms:kme, jms:jme, max_chem) ! chem array + integer i,j,k,l,m,n,p + real hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk + integer ntype_aer, nsize_aer(maxd_atype),ncomp_aer(maxd_atype), nphase_aer + integer massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), & + waterptr_aer( maxd_asize, maxd_atype ), & + numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), & + ai_phase, cw_phase real dlo_sect( maxd_asize, maxd_atype ), & ! minimum size of section (cm) dhi_sect( maxd_asize, maxd_atype ), & ! maximum size of section (cm) - sigmag_aer(maxd_asize, maxd_atype), & ! geometric standard deviation of aerosol size dist - dgnum_aer(maxd_asize, maxd_atype), & ! median diameter (cm) of number distrib of mode - dens_aer( maxd_acomp, maxd_atype), & ! density (g/cm3) of material - mw_aer( maxd_acomp, maxd_atype), & ! molecular weight (g/mole) - dpvolmean_aer(maxd_asize, maxd_atype) ! mean-volume diameter (cm) of mode + sigmag_aer(maxd_asize, maxd_atype), & ! geometric standard deviation of aerosol size dist + dgnum_aer(maxd_asize, maxd_atype), & ! median diameter (cm) of number distrib of mode + dens_aer( maxd_acomp, maxd_atype), & ! density (g/cm3) of material + mw_aer( maxd_acomp, maxd_atype), & ! molecular weight (g/mole) + dpvolmean_aer(maxd_asize, maxd_atype) ! mean-volume diameter (cm) of mode ! terminology: (pi/6) * (mean-volume diameter)**3 == ! (volume mixing ratio of section/mode)/(number mixing ratio) - real, dimension(ims:ime,kms:kme,jms:jme) :: & - ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat - integer idrydep_onoff - real, dimension(ims:ime,kms:kme,jms:jme) :: t_phy - integer msectional - - - integer ptr - real maer + real, dimension(ims:ime,kms:kme,jms:jme) :: & + ccn1,ccn2,ccn3,ccn4,ccn5,ccn6,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 + integer idrydep_onoff + real, dimension(ims:ime,kms:kme,jms:jme) :: t_phy + integer msectional + + + integer ptr + real maer ! if(naer.lt.1.)then -! naer=1000.e6 ! #/kg default value +! naer=1000.e6 ! #/kg default value ! endif IF ( (naer.lt.1.) .OR. ( PRESENT(chem_opt) .AND. (chem_opt.eq.401))) THEN naer = ccn_conc !CCN value set in namelist ENDIF - ai_phase=1 - cw_phase=2 - idrydep_onoff = 0 - msectional = 0 + ai_phase=1 + cw_phase=2 + idrydep_onoff = 0 + msectional = 0 - t_phy(its:ite,kts:kte,jts:jte)=th_phy(its:ite,kts:kte,jts:jte)*pi_phy(its:ite,kts:kte,jts:jte) + t_phy(its:ite,kts:kte,jts:jte)=th_phy(its:ite,kts:kte,jts:jte)*pi_phy(its:ite,kts:kte,jts:jte) ntype_aer=maxd_atype do n=1,ntype_aer nsize_aer(n)=maxd_asize - ncomp_aer(n)=maxd_acomp + ncomp_aer(n)=maxd_acomp end do nphase_aer=maxd_aphase @@ -131,14 +147,14 @@ subroutine prescribe_aerosol_mixactivate ( & do m=1,nsize_aer(n) dlo_sect( m,n )=0.01e-4 ! minimum size of section (cm) dhi_sect( m,n )=0.5e-4 ! maximum size of section (cm) - sigmag_aer(m,n)=2. ! geometric standard deviation of aerosol size dist - dgnum_aer(m,n)=0.1e-4 ! median diameter (cm) of number distrib of mode - dpvolmean_aer(m,n) = dgnum_aer(m,n) * exp( 1.5 * (log(sigmag_aer(m,n)))**2 ) - end do - do l=1,ncomp_aer(n) - dens_aer( l, n)=1.0 ! density (g/cm3) of material - mw_aer( l, n)=132. ! molecular weight (g/mole) - end do + sigmag_aer(m,n)=2. ! geometric standard deviation of aerosol size dist + dgnum_aer(m,n)=0.1e-4 ! median diameter (cm) of number distrib of mode + dpvolmean_aer(m,n) = dgnum_aer(m,n) * exp( 1.5 * (log(sigmag_aer(m,n)))**2 ) + end do + do l=1,ncomp_aer(n) + dens_aer( l, n)=1.0 ! density (g/cm3) of material + mw_aer( l, n)=132. ! molecular weight (g/mole) + end do end do ptr=0 do p=1,nphase_aer @@ -162,16 +178,16 @@ subroutine prescribe_aerosol_mixactivate ( & do p=1,maxd_aphase do n=1,ntype_aer do m=1,nsize_aer(n) - do l=1,ncomp_aer(n) + do l=1,ncomp_aer(n) ptr=ptr+1 - if(ptr.gt.max_chem)then - write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate' - call wrf_error_fatal("1") - endif - massptr_aer(l, m, n, p)=ptr + if(ptr.gt.max_chem)then + write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate' + call wrf_error_fatal("1") + endif + massptr_aer(l, m, n, p)=ptr ! maer is ug/kg-air; naer is #/kg-air; dgnum is cm; dens_aer is g/cm3 ! 1.e6 factor converts g to ug - maer= 1.0e6 * naer * dens_aer(l,n) * ( (3.1416/6.) * & + maer= 1.0e6 * naer * dens_aer(l,n) * ( (3.1416/6.) * & (dgnum_aer(m,n)**3) * exp( 4.5*((log(sigmag_aer(m,n)))**2) ) ) if(p.eq.ai_phase)then IF ( present( cn ) ) THEN @@ -185,21 +201,21 @@ subroutine prescribe_aerosol_mixactivate ( & chem(its:ite,kts:kte,jts:jte,ptr)=0. endif end do - end do ! size - end do ! type - end do ! phase + end do ! size + end do ! type + end do ! phase do n=1,ntype_aer do m=1,nsize_aer(n) ptr=ptr+1 - if(ptr.gt.max_chem)then - write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate' - call wrf_error_fatal("1") - endif -!wig waterptr_aer(m, n)=ptr - waterptr_aer(m, n)=-1 - end do ! size - end do ! type - ddvel(its:ite,jts:jte,:)=0. + if(ptr.gt.max_chem)then + write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate' + call wrf_error_fatal("1") + endif +!wig waterptr_aer(m, n)=ptr + waterptr_aer(m, n)=-1 + end do ! size + end do ! type + ddvel(its:ite,jts:jte,:)=0. hygro(its:ite,kts:kte,jts:jte,:,:) = 0.5 ! 06-nov-2005 rce - grid_id & ktau added to arg list @@ -215,10 +231,24 @@ subroutine prescribe_aerosol_mixactivate ( & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & rho_phy, z, dz8w, p_at_w, t_at_w, exch_h, & - cldfra, cldfra_old, qsrflx, & + cldfra, cldfra_old, qsrflx,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & grid_id, ktau, dtstep, & - F_QC=f_qc, F_QI=f_qi ) + F_QC=f_qc, F_QI=f_qi, CONFIG_FLAGS=config_flags ) ! ERM : If CCN field was passed in, then copy back the new field, which is the first @@ -246,10 +276,24 @@ subroutine mixactivate( msectional, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & rho, zm, dz8w, p_at_w, t_at_w, kvh, & - cldfra, cldfra_old, qsrflx, & + cldfra, cldfra_old, qsrflx,& + fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24,& + fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44,& + na11,na12,na13,na14,na21,na22,na23,na24,& + na31,na32,na33,na34,na41,na42,na43,na44,& + vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24,& + vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44,& + dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24,& + dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44,& + dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24,& + dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44,& + sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24,& + sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44,& + hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24,& + hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44,& ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, & grid_id, ktau, dtstep, & - f_qc, f_qi ) + f_qc, f_qi, config_flags ) ! vertical diffusion and nucleation of cloud droplets @@ -258,6 +302,7 @@ subroutine mixactivate( msectional, & USE module_model_constants, only: g, rhowater, xlv, cp, rvovrd, r_d, r_v, mwdry, ep_2 USE module_radiation_driver, only: cal_cldfra2 + USE module_configure, only: grid_config_rec_type implicit none @@ -316,9 +361,24 @@ subroutine mixactivate( msectional, & real, intent(in) :: hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk hygroscopicity & REAL, intent(out), DIMENSION( ims:ime, jms:jme, num_chem ) :: qsrflx ! dry deposition rate for aerosol + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: fn11,fn12,fn13,fn14,fn21,fn22,fn23,fn24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: fn31,fn32,fn33,fn34,fn41,fn42,fn43,fn44 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: na11,na12,na13,na14,na21,na22,na23,na24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: na31,na32,na33,na34,na41,na42,na43,na44 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: vo11,vo12,vo13,vo14,vo21,vo22,vo23,vo24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: vo31,vo32,vo33,vo34,vo41,vo42,vo43,vo44 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: dl11,dl12,dl13,dl14,dl21,dl22,dl23,dl24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: dl31,dl32,dl33,dl34,dl41,dl42,dl43,dl44 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: dh11,dh12,dh13,dh14,dh21,dh22,dh23,dh24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: dh31,dh32,dh33,dh34,dh41,dh42,dh43,dh44 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: sg11,sg12,sg13,sg14,sg21,sg22,sg23,sg24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: sg31,sg32,sg33,sg34,sg41,sg42,sg43,sg44 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: hg11,hg12,hg13,hg14,hg21,hg22,hg23,hg24 + real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: hg31,hg32,hg33,hg34,hg41,hg42,hg43,hg44 real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, & ! droplet number source (#/kg/s) ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat + type(grid_config_rec_type), intent(in) :: config_flags !--------------------Local storage------------------------------------- ! @@ -357,6 +417,7 @@ subroutine mixactivate( msectional, & real fac_srflx real depvel_drop, depvel_tmp real, parameter :: depvel_uplimit = 1.0 ! upper limit for dep vels (m/s) + real, parameter :: eps = 0.1 real :: surfrate(num_chem) ! surface exchange rate (/s) real surfratemax ! max surfrate for all species treated here real surfrate_drop ! surfade exchange rate for droplelts @@ -510,15 +571,15 @@ subroutine mixactivate( msectional, & ! initialization for current i ......................................... do k=kts+1,kte - zs(k)=1./(zm(i,k,j)-zm(i,k-1,j)) - enddo - zs(kts)=zs(kts+1) + zs(k)=1./(zm(i,k,j)-zm(i,k-1,j)) + enddo + zs(kts)=zs(kts+1) zs(kte+1)=0. do k=kts,kte -!!$ if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then -!!$! call wrf_error_fatal("1") -!!$ endif +!!$ if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then +!!$! call wrf_error_fatal("1") +!!$ endif if(f_qi)then qcld=qc(i,k,j)+qi(i,k,j) else @@ -536,7 +597,7 @@ subroutine mixactivate( msectional, & lcldfra_old(k)=0. endif qndrop(k)=qndrop3d(i,k,j) -! qndrop(k)=1.e5 +! qndrop(k)=1.e5 cs(k)=rho(i,k,j) ! air density (kg/m3) dz(k)=dz8w(i,k,j) do n=1,ntype_aer @@ -710,11 +771,169 @@ subroutine mixactivate( msectional, & enddo ! 06-nov-2005 rce - grid_id & ktau added to arg list - call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), & - msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & - naerosol, vaerosol, & - dlo_sect,dhi_sect,sigmag_aer,hygro_aer, & - fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k ) +! Check namelist flag to decide between ML activation or traditional activation + if (config_flags%use_ml_activate == 1) then + call activate_ml(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), qv(i,k,j), & + msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & + naerosol, vaerosol, & + dlo_sect,dhi_sect,sigmag_aer,hygro_aer, & + fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k ) + else + call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), & + msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & + naerosol, vaerosol, & + dlo_sect,dhi_sect,sigmag_aer,hygro_aer, & + fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k ) + endif + do n = 1,ntype_aer + do m = 1,nsize_aer(n) + if ((abs(m - 1.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn11(i,k,j)=fn(m,n) + na11(i,k,j)=naerosol(m,n) + vo11(i,k,j)=vaerosol(m,n) + dl11(i,k,j)=dlo_sect(m,n) + dh11(i,k,j)=dhi_sect(m,n) + sg11(i,k,j)=sigmag_aer(m,n) + hg11(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 1.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn12(i,k,j)=fn(m,n) + na12(i,k,j)=naerosol(m,n) + vo12(i,k,j)=vaerosol(m,n) + dl12(i,k,j)=dlo_sect(m,n) + dh12(i,k,j)=dhi_sect(m,n) + sg12(i,k,j)=sigmag_aer(m,n) + hg12(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 1.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn13(i,k,j)=fn(m,n) + na13(i,k,j)=naerosol(m,n) + vo13(i,k,j)=vaerosol(m,n) + dl13(i,k,j)=dlo_sect(m,n) + dh13(i,k,j)=dhi_sect(m,n) + sg13(i,k,j)=sigmag_aer(m,n) + hg13(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 1.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn14(i,k,j)=fn(m,n) + na14(i,k,j)=naerosol(m,n) + vo14(i,k,j)=vaerosol(m,n) + dl14(i,k,j)=dlo_sect(m,n) + dh14(i,k,j)=dhi_sect(m,n) + sg14(i,k,j)=sigmag_aer(m,n) + hg14(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn21(i,k,j)=fn(m,n) + na21(i,k,j)=naerosol(m,n) + vo21(i,k,j)=vaerosol(m,n) + dl21(i,k,j)=dlo_sect(m,n) + dh21(i,k,j)=dhi_sect(m,n) + sg21(i,k,j)=sigmag_aer(m,n) + hg21(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn22(i,k,j)=fn(m,n) + na22(i,k,j)=naerosol(m,n) + vo22(i,k,j)=vaerosol(m,n) + dl22(i,k,j)=dlo_sect(m,n) + dh22(i,k,j)=dhi_sect(m,n) + sg22(i,k,j)=sigmag_aer(m,n) + hg22(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn23(i,k,j)=fn(m,n) + na23(i,k,j)=naerosol(m,n) + vo23(i,k,j)=vaerosol(m,n) + dl23(i,k,j)=dlo_sect(m,n) + dh23(i,k,j)=dhi_sect(m,n) + sg23(i,k,j)=sigmag_aer(m,n) + hg23(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn24(i,k,j)=fn(m,n) + na24(i,k,j)=naerosol(m,n) + vo24(i,k,j)=vaerosol(m,n) + dl24(i,k,j)=dlo_sect(m,n) + dh24(i,k,j)=dhi_sect(m,n) + sg24(i,k,j)=sigmag_aer(m,n) + hg24(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn31(i,k,j)=fn(m,n) + na31(i,k,j)=naerosol(m,n) + vo31(i,k,j)=vaerosol(m,n) + dl31(i,k,j)=dlo_sect(m,n) + dh31(i,k,j)=dhi_sect(m,n) + sg31(i,k,j)=sigmag_aer(m,n) + hg31(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn32(i,k,j)=fn(m,n) + na32(i,k,j)=naerosol(m,n) + vo32(i,k,j)=vaerosol(m,n) + dl32(i,k,j)=dlo_sect(m,n) + dh32(i,k,j)=dhi_sect(m,n) + sg32(i,k,j)=sigmag_aer(m,n) + hg32(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn33(i,k,j)=fn(m,n) + na33(i,k,j)=naerosol(m,n) + vo33(i,k,j)=vaerosol(m,n) + dl33(i,k,j)=dlo_sect(m,n) + dh33(i,k,j)=dhi_sect(m,n) + sg33(i,k,j)=sigmag_aer(m,n) + hg33(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn34(i,k,j)=fn(m,n) + na34(i,k,j)=naerosol(m,n) + vo34(i,k,j)=vaerosol(m,n) + dl34(i,k,j)=dlo_sect(m,n) + dh34(i,k,j)=dhi_sect(m,n) + sg34(i,k,j)=sigmag_aer(m,n) + hg34(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn41(i,k,j)=fn(m,n) + na41(i,k,j)=naerosol(m,n) + vo41(i,k,j)=vaerosol(m,n) + dl41(i,k,j)=dlo_sect(m,n) + dh41(i,k,j)=dhi_sect(m,n) + sg41(i,k,j)=sigmag_aer(m,n) + hg41(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn42(i,k,j)=fn(m,n) + na42(i,k,j)=naerosol(m,n) + vo42(i,k,j)=vaerosol(m,n) + dl42(i,k,j)=dlo_sect(m,n) + dh42(i,k,j)=dhi_sect(m,n) + sg42(i,k,j)=sigmag_aer(m,n) + hg42(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn43(i,k,j)=fn(m,n) + na43(i,k,j)=naerosol(m,n) + vo43(i,k,j)=vaerosol(m,n) + dl43(i,k,j)=dlo_sect(m,n) + dh43(i,k,j)=dhi_sect(m,n) + sg43(i,k,j)=sigmag_aer(m,n) + hg43(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn44(i,k,j)=fn(m,n) + na44(i,k,j)=naerosol(m,n) + vo44(i,k,j)=vaerosol(m,n) + dl44(i,k,j)=dlo_sect(m,n) + dh44(i,k,j)=dhi_sect(m,n) + sg44(i,k,j)=sigmag_aer(m,n) + hg44(i,k,j)=hygro_aer(m,n) + endif + enddo + enddo + do n = 1,ntype_aer do m = 1,nsize_aer(n) @@ -842,11 +1061,169 @@ subroutine mixactivate( msectional, & enddo ! print *,'old cloud wbar,wmix=',wbar,wmix - call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), & - msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & - naerosol, vaerosol, & - dlo_sect,dhi_sect, sigmag_aer,hygro_aer, & - fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k ) +! Check namelist flag to decide between ML activation or traditional activation + if (config_flags%use_ml_activate == 1) then + call activate_ml(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), qv(i,k,j),& + msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & + naerosol, vaerosol, & + dlo_sect,dhi_sect, sigmag_aer,hygro_aer, & + fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k ) + else + call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), & + msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & + naerosol, vaerosol, & + dlo_sect,dhi_sect, sigmag_aer,hygro_aer, & + fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k ) + endif + do n=1,ntype_aer + do m=1,nsize_aer(n) + if ((abs(m - 1.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn11(i,k,j)=fn(m,n) + na11(i,k,j)=naerosol(m,n) + vo11(i,k,j)=vaerosol(m,n) + dl11(i,k,j)=dlo_sect(m,n) + dh11(i,k,j)=dhi_sect(m,n) + sg11(i,k,j)=sigmag_aer(m,n) + hg11(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 1.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn12(i,k,j)=fn(m,n) + na12(i,k,j)=naerosol(m,n) + vo12(i,k,j)=vaerosol(m,n) + dl12(i,k,j)=dlo_sect(m,n) + dh12(i,k,j)=dhi_sect(m,n) + sg12(i,k,j)=sigmag_aer(m,n) + hg12(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 1.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn13(i,k,j)=fn(m,n) + na13(i,k,j)=naerosol(m,n) + vo13(i,k,j)=vaerosol(m,n) + dl13(i,k,j)=dlo_sect(m,n) + dh13(i,k,j)=dhi_sect(m,n) + sg13(i,k,j)=sigmag_aer(m,n) + hg13(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 1.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn14(i,k,j)=fn(m,n) + na14(i,k,j)=naerosol(m,n) + vo14(i,k,j)=vaerosol(m,n) + dl14(i,k,j)=dlo_sect(m,n) + dh14(i,k,j)=dhi_sect(m,n) + sg14(i,k,j)=sigmag_aer(m,n) + hg14(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn21(i,k,j)=fn(m,n) + na21(i,k,j)=naerosol(m,n) + vo21(i,k,j)=vaerosol(m,n) + dl21(i,k,j)=dlo_sect(m,n) + dh21(i,k,j)=dhi_sect(m,n) + sg21(i,k,j)=sigmag_aer(m,n) + hg21(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn22(i,k,j)=fn(m,n) + na22(i,k,j)=naerosol(m,n) + vo22(i,k,j)=vaerosol(m,n) + dl22(i,k,j)=dlo_sect(m,n) + dh22(i,k,j)=dhi_sect(m,n) + sg22(i,k,j)=sigmag_aer(m,n) + hg22(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn23(i,k,j)=fn(m,n) + na23(i,k,j)=naerosol(m,n) + vo23(i,k,j)=vaerosol(m,n) + dl23(i,k,j)=dlo_sect(m,n) + dh23(i,k,j)=dhi_sect(m,n) + sg23(i,k,j)=sigmag_aer(m,n) + hg23(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 2.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn24(i,k,j)=fn(m,n) + na24(i,k,j)=naerosol(m,n) + vo24(i,k,j)=vaerosol(m,n) + dl24(i,k,j)=dlo_sect(m,n) + dh24(i,k,j)=dhi_sect(m,n) + sg24(i,k,j)=sigmag_aer(m,n) + hg24(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn31(i,k,j)=fn(m,n) + na31(i,k,j)=naerosol(m,n) + vo31(i,k,j)=vaerosol(m,n) + dl31(i,k,j)=dlo_sect(m,n) + dh31(i,k,j)=dhi_sect(m,n) + sg31(i,k,j)=sigmag_aer(m,n) + hg31(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn32(i,k,j)=fn(m,n) + na32(i,k,j)=naerosol(m,n) + vo32(i,k,j)=vaerosol(m,n) + dl32(i,k,j)=dlo_sect(m,n) + dh32(i,k,j)=dhi_sect(m,n) + sg32(i,k,j)=sigmag_aer(m,n) + hg32(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn33(i,k,j)=fn(m,n) + na33(i,k,j)=naerosol(m,n) + vo33(i,k,j)=vaerosol(m,n) + dl33(i,k,j)=dlo_sect(m,n) + dh33(i,k,j)=dhi_sect(m,n) + sg33(i,k,j)=sigmag_aer(m,n) + hg33(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 3.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn34(i,k,j)=fn(m,n) + na34(i,k,j)=naerosol(m,n) + vo34(i,k,j)=vaerosol(m,n) + dl34(i,k,j)=dlo_sect(m,n) + dh34(i,k,j)=dhi_sect(m,n) + sg34(i,k,j)=sigmag_aer(m,n) + hg34(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 1.0) < eps)) then + fn41(i,k,j)=fn(m,n) + na41(i,k,j)=naerosol(m,n) + vo41(i,k,j)=vaerosol(m,n) + dl41(i,k,j)=dlo_sect(m,n) + dh41(i,k,j)=dhi_sect(m,n) + sg41(i,k,j)=sigmag_aer(m,n) + hg41(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 2.0) < eps)) then + fn42(i,k,j)=fn(m,n) + na42(i,k,j)=naerosol(m,n) + vo42(i,k,j)=vaerosol(m,n) + dl42(i,k,j)=dlo_sect(m,n) + dh42(i,k,j)=dhi_sect(m,n) + sg42(i,k,j)=sigmag_aer(m,n) + hg42(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 3.0) < eps)) then + fn43(i,k,j)=fn(m,n) + na43(i,k,j)=naerosol(m,n) + vo43(i,k,j)=vaerosol(m,n) + dl43(i,k,j)=dlo_sect(m,n) + dh43(i,k,j)=dhi_sect(m,n) + sg43(i,k,j)=sigmag_aer(m,n) + hg43(i,k,j)=hygro_aer(m,n) + endif + if ((abs(m - 4.0) < eps) .and. (abs(n - 4.0) < eps)) then + fn44(i,k,j)=fn(m,n) + na44(i,k,j)=naerosol(m,n) + vo44(i,k,j)=vaerosol(m,n) + dl44(i,k,j)=dlo_sect(m,n) + dh44(i,k,j)=dhi_sect(m,n) + sg44(i,k,j)=sigmag_aer(m,n) + hg44(i,k,j)=hygro_aer(m,n) + endif + enddo + enddo + ! rce-comment ! the activation-fraction fluxes (fluxn, fluxm) from subr activate assume that @@ -2590,7 +2967,1535 @@ subroutine activate(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, & return end subroutine activate +!---------------------------------------------------------------------- +! 06-nov-2005 rce - grid_id & ktau added to arg list + subroutine activate_ml(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, qv, & + msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, & + na, volc, dlo_sect,dhi_sect,sigman, hygro, & + fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, & + grid_id, ktau, ii, jj, kk,smax_prescribed )!BSINGH - Added smax_prescribed for WRFCuP + +! calculates number, surface, and mass fraction of aerosols activated as CCN +! calculates flux of cloud droplets, surface area, and aerosol mass into cloud +! assumes an internal mixture within each of aerosol mode. +! A sectional treatment within each type is assumed if ntype_aer >7. +! A gaussiam spectrum of updrafts can be treated. + +! mks units + +! Abdul-Razzak and Ghan, A parameterization of aerosol activation. +! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844. + use module_wrf_error , only : wrf_err_message + USE module_model_constants, only: g,rhowater, xlv, cp, rvovrd, r_d, r_v, & + mwdry,svp1,svp2,svp3,ep_2 + + ! Liran add ML start ==================================== + use, intrinsic :: iso_c_binding + use, intrinsic :: iso_fortran_env, only: real32 + use ftorch, only: & + torch_model, torch_tensor, & + torch_model_load, torch_model_forward, & + torch_kCPU, torch_kCUDA, & + torch_tensor_from_array_real32_2d_default_layout ! generic present in your build + ! Liran add ML end ==================================== + + implicit none + + +! input. + + integer,intent(in) :: maxd_atype ! dimension of types + integer,intent(in) :: maxd_asize ! dimension of sizes + integer,intent(in) :: ntype_aer ! number of types + integer,intent(in) :: nsize_aer(maxd_atype) ! number of sizes for type + integer,intent(in) :: msectional ! 1 for sectional, 0 for modal + integer,intent(in) :: grid_id ! WRF grid%id + integer,intent(in) :: ktau ! WRF time step count + integer,intent(in) :: ii, jj, kk ! i,j,k of current grid cell + real,intent(in) :: wbar ! grid cell mean vertical velocity (m/s) + real,intent(in) :: sigw ! subgrid standard deviation of vertical vel (m/s) + real,intent(in) :: wdiab ! diabatic vertical velocity (0 if adiabatic) + real,intent(in) :: wminf ! minimum updraft velocity for integration (m/s) + real,intent(in) :: wmaxf ! maximum updraft velocity for integration (m/s) + real,intent(in) :: tair ! air temperature (K) + real,intent(in) :: rhoair ! air density (kg/m3) + real,intent(in) :: qv ! water vapor mixing ratio (g/g) + real,intent(in) :: na(maxd_asize,maxd_atype) ! aerosol number concentration (/m3) + real,intent(in) :: sigman(maxd_asize,maxd_atype) ! geometric standard deviation of aerosol size distribution + real,intent(in) :: hygro(maxd_asize,maxd_atype) ! bulk hygroscopicity of aerosol mode + real,intent(in) :: volc(maxd_asize,maxd_atype) ! total aerosol volume concentration (m3/m3) + real,intent(in) :: dlo_sect( maxd_asize, maxd_atype ), & ! minimum size of section (cm) + dhi_sect( maxd_asize, maxd_atype ) ! maximum size of section (cm) + real,intent(in),optional :: smax_prescribed ! prescribed max. supersaturation for secondary activation !BSINGH - Added for WRFCuP + +! output + + real,intent(inout) :: fn(maxd_asize,maxd_atype) ! number fraction of aerosols activated + real,intent(inout) :: fs(maxd_asize,maxd_atype) ! surface fraction of aerosols activated + real,intent(inout) :: fm(maxd_asize,maxd_atype) ! mass fraction of aerosols activated + real,intent(inout) :: fluxn(maxd_asize,maxd_atype) ! flux of activated aerosol number fraction into cloud (m/s) + real,intent(inout) :: fluxs(maxd_asize,maxd_atype) ! flux of activated aerosol surface fraction (m/s) + real,intent(inout) :: fluxm(maxd_asize,maxd_atype) ! flux of activated aerosol mass fraction into cloud (m/s) + real,intent(inout) :: flux_fullact ! flux when activation fraction = 100% (m/s) + +! local + +!!$ external erf,erfc +!!$ real erf,erfc +! external qsat_water + integer, parameter:: nx=200 + integer iquasisect_option, isectional + real integ,integf + real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m) + real, parameter :: p0 = 1013.25e2 ! reference pressure (Pa) + real, parameter :: t0 = 273.15 ! reference temperature (K) + real ylo(maxd_asize,maxd_atype),yhi(maxd_asize,maxd_atype) ! 1-particle volume at section interfaces + real ymean(maxd_asize,maxd_atype) ! 1-particle volume at r=rmean + real ycut, lnycut, betayy, betayy2, gammayy, phiyy + real surfc(maxd_asize,maxd_atype) ! surface concentration (m2/m3) + real sign(maxd_asize,maxd_atype) ! geometric standard deviation of size distribution + real alnsign(maxd_asize,maxd_atype) ! natl log of geometric standard dev of aerosol + real am(maxd_asize,maxd_atype) ! number mode radius of dry aerosol (m) + real lnhygro(maxd_asize,maxd_atype) ! ln(b) + real f1(maxd_asize,maxd_atype) ! array to hold parameter for maxsat + real pres ! pressure (Pa) + real path ! mean free path (m) + real diff ! diffusivity (m2/s) + real conduct ! thermal conductivity (Joule/m/sec/deg) + real diff0,conduct0 + real es ! saturation vapor pressure + real qs ! water vapor saturation mixing ratio + real dqsdt ! change in qs with temperature + real dqsdp ! change in qs with pressure + real gg ! thermodynamic function (m2/s) + real sqrtg ! sqrt(gg) + real sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius + real lnsm(maxd_asize,maxd_atype) ! ln( sm ) + real zeta, eta(maxd_asize,maxd_atype) + real lnsmax ! ln(smax) + real alpha + real gamma + real beta + real gaus + logical :: top ! true if cloud top, false if cloud base or new cloud + real asub(maxd_asize,maxd_atype),bsub(maxd_asize,maxd_atype) ! coefficients of submode size distribution N=a+bx + real totn(maxd_atype) ! total aerosol number concentration + real aten ! surface tension parameter + real gmrad(maxd_atype) ! geometric mean radius + real gmradsq(maxd_atype) ! geometric mean of radius squared + real gmlnsig(maxd_atype) ! geometric standard deviation + real gmsm(maxd_atype) ! critical supersaturation at radius gmrad + real sumflxn(maxd_asize,maxd_atype) + real sumflxs(maxd_asize,maxd_atype) + real sumflxm(maxd_asize,maxd_atype) + real sumflx_fullact + real sumfn(maxd_asize,maxd_atype) + real sumfs(maxd_asize,maxd_atype) + real sumfm(maxd_asize,maxd_atype) + real sumns(maxd_atype) + real fnold(maxd_asize,maxd_atype) ! number fraction activated + real fsold(maxd_asize,maxd_atype) ! surface fraction activated + real fmold(maxd_asize,maxd_atype) ! mass fraction activated + real wold,gold + real alogten,alog2,alog3,alogaten + real alogam + real rlo(maxd_asize,maxd_atype), rhi(maxd_asize,maxd_atype) + real rmean(maxd_asize,maxd_atype) + ! mean radius (m) for the section (not used with modal) + ! calculated from current volume & number + real ccc + real dumaa,dumbb + real wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb + real dfmin,dfmax,fnew,fold,fnmin,fnbar,fsbar,fmbar + real alw,sqrtalw + real smax + real x,arg + real xmincoeff,xcut + real z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf + real etafactor1,etafactor2(maxd_asize,maxd_atype),etafactor2max + integer m,n,nw,nwmax + +! numerical integration parameters + real, parameter :: eps = 0.3 + real, parameter :: fmax = 0.99 + real, parameter :: sds = 3. + + ! Liran add ML start ==================================== + ! Kinds for FTorch interop + + ! FTorch variables + logical :: model_loaded = .false. + type(torch_model), save :: ml_model + type(torch_tensor) :: model_input_tensor, model_output_tensor + type(torch_tensor) :: tin, tout, in_t(1), out_t(1) ! arrays-of-tensors for forward() + integer(c_int) :: in_shape(2), out_shape(2) ! shapes for 1D arrays + + ! Dimensions (adjust as needed) + integer, parameter :: in_dim = 16 + integer, parameter :: out_dim = 4 + integer :: dev_idx + ! Host arrays (match real32 API; keep contiguous/target) + real(real32), allocatable :: input_array(:,:), output_array(:,:) + real(real32) :: tmp_swap + ! Misc + integer :: idx + character(len=256) :: msg + real :: a1_ml, b1_ml, AA_ml, BB_ml, CC_ml + real :: disc_ml, phi_ml, phi1_ml, phi2_ml + real :: phiyy_ml, betayy_ml, dumaa_fs_ml, dumbb_fs_ml, dumaa_fm_ml + real :: x_ml, rhs_ml, sigma_ml, fval_ml, dfval_ml + integer :: it_ml + ! Liran add ML end ==================================== + +! mathematical constants + real third, twothird, sixth, zero, one, two, three + + real, parameter :: sq2 = 1.4142135624 + real, parameter :: sqpi = 1.7724538509 + real, parameter :: pi = 3.1415926536 + +! integer, save :: ndist(nx) ! accumulates frequency distribution of integration bins required +! data ndist/nx*0/ + +! for nsize_aer>7, a sectional approach is used and isectional = iquasisect_option +! activation fractions (fn,fs,fm) are computed as follows +! iquasisect_option = 1,3 - each section treated as a narrow lognormal +! iquasisect_option = 2,4 - within-section dn/dx = a + b*x, x = ln(r) +! smax is computed as follows (when explicit activation is OFF) +! iquasisect_option = 1,2 - razzak-ghan modal parameterization with +! single mode having same ntot, dgnum, sigmag as the combined sections +! iquasisect_option = 3,4 - razzak-ghan sectional parameterization +! for nsize_aer=<9, a modal approach is used and isectional = 0 + +! rce 08-jul-2005 +! if either (na(n,m) < nsmall) or (volc(n,m) < vsmall) +! then treat bin/mode (n,m) as being empty, and set its fn/fs/fm=0.0 +! (for single precision, gradual underflow starts around 1.0e-38, +! and strange things can happen when in that region) + real, parameter :: nsmall = 1.0e-20 ! aer number conc in #/m3 + real, parameter :: vsmall = 1.0e-37 ! aer volume conc in m3/m3 + logical bin_is_empty(maxd_asize,maxd_atype), all_bins_empty + logical bin_is_narrow(maxd_asize,maxd_atype) + + integer idiagaa, ipass_nwloop + integer idiag_dndy_neg, idiag_fnsm_prob + +! The flag for cloud top is no longer used so set it to false. This is an +! antiquated feature related to radiation enhancing mass fluxes at cloud +! top. It is currently, as of Feb. 2009, set to false in the CAM version +! as well. + top = .false. + +!....................................................................... +! +! start calc. of modal or sectional activation properties (start of section 1) +! +!....................................................................... + idiag_dndy_neg = 1 ! set this to 0 to turn off + ! warnings about dn/dy < 0 + idiag_fnsm_prob = 1 ! set this to 0 to turn off + ! warnings about fn/fs/fm misbehavior + + iquasisect_option = 2 + if(msectional.gt.0)then + isectional = iquasisect_option + else + isectional = 0 + endif + !BSINGH - For WRFCuP + if ( present( smax_prescribed ) ) then + if (smax_prescribed <= 0.0) then + do n=1,ntype_aer + do m=1,nsize_aer(n) + fluxn(m,n)=0. + fn(m,n)=0. + fluxs(m,n)=0. + fs(m,n)=0. + fluxm(m,n)=0. + fm(m,n)=0. + end do + end do + flux_fullact=0. + return + end if + end if + + !BSINGH - ENDS + + do n=1,ntype_aer +! print *,'ntype_aer,n,nsize_aer(n)=',ntype_aer,n,nsize_aer(n) + + if(ntype_aer.eq.1.and.nsize_aer(n).eq.1.and.na(1,1).lt.1.e-20)then + fn(1,1)=0. + fs(1,1)=0. + fm(1,1)=0. + fluxn(1,1)=0. + fluxs(1,1)=0. + fluxm(1,1)=0. + flux_fullact=0. + return + endif + enddo + + zero = 0.0 + one = 1.0 + two = 2.0 + three = 3.0 + third = 1.0/3.0 + twothird = 2.0/3.0 !wig, 1-Mar-2009: Corrected value from 2/6 + sixth = 1.0/6.0 + + pres=r_d*rhoair*tair + diff0=0.211e-4*(p0/pres)*(tair/t0)**1.94 + conduct0=(5.69+0.017*(tair-t0))*4.186e2*1.e-5 ! convert to J/m/s/deg + es=1000.*svp1*exp( svp2*(tair-t0)/(tair-svp3) ) + qs=ep_2*es/(pres-es) + dqsdt=xlv/(r_v*tair*tair)*qs + alpha=g*(xlv/(cp*r_v*tair*tair)-1./(r_d*tair)) + gamma=(1+xlv/cp*dqsdt)/(rhoair*qs) + gg=1./(rhowater/(diff0*rhoair*qs)+xlv*rhowater/(conduct0*tair)*(xlv/(r_v*tair)-1.)) + sqrtg=sqrt(gg) + beta=4.*pi*rhowater*gg*gamma + aten=2.*surften/(r_v*tair*rhowater) + alogaten=log(aten) + alog2=log(two) + alog3=log(three) + ccc=4.*pi*third + etafactor2max=1.e10/(alpha*wmaxf)**1.5 ! this should make eta big if na is very small. + + all_bins_empty = .true. + do n=1,ntype_aer + totn(n)=0. + gmrad(n)=0. + gmradsq(n)=0. + sumns(n)=0. + do m=1,nsize_aer(n) + alnsign(m,n)=log(sigman(m,n)) +! internal mixture of aerosols + + bin_is_empty(m,n) = .true. + if (volc(m,n).gt.vsmall .and. na(m,n).gt.nsmall) then + bin_is_empty(m,n) = .false. + all_bins_empty = .false. + lnhygro(m,n)=log(hygro(m,n)) +! number mode radius (m,n) +! write(6,*)'alnsign,volc,na=',alnsign(m,n),volc(m,n),na(m,n) + am(m,n)=exp(-1.5*alnsign(m,n)*alnsign(m,n))* & + (3.*volc(m,n)/(4.*pi*na(m,n)))**third + + if (isectional .gt. 0) then +! sectional model. +! need to use bulk properties because parameterization doesn't +! work well for narrow bins. + totn(n)=totn(n)+na(m,n) + alogam=log(am(m,n)) + gmrad(n)=gmrad(n)+na(m,n)*alogam + gmradsq(n)=gmradsq(n)+na(m,n)*alogam*alogam + endif + etafactor2(m,n)=1./(na(m,n)*beta*sqrtg) + + if(hygro(m,n).gt.1.e-10)then + sm(m,n)=2.*aten/(3.*am(m,n))*sqrt(aten/(3.*hygro(m,n)*am(m,n))) + else + sm(m,n)=100. + endif +! write(6,*)'sm,hygro,am=',sm(m,n),hygro(m,n),am(m,n) + else + sm(m,n)=1. + etafactor2(m,n)=etafactor2max ! this should make eta big if na is very small. + + endif + lnsm(m,n)=log(sm(m,n)) + if ((isectional .eq. 3) .or. (isectional .eq. 4)) then + sumns(n)=sumns(n)+na(m,n)/sm(m,n)**twothird + endif +! write(6,'(a,i4,6g12.2)')'m,na,am,hygro,lnhygro,sm,lnsm=',m,na(m,n),am(m,n),hygro(m,n),lnhygro(m,n),sm(m,n),lnsm(m,n) + end do ! size + end do ! type + +! if all bins are empty, set all activation fractions to zero and exit + if ( all_bins_empty ) then + do n=1,ntype_aer + do m=1,nsize_aer(n) + fluxn(m,n)=0. + fn(m,n)=0. + fluxs(m,n)=0. + fs(m,n)=0. + fluxm(m,n)=0. + fm(m,n)=0. + end do + end do + flux_fullact=0. + return + endif + + + + if (isectional .le. 0) then + ! Initialize maxsat at this cell and timestep for the + ! modal setup (the sectional case is handled below). + call maxsat_init(maxd_atype, ntype_aer, & + maxd_asize, nsize_aer, alnsign, f1) + + goto 30000 + end if + + do n=1,ntype_aer + !wig 19-Oct-2006: Add zero trap based May 2006 e-mail from + !Ghan. Transport can clear out a cell leading to + !inconsistencies with the mass. + gmrad(n)=gmrad(n)/max(totn(n),1e-20) + gmlnsig=gmradsq(n)/totn(n)-gmrad(n)*gmrad(n) ! [ln(sigmag)]**2 + gmlnsig(n)=sqrt( max( 1.e-4, gmlnsig(n) ) ) + gmrad(n)=exp(gmrad(n)) + if ((isectional .eq. 3) .or. (isectional .eq. 4)) then + gmsm(n)=totn(n)/sumns(n) + gmsm(n)=gmsm(n)*gmsm(n)*gmsm(n) + gmsm(n)=sqrt(gmsm(n)) + else +! gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(1,n)*gmrad(n))) + gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(nsize_aer(n),n)*gmrad(n))) + endif + enddo + + ! Initialize maxsat at this cell and timestep for the + ! sectional setup (the modal case is handled above)... + call maxsat_init(maxd_atype, ntype_aer, & + maxd_asize, (/1/), gmlnsig, f1) + +!....................................................................... +! calculate sectional "sub-bin" size distribution +! +! dn/dy = nt*( a + b*y ) for ylo < y < yhi +! +! nt = na(m,n) = number mixing ratio of the bin +! y = v/vhi +! v = (4pi/3)*r**3 = particle volume +! vhi = v at r=rhi (upper bin boundary) +! ylo = y at lower bin boundary = vlo/vhi = (rlo/rhi)**3 +! yhi = y at upper bin boundary = 1.0 +! +! dv/dy = v * dn/dy = nt*vhi*( a*y + b*y*y ) +! +!....................................................................... +! 02-may-2006 - this dn/dy replaces the previous +! dn/dx = a + b*x where l = ln(r) +! the old dn/dx was overly complicated for cases of rmean near rlo or rhi +! the new dn/dy is consistent with that used in the movesect routine, +! which does continuous growth by condensation and aqueous chemistry +!....................................................................... + do 25002 n = 1,ntype_aer + do 25000 m = 1,nsize_aer(n) + +! convert from diameter in cm to radius in m + rlo(m,n) = 0.5*0.01*dlo_sect(m,n) + rhi(m,n) = 0.5*0.01*dhi_sect(m,n) + ylo(m,n) = (rlo(m,n)/rhi(m,n))**3 + yhi(m,n) = 1.0 + +! 04-nov-2005 - extremely narrow bins will be treated using 0/1 activation +! this is to avoid potential numerical problems + bin_is_narrow(m,n) = .false. + if ((rhi(m,n)/rlo(m,n)) .le. 1.01) bin_is_narrow(m,n) = .true. + +! rmean is mass mean radius for the bin; xmean = log(rmean) +! just use section midpoint if bin is empty + if ( bin_is_empty(m,n) ) then + rmean(m,n) = sqrt(rlo(m,n)*rhi(m,n)) + ymean(m,n) = (rmean(m,n)/rhi(m,n))**3 + goto 25000 + end if + + rmean(m,n) = (volc(m,n)/(ccc*na(m,n)))**third + rmean(m,n) = max( rlo(m,n), min( rhi(m,n), rmean(m,n) ) ) + ymean(m,n) = (rmean(m,n)/rhi(m,n))**3 + if ( bin_is_narrow(m,n) ) goto 25000 + +! if rmean is extremely close to either rlo or rhi, +! treat the bin as extremely narrow + if ((rhi(m,n)/rmean(m,n)) .le. 1.01) then + bin_is_narrow(m,n) = .true. + rlo(m,n) = min( rmean(m,n), (rhi(m,n)/1.01) ) + ylo(m,n) = (rlo(m,n)/rhi(m,n))**3 + goto 25000 + else if ((rmean(m,n)/rlo(m,n)) .le. 1.01) then + bin_is_narrow(m,n) = .true. + rhi(m,n) = max( rmean(m,n), (rlo(m,n)*1.01) ) + ylo(m,n) = (rlo(m,n)/rhi(m,n))**3 + ymean(m,n) = (rmean(m,n)/rhi(m,n))**3 + goto 25000 + endif + +! if rmean is somewhat close to either rlo or rhi, then dn/dy will be +! negative near the upper or lower bin boundary +! in these cases, assume that all the particles are in a subset of the full bin, +! and adjust rlo or rhi so that rmean will be near the center of this subset +! note that the bin is made narrower LOCALLY/TEMPORARILY, +! just for the purposes of the activation calculation + gammayy = (ymean(m,n)-ylo(m,n)) / (yhi(m,n)-ylo(m,n)) + if (gammayy .lt. 0.34) then + dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*(gammayy/0.34) + rhi(m,n) = rhi(m,n)*(dumaa**third) + ylo(m,n) = (rlo(m,n)/rhi(m,n))**3 + ymean(m,n) = (rmean(m,n)/rhi(m,n))**3 + else if (gammayy .ge. 0.66) then + dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*((gammayy-0.66)/0.34) + ylo(m,n) = dumaa + rlo(m,n) = rhi(m,n)*(dumaa**third) + end if + if ((rhi(m,n)/rlo(m,n)) .le. 1.01) then + bin_is_narrow(m,n) = .true. + goto 25000 + end if + + betayy = ylo(m,n)/yhi(m,n) + betayy2 = betayy*betayy + bsub(m,n) = (12.0*ymean(m,n) - 6.0*(1.0+betayy)) / & + (4.0*(1.0-betayy2*betayy) - 3.0*(1.0-betayy2)*(1.0+betayy)) + asub(m,n) = (1.0 - bsub(m,n)*(1.0-betayy2)*0.5) / (1.0-betayy) + + if ( asub(m,n)+bsub(m,n)*ylo(m,n) .lt. 0. ) then + if (idiag_dndy_neg .gt. 0) then + print *,'dndy<0 at lower boundary' + print *,'n,m=',n,m + print *,'na=',na(m,n),' volc=',volc(m,n) + print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc)) + print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n) + print *,'dlo_sect/2,dhi_sect/2=', & + (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n)) + print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n) + print *,'asub+bsub*ylo=', & + (asub(m,n)+bsub(m,n)*ylo(m,n)) + print *,'subr activate error 11 - i,j,k =', ii, jj, kk + endif + endif + if ( asub(m,n)+bsub(m,n)*yhi(m,n) .lt. 0. ) then + if (idiag_dndy_neg .gt. 0) then + print *,'dndy<0 at upper boundary' + print *,'n,m=',n,m + print *,'na=',na(m,n),' volc=',volc(m,n) + print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc)) + print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n) + print *,'dlo_sect/2,dhi_sect/2=', & + (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n)) + print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n) + print *,'asub+bsub*yhi=', & + (asub(m,n)+bsub(m,n)*yhi(m,n)) + print *,'subr activate error 12 - i,j,k =', ii, jj, kk + endif + endif + +25000 continue ! m=1,nsize_aer(n) +25002 continue ! n=1,ntype_aer + + +30000 continue +!....................................................................... +! +! end calc. of modal or sectional activation properties (end of section 1) +! +!....................................................................... + + + +! sjg 7-16-98 upward +! print *,'wbar,sigw=',wbar,sigw + + if(sigw.le.1.e-5) goto 50000 + +!....................................................................... +! +! start calc. of activation fractions/fluxes +! for spectrum of updrafts (start of section 2) +! +!....................................................................... + ipass_nwloop = 1 + idiagaa = 0 +! 06-nov-2005 rce - set idiagaa=1 for testing/debugging +! if ((grid_id.eq.1) .and. (ktau.eq.167) .and. & +! (ii.eq.24) .and. (jj.eq. 1) .and. (kk.eq.14)) idiagaa = 1 + +40000 continue + if(top)then + wmax=0. + wmin=min(zero,-wdiab) + else + wmax=min(wmaxf,wbar+sds*sigw) + wmin=max(wminf,-wdiab) + endif + wmin=max(wmin,wbar-sds*sigw) + w=wmin + dwmax=eps*sigw + dw=dwmax + dfmax=0.2 + dfmin=0.1 + if(wmax.le.w)then + do n=1,ntype_aer + do m=1,nsize_aer(n) + fluxn(m,n)=0. + fn(m,n)=0. + fluxs(m,n)=0. + fs(m,n)=0. + fluxm(m,n)=0. + fm(m,n)=0. + end do + end do + flux_fullact=0. + return + endif + do n=1,ntype_aer + do m=1,nsize_aer(n) + sumflxn(m,n)=0. + sumfn(m,n)=0. + fnold(m,n)=0. + sumflxs(m,n)=0. + sumfs(m,n)=0. + fsold(m,n)=0. + sumflxm(m,n)=0. + sumfm(m,n)=0. + fmold(m,n)=0. + enddo + enddo + sumflx_fullact=0. + + fold=0 + gold=0 +! 06-nov-2005 rce - set wold=w here +! wold=0 + wold=w + + +! 06-nov-2005 rce - define nwmax; calc dwmin from nwmax + nwmax = 200 +! dwmin = min( dwmax, 0.01 ) + dwmin = (wmax - wmin)/(nwmax-1) + dwmin = min( dwmax, dwmin ) + dwmin = max( 0.01, dwmin ) + +! +! loop over updrafts, incrementing sums as you go +! the "200" is (arbitrary) upper limit for number of updrafts +! if integration finishes before this, OK; otherwise, ERROR +! + if (idiagaa.gt.0) then + write(*,94700) ktau, grid_id, ii, jj, kk, nwmax + write(*,94710) 'wbar,sigw,wdiab=', wbar, sigw, wdiab + write(*,94710) 'wmin,wmax,dwmin,dwmax=', wmin, wmax, dwmin, dwmax + write(*,94720) -1, w, wold, dw + end if +94700 format( / 'activate 47000 - ktau,id,ii,jj,kk,nwmax=', 6i5 ) +94710 format( 'activate 47000 - ', a, 6(1x,f11.5) ) +94720 format( 'activate 47000 - nw,w,wold,dw=', i5, 3(1x,f11.5) ) + + do 47000 nw = 1, nwmax +41000 wnuc=w+wdiab + + if (idiagaa.gt.0) write(*,94720) nw, w, wold, dw + +! write(6,*)'wnuc=',wnuc + alw=alpha*wnuc + sqrtalw=sqrt(alw) + zeta=2.*sqrtalw*aten/(3.*sqrtg) + etafactor1=2.*alw*sqrtalw + if (isectional .gt. 0) then +! sectional model. +! use bulk properties + + do n=1,ntype_aer + if(totn(n).gt.1.e-10)then + eta(1,n)=etafactor1/(totn(n)*beta*sqrtg) + else + eta(1,n)=1.e10 + endif + enddo + !BSINGH - For WRFCuP scheme + ! use smax_prescribed if it is present; otherwise get smax from subr maxsat + if ( present( smax_prescribed ) ) then + smax = smax_prescribed + else + !BSINGH -ENDS + + call maxsat(zeta,eta,maxd_atype,ntype_aer, & + maxd_asize,(/1/),gmsm,gmlnsig,f1,smax) + endif + lnsmax=log(smax) + x=2*(log(gmsm(1))-lnsmax)/(3*sq2*gmlnsig(1)) + fnew=0.5*(1.-ERF_ALT(x)) + + else + + do n=1,ntype_aer + do m=1,nsize_aer(n) + eta(m,n)=etafactor1*etafactor2(m,n) + enddo + enddo + !BSINGH - For WRFCuP scheme + if ( present( smax_prescribed ) ) then + smax = smax_prescribed + else + !BSINGH -ENDS + call maxsat(zeta,eta,maxd_atype,ntype_aer, & + maxd_asize,nsize_aer,sm,alnsign,f1,smax) + endif +! write(6,*)'w,smax=',w,smax + + lnsmax=log(smax) + + x=2*(lnsm(nsize_aer(1),1)-lnsmax)/(3*sq2*alnsign(nsize_aer(1),1)) + fnew=0.5*(1.-ERF_ALT(x)) + + endif + + dwnew = dw +! 06-nov-2005 rce - "n" here should be "nw" (?) +! if(fnew-fold.gt.dfmax.and.n.gt.1)then + if(fnew-fold.gt.dfmax.and.nw.gt.1)then +! reduce updraft increment for greater accuracy in integration + if (dw .gt. 1.01*dwmin) then + dw=0.7*dw + dw=max(dw,dwmin) + w=wold+dw + go to 41000 + else + dwnew = dwmin + endif + endif + + if(fnew-fold.lt.dfmin)then +! increase updraft increment to accelerate integration + dwnew=min(1.5*dw,dwmax) + endif + fold=fnew + + z=(w-wbar)/(sigw*sq2) + gaus=exp(-z*z) + fnmin=1. + xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3 +! write(6,*)'xmincoeff=',xmincoeff + + + do 44002 n=1,ntype_aer + do 44000 m=1,nsize_aer(n) + if ( bin_is_empty(m,n) ) then + fn(m,n)=0. + fs(m,n)=0. + fm(m,n)=0. + else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then +! sectional +! within-section dn/dx = a + b*x + xcut=xmincoeff-third*lnhygro(m,n) +! ycut=(exp(xcut)/rhi(m,n))**3 +! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20 +! if (ycut > yhi), then actual value of ycut is unimportant, +! so do the following to avoid overflow + lnycut = 3.0 * ( xcut - log(rhi(m,n)) ) + lnycut = min( lnycut, log(yhi(m,n)*1.0e5) ) + ycut=exp(lnycut) +! write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n) +! if(lnsmax.lt.lnsmn(m,n))then + if(ycut.gt.yhi(m,n))then + fn(m,n)=0. + fs(m,n)=0. + fm(m,n)=0. + elseif(ycut.lt.ylo(m,n))then + fn(m,n)=1. + fs(m,n)=1. + fm(m,n)=1. + elseif ( bin_is_narrow(m,n) ) then +! 04-nov-2005 rce - for extremely narrow bins, +! do zero activation if xcut>xmean, 100% activation otherwise + if (ycut.gt.ymean(m,n)) then + fn(m,n)=0. + fs(m,n)=0. + fm(m,n)=0. + else + fn(m,n)=1. + fs(m,n)=1. + fm(m,n)=1. + endif + else + phiyy=ycut/yhi(m,n) + fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy) + if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then + if (idiag_fnsm_prob .gt. 0) then + print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21' + print *,'na,volc =', na(m,n), volc(m,n) + print *,'asub,bsub =', asub(m,n), bsub(m,n) + print *,'yhi,ycut =', yhi(m,n), ycut + endif + endif + + if (fn(m,n) .le. zero) then +! 10-nov-2005 rce - if fn=0, then fs & fm must be 0 + fn(m,n)=zero + fs(m,n)=zero + fm(m,n)=zero + else if (fn(m,n) .ge. one) then +! 10-nov-2005 rce - if fn=1, then fs & fm must be 1 + fn(m,n)=one + fs(m,n)=one + fm(m,n)=one + else +! 10-nov-2005 rce - otherwise, calc fm and check it + fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) + & + third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy)) + if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then + if (idiag_fnsm_prob .gt. 0) then + print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22' + print *,'na,volc,fn =', na(m,n), volc(m,n), fn(m,n) + print *,'asub,bsub =', asub(m,n), bsub(m,n) + print *,'yhi,ycut =', yhi(m,n), ycut + endif + endif + if (fm(m,n) .le. fn(m,n)) then +! 10-nov-2005 rce - if fm=fn, then fs must =fn + fm(m,n)=fn(m,n) + fs(m,n)=fn(m,n) + else if (fm(m,n) .ge. one) then +! 10-nov-2005 rce - if fm=1, then fs & fn must be 1 + fm(m,n)=one + fs(m,n)=one + fn(m,n)=one + else +! 10-nov-2005 rce - these two checks assure that the mean size +! of the activated & interstitial particles will be between rlo & rhi + dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n)) + fm(m,n) = min( fm(m,n), dumaa ) + dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n)) + fm(m,n) = min( fm(m,n), dumaa ) +! 10-nov-2005 rce - now calculate fs and bound it by fn, fm + betayy = ylo(m,n)/yhi(m,n) + dumaa = phiyy**twothird + dumbb = betayy**twothird + fs(m,n) = & + (asub(m,n)*(1.0-phiyy*dumaa) + & + 0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) / & + (asub(m,n)*(1.0-betayy*dumbb) + & + 0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb)) + fs(m,n)=max(fs(m,n),fn(m,n)) + fs(m,n)=min(fs(m,n),fm(m,n)) + endif + endif + endif + + else +! modal + x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n)) + fn(m,n)=0.5*(1.-ERF_ALT(x)) + arg=x-sq2*alnsign(m,n) + fs(m,n)=0.5*(1.-ERF_ALT(arg)) + arg=x-1.5*sq2*alnsign(m,n) + fm(m,n)=0.5*(1.-ERF_ALT(arg)) +! print *,'w,x,fn,fs,fm=',w,x,fn(m,n),fs(m,n),fm(m,n) + endif + +! fn(m,n)=1. !test +! fs(m,n)=1. +! fm(m,n)=1. + fnmin=min(fn(m,n),fnmin) +! integration is second order accurate +! assumes linear variation of f*gaus with w + wb=(w+wold) + fnbar=(fn(m,n)*gaus+fnold(m,n)*gold) + fsbar=(fs(m,n)*gaus+fsold(m,n)*gold) + fmbar=(fm(m,n)*gaus+fmold(m,n)*gold) + if((top.and.w.lt.0.).or.(.not.top.and.w.gt.0.))then + sumflxn(m,n)=sumflxn(m,n)+sixth*(wb*fnbar & + +(fn(m,n)*gaus*w+fnold(m,n)*gold*wold))*dw + sumflxs(m,n)=sumflxs(m,n)+sixth*(wb*fsbar & + +(fs(m,n)*gaus*w+fsold(m,n)*gold*wold))*dw + sumflxm(m,n)=sumflxm(m,n)+sixth*(wb*fmbar & + +(fm(m,n)*gaus*w+fmold(m,n)*gold*wold))*dw + endif + sumfn(m,n)=sumfn(m,n)+0.5*fnbar*dw +! write(6,'(a,9g10.2)')'lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw=', & +! lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw + fnold(m,n)=fn(m,n) + sumfs(m,n)=sumfs(m,n)+0.5*fsbar*dw + fsold(m,n)=fs(m,n) + sumfm(m,n)=sumfm(m,n)+0.5*fmbar*dw + fmold(m,n)=fm(m,n) + +44000 continue ! m=1,nsize_aer(n) +44002 continue ! n=1,ntype_aer + +! same form as sumflxm(m,n) but replace the fm/fmold(m,n) with 1.0 + sumflx_fullact = sumflx_fullact & + + sixth*(wb*(gaus+gold) + (gaus*w + gold*wold))*dw +! sumg=sumg+0.5*(gaus+gold)*dw + gold=gaus + wold=w + dw=dwnew + + if(nw.gt.1.and.(w.gt.wmax.or.fnmin.gt.fmax))go to 48000 + w=w+dw + +47000 continue ! nw = 1, nwmax + + + print *,'do loop is too short in activate' + print *,'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw + print *,'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab + print *,'wnuc=',wnuc + do n=1,ntype_aer + print *,'ntype=',n + print *,'na=',(na(m,n),m=1,nsize_aer(n)) + print *,'fn=',(fn(m,n),m=1,nsize_aer(n)) + end do +! dump all subr parameters to allow testing with standalone code +! (build a driver that will read input and call activate) + print *,'top,wbar,sigw,wdiab,tair,rhoair,ntype_aer=' + print *, top,wbar,sigw,wdiab,tair,rhoair,ntype_aer + print *,'na=' + print *, na + print *,'volc=' + print *, volc + print *,'sigman=' + print *, sigman + print *,'hygro=' + print *, hygro + + print *,'subr activate error 31 - i,j,k =', ii, jj, kk +! 06-nov-2005 rce - if integration fails, repeat it once with additional diagnostics + if (ipass_nwloop .eq. 1) then + ipass_nwloop = 2 + idiagaa = 2 + goto 40000 + end if + call wrf_error_fatal("STOP: activate before 48000") + +48000 continue + + +! ndist(n)=ndist(n)+1 + if(.not.top.and.w.lt.wmaxf)then + +! contribution from all updrafts stronger than wmax +! assuming constant f (close to fmax) + wnuc=w+wdiab + + z1=(w-wbar)/(sigw*sq2) + z2=(wmaxf-wbar)/(sigw*sq2) + integ=sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(z1)-ERFC_NUM_RECIPES(z2)) +! consider only upward flow into cloud base when estimating flux + wf1=max(w,zero) + zf1=(wf1-wbar)/(sigw*sq2) + gf1=exp(-zf1*zf1) + wf2=max(wmaxf,zero) + zf2=(wf2-wbar)/(sigw*sq2) + gf2=exp(-zf2*zf2) + gf=(gf1-gf2) + integf=wbar*sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(zf1)-ERFC_NUM_RECIPES(zf2))+sigw*sigw*gf + + do n=1,ntype_aer + do m=1,nsize_aer(n) + sumflxn(m,n)=sumflxn(m,n)+integf*fn(m,n) + sumfn(m,n)=sumfn(m,n)+fn(m,n)*integ + sumflxs(m,n)=sumflxs(m,n)+integf*fs(m,n) + sumfs(m,n)=sumfs(m,n)+fs(m,n)*integ + sumflxm(m,n)=sumflxm(m,n)+integf*fm(m,n) + sumfm(m,n)=sumfm(m,n)+fm(m,n)*integ + end do + end do +! same form as sumflxm(m,n) but replace the fm(m,n) with 1.0 + sumflx_fullact = sumflx_fullact + integf +! sumg=sumg+integ + endif + + + do n=1,ntype_aer + do m=1,nsize_aer(n) + +! fn(m,n)=sumfn(m,n)/(sumg) + fn(m,n)=sumfn(m,n)/(sq2*sqpi*sigw) + fluxn(m,n)=sumflxn(m,n)/(sq2*sqpi*sigw) + if(fn(m,n).gt.1.01)then + if (idiag_fnsm_prob .gt. 0) then + print *,'fn=',fn(m,n),' > 1 - activate err41' + print *,'w,m,n,na,am=',w,m,n,na(m,n),am(m,n) + print *,'integ,sumfn,sigw=',integ,sumfn(m,n),sigw + print *,'subr activate error - i,j,k =', ii, jj, kk + endif + fluxn(m,n) = fluxn(m,n)/fn(m,n) + endif + + fs(m,n)=sumfs(m,n)/(sq2*sqpi*sigw) + fluxs(m,n)=sumflxs(m,n)/(sq2*sqpi*sigw) + if(fs(m,n).gt.1.01)then + if (idiag_fnsm_prob .gt. 0) then + print *,'fs=',fs(m,n),' > 1 - activate err42' + print *,'m,n,isectional=',m,n,isectional + print *,'alnsign(m,n)=',alnsign(m,n) + print *,'rcut,rlo(m,n),rhi(m,n)',exp(xcut),rlo(m,n),rhi(m,n) + print *,'w,m,na,am=',w,m,na(m,n),am(m,n) + print *,'integ,sumfs,sigw=',integ,sumfs(m,n),sigw + endif + fluxs(m,n) = fluxs(m,n)/fs(m,n) + endif + +! fm(m,n)=sumfm(m,n)/(sumg) + fm(m,n)=sumfm(m,n)/(sq2*sqpi*sigw) + fluxm(m,n)=sumflxm(m,n)/(sq2*sqpi*sigw) + if(fm(m,n).gt.1.01)then + if (idiag_fnsm_prob .gt. 0) then + print *,'fm(',m,n,')=',fm(m,n),' > 1 - activate err43' + endif + fluxm(m,n) = fluxm(m,n)/fm(m,n) + endif + + end do + end do +! same form as fluxm(m,n) + flux_fullact = sumflx_fullact/(sq2*sqpi*sigw) + + goto 60000 +!....................................................................... +! +! end calc. of activation fractions/fluxes +! for spectrum of updrafts (end of section 2) +! +!....................................................................... + +!....................................................................... +! +! start calc. of activation fractions/fluxes +! for (single) uniform updraft (start of section 3) +! +!....................................................................... +50000 continue + + wnuc=wbar+wdiab +! write(6,*)'uniform updraft =',wnuc + +! 04-nov-2005 rce - moved the code for "wnuc.le.0" code to here + if(wnuc.le.0.)then + do n=1,ntype_aer + do m=1,nsize_aer(n) + fn(m,n)=0 + fluxn(m,n)=0 + fs(m,n)=0 + fluxs(m,n)=0 + fm(m,n)=0 + fluxm(m,n)=0 + end do + end do + flux_fullact=0. + return + endif + + w=wbar + alw=alpha*wnuc + sqrtalw=sqrt(alw) + zeta=2.*sqrtalw*aten/(3.*sqrtg) + + if (isectional .gt. 0) then +! sectional model. +! use bulk properties + do n=1,ntype_aer + if(totn(n).gt.1.e-10)then + eta(1,n)=2*alw*sqrtalw/(totn(n)*beta*sqrtg) + else + eta(1,n)=1.e10 + endif + end do + !BSINGH - For WRFCuP + ! use smax_prescribed if it is present; otherwise get smax from subr maxsat + if ( present( smax_prescribed ) ) then + smax = smax_prescribed + else + !BSINGH -ENDS + call maxsat(zeta,eta,maxd_atype,ntype_aer, & + maxd_asize,(/1/),gmsm,gmlnsig,f1,smax) + endif + + else + + do n=1,ntype_aer + do m=1,nsize_aer(n) + if(na(m,n).gt.1.e-10)then + eta(m,n)=2*alw*sqrtalw/(na(m,n)*beta*sqrtg) + else + eta(m,n)=1.e10 + endif + end do + end do + !BSINGH - For WRFCuP + ! use smax_prescribed if it is present; otherwise get smax from subr maxsat + if ( present( smax_prescribed ) ) then + smax = smax_prescribed + else + !BSINGH -ENDS + + call maxsat(zeta,eta,maxd_atype,ntype_aer, & + maxd_asize,nsize_aer,sm,alnsign,f1,smax) + endif + + endif + + lnsmax=log(smax) + xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3 + + do 55002 n=1,ntype_aer + do 55000 m=1,nsize_aer(n) + +! 04-nov-2005 rce - check for bin_is_empty here too, just like earlier + if ( bin_is_empty(m,n) ) then + fn(m,n)=0. + fs(m,n)=0. + fm(m,n)=0. + + else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then +! sectional +! within-section dn/dx = a + b*x + xcut=xmincoeff-third*lnhygro(m,n) +! ycut=(exp(xcut)/rhi(m,n))**3 +! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20 +! if (ycut > yhi), then actual value of ycut is unimportant, +! so do the following to avoid overflow + lnycut = 3.0 * ( xcut - log(rhi(m,n)) ) + lnycut = min( lnycut, log(yhi(m,n)*1.0e5) ) + ycut=exp(lnycut) +! write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n) +! if(lnsmax.lt.lnsmn(m,n))then + if(ycut.gt.yhi(m,n))then + fn(m,n)=0. + fs(m,n)=0. + fm(m,n)=0. +! elseif(lnsmax.gt.lnsmx(m,n))then + elseif(ycut.lt.ylo(m,n))then + fn(m,n)=1. + fs(m,n)=1. + fm(m,n)=1. + elseif ( bin_is_narrow(m,n) ) then +! 04-nov-2005 rce - for extremely narrow bins, +! do zero activation if xcut>xmean, 100% activation otherwise + if (ycut.gt.ymean(m,n)) then + fn(m,n)=0. + fs(m,n)=0. + fm(m,n)=0. + else + fn(m,n)=1. + fs(m,n)=1. + fm(m,n)=1. + endif + else + phiyy=ycut/yhi(m,n) + fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy) + if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then + if (idiag_fnsm_prob .gt. 0) then + print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21' + print *,'na,volc =', na(m,n), volc(m,n) + print *,'asub,bsub =', asub(m,n), bsub(m,n) + print *,'yhi,ycut =', yhi(m,n), ycut + endif + endif + + if (fn(m,n) .le. zero) then +! 10-nov-2005 rce - if fn=0, then fs & fm must be 0 + fn(m,n)=zero + fs(m,n)=zero + fm(m,n)=zero + else if (fn(m,n) .ge. one) then +! 10-nov-2005 rce - if fn=1, then fs & fm must be 1 + fn(m,n)=one + fs(m,n)=one + fm(m,n)=one + else +! 10-nov-2005 rce - otherwise, calc fm and check it + fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) + & + third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy)) + if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then + if (idiag_fnsm_prob .gt. 0) then + print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22' + print *,'na,volc,fn =', na(m,n), volc(m,n), fn(m,n) + print *,'asub,bsub =', asub(m,n), bsub(m,n) + print *,'yhi,ycut =', yhi(m,n), ycut + endif + endif + if (fm(m,n) .le. fn(m,n)) then +! 10-nov-2005 rce - if fm=fn, then fs must =fn + fm(m,n)=fn(m,n) + fs(m,n)=fn(m,n) + else if (fm(m,n) .ge. one) then +! 10-nov-2005 rce - if fm=1, then fs & fn must be 1 + fm(m,n)=one + fs(m,n)=one + fn(m,n)=one + else +! 10-nov-2005 rce - these two checks assure that the mean size +! of the activated & interstitial particles will be between rlo & rhi + dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n)) + fm(m,n) = min( fm(m,n), dumaa ) + dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n)) + fm(m,n) = min( fm(m,n), dumaa ) +! 10-nov-2005 rce - now calculate fs and bound it by fn, fm + betayy = ylo(m,n)/yhi(m,n) + dumaa = phiyy**twothird + dumbb = betayy**twothird + fs(m,n) = & + (asub(m,n)*(1.0-phiyy*dumaa) + & + 0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) / & + (asub(m,n)*(1.0-betayy*dumbb) + & + 0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb)) + fs(m,n)=max(fs(m,n),fn(m,n)) + fs(m,n)=min(fs(m,n),fm(m,n)) + endif + endif + + endif + + else +! modal + x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n)) + fn(m,n)=0.5*(1.-ERF_ALT(x)) + arg=x-sq2*alnsign(m,n) + fs(m,n)=0.5*(1.-ERF_ALT(arg)) + arg=x-1.5*sq2*alnsign(m,n) + fm(m,n)=0.5*(1.-ERF_ALT(arg)) + endif + +! fn(m,n)=1. ! test +! fs(m,n)=1. +! fm(m,n)=1. + if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then + fluxn(m,n)=fn(m,n)*w + fluxs(m,n)=fs(m,n)*w + fluxm(m,n)=fm(m,n)*w + else + fluxn(m,n)=0 + fluxs(m,n)=0 + fluxm(m,n)=0 + endif + +55000 continue ! m=1,nsize_aer(n) +55002 continue ! n=1,ntype_aer + + if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then + flux_fullact = w + else + flux_fullact = 0.0 + endif + +! 04-nov-2005 rce - moved the code for "wnuc.le.0" from here +! to near the start the uniform undraft section + +!....................................................................... +! +! end calc. of activation fractions/fluxes +! for (single) uniform updraft (end of section 3) +! +!....................................................................... + + + +60000 continue + + + +! Liran add ML start ==================================== +! Allocate once +if (.not. allocated(input_array)) allocate(input_array(1, in_dim)) +if (.not. allocated(output_array)) allocate(output_array(1, out_dim)) +input_array = 0.0_real32 +output_array = 0.0_real32 + +! Try to load the model + if (.not. model_loaded) then + ! Load model on CPU (use enum, not string) + CALL wrf_debug(15,'Load model on CPU 00') + call torch_model_load( ml_model, & + '/pscratch/sd/y/yanxia/PINN_Parcel_Model/artifacts/deployment/mlp_wrapper.pt', & + device_type = torch_kCPU, device_index = 0 ) + model_loaded = .true. + CALL wrf_debug(15,'ML activation model loaded successfully') + else + idx = 1 + input_array(1,idx) = tair; idx = idx + 1 + input_array(1,idx) = pres/100; idx = idx + 1 + input_array(1,idx) = qv/qs*100; idx = idx + 1 + input_array(1,idx) = wbar*100; idx = idx + 1 + do n=1,2; do m=1,2; input_array(1,idx) = na(m,n)*1e-6; idx = idx + 1; end do; end do + do n=1,2; do m=1,2 + if (na(m,n) > 0.0) then + !input_array(1,idx) = real( ((volc(m,n)/(ccc*na(m,n)))**third)*1e6, kind=real32 ) + input_array(1,idx) = real(exp(-1.5*alnsign(m,n)*alnsign(m,n))* (3.*volc(m,n)/(4.*pi*na(m,n)))**third*1e6, kind=real32 ) + else + input_array(1,idx) = 0.0_real32 + end if + idx = idx + 1 + end do; end do + + do n=1,2; do m=1,2; input_array(1,idx) = hygro(m,n); idx = idx + 1; end do; end do + idx = 1 + do n=1,2; do m=1,2; output_array(1,idx) = 0.0; idx = idx + 1; end do; end do + + ! 9 <-> 10 + tmp_swap = input_array(1,9) + input_array(1,9) = input_array(1,10) + input_array(1,10) = tmp_swap + + ! 5 <-> 6 + tmp_swap = input_array(1,5) + input_array(1,5) = input_array(1,6) + input_array(1,6) = tmp_swap + + ! 13 <-> 14 + tmp_swap = input_array(1,13) + input_array(1,13) = input_array(1,14) + input_array(1,14) = tmp_swap + + ! Rank-2 shapes (features, batch) + in_shape = [ int(in_dim, c_int), 1_c_int ] + out_shape = [ int(out_dim, c_int), 1_c_int ] + + WRITE( msg,* ) 'tair: ', tair + CALL wrf_debug ( 15 , TRIM(msg) ) + WRITE( msg,* ) 'pres: ', pres + CALL wrf_debug ( 15 , TRIM(msg) ) + WRITE( msg,* ) 'qv/qs: ', qv/qs + CALL wrf_debug ( 15 , TRIM(msg) ) + WRITE( msg,* ) 'wbar: ', wbar + CALL wrf_debug ( 15 , TRIM(msg) ) + + dev_idx = int(0, kind=c_int) + ! Input/Output tensor from host buffer (CPU) + call torch_tensor_from_array_real32_2d_default_layout( & + tin, input_array, torch_kCPU, requires_grad=.false. ) + + CALL wrf_debug(15,'tin finished') + ! Forward expects arrays-of-tensors, even for single in/out + in_t(1) = tin + WRITE(msg,*) 'tin shape = ', in_shape(1), ' x ', in_shape(2) + CALL wrf_debug(15, TRIM(msg)) + ! e.g., a couple of input features + WRITE(msg,*) 'input_array(1,1)=tair (K)=', input_array(1,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(2,1)=pres (hPa)=', input_array(2,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(3,1)=RH (%)=', input_array(3,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(4,1)=wbar (cm/s)=', input_array(4,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(5,1)=na (#/cm3)', input_array(5,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(6,1)=na (#/cm3)', input_array(6,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(7,1)=na (#/cm3)', input_array(7,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(8,1)=na (#/cm3)', input_array(8,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(9,1)=rmean (um)', input_array(9,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(10,1)=rmean (um)', input_array(10,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(11,1)=rmean (um)', input_array(11,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(12,1)=rmean (um)', input_array(12,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(13,1)=hygro', input_array(13,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(14,1)=hygro', input_array(14,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(15,1)=hygro', input_array(15,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'input_array(16,1)=hygro', input_array(16,1) + CALL wrf_debug(15, TRIM(msg)) + call torch_tensor_from_array_real32_2d_default_layout( & + tout, output_array, torch_kCPU, requires_grad=.false. ) + CALL wrf_debug(15,'tout finished') + out_t(1) = tout + CALL wrf_debug(15,'tout finished2') + CALL wrf_debug(15,'torch_model_forward start') + call torch_model_forward(ml_model, in_t, out_t) + CALL wrf_debug(15,'torch_model_forward finished') + ! ---------------- print emulator outputs ---------------- + do n = 1, out_dim + write(msg,'(A,I0,A,ES12.5)') 'ML output_array(1,', n, ')=', output_array(1,n) + call wrf_debug(15, trim(msg)) + end do + WRITE(msg,*) 'fn(1,1)=', fn(1,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'fn(2,1)=', fn(2,1) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'fn(1,2)=', fn(1,2) + CALL wrf_debug(15, TRIM(msg)) + WRITE(msg,*) 'fn(2,2)=', fn(2,2) + CALL wrf_debug(15, TRIM(msg)) + + ! ===== Map emulator outputs -> fn(m,n) in bin order (idx = 1..out_dim) ===== + idx = 1 + do n = 1, 2 + do m = 1, 2 + fn(m,n) = max( 0.0, min( 1.0, real(output_array(1,idx)) ) ) + idx = idx + 1 + end do + end do + + ! ===== Recompute fs, fm consistently with the new fn ===== + + if ( (isectional == 2) .or. (isectional == 4) ) then + ! ---------------- SECTIONAL branch ---------------- + do n = 1, ntype_aer + do m = 1, nsize_aer(n) + + if ( bin_is_empty(m,n) ) then + fn(m,n) = 0.0 + fs(m,n) = 0.0 + fm(m,n) = 0.0 + cycle + end if + + if (fn(m,n) <= 0.0) then + fn(m,n) = 0.0 + fs(m,n) = 0.0 + fm(m,n) = 0.0 + cycle + else if (fn(m,n) >= 1.0) then + fn(m,n) = 1.0 + fs(m,n) = 1.0 + fm(m,n) = 1.0 + cycle + end if + + ! Recover phi_ml from: fn = a*(1-phi_ml) + 0.5*b*(1-phi_ml^2) + a1_ml = asub(m,n) + b1_ml = bsub(m,n) + AA_ml = 0.5*b1_ml + BB_ml = a1_ml + CC_ml = fn(m,n) - a1_ml - 0.5*b1_ml + + if (abs(AA_ml) < 1.0e-30) then + ! Degenerate: fn ≈ a*(1-phi_ml) => phi_ml ≈ 1 - fn/a + if (a1_ml > 1.0e-30) then + phi_ml = 1.0 - fn(m,n)/a1_ml + else + phi_ml = 1.0 + end if + else + disc_ml = BB_ml*BB_ml - 4.0*AA_ml*CC_ml + if (disc_ml < 0.0) disc_ml = 0.0 + phi1_ml = (-BB_ml + sqrt(disc_ml)) / (2.0*AA_ml) + phi2_ml = (-BB_ml - sqrt(disc_ml)) / (2.0*AA_ml) + phi_ml = 0.5 + if ( (phi1_ml >= 0.0) .and. (phi1_ml <= 1.0) ) then + phi_ml = phi1_ml + else if ( (phi2_ml >= 0.0) .and. (phi2_ml <= 1.0) ) then + phi_ml = phi2_ml + end if + end if + if (phi_ml < 0.0) phi_ml = 0.0 + if (phi_ml > 1.0) phi_ml = 1.0 + + ! Recompute fm (mass fraction) + phiyy_ml = phi_ml + fm(m,n) = ( yhi(m,n)/ymean(m,n) ) * ( 0.5*asub(m,n)*(1.0 - phiyy_ml*phiyy_ml) + & + third*bsub(m,n)*(1.0 - phiyy_ml*phiyy_ml*phiyy_ml) ) + + if (fm(m,n) < fn(m,n)) fm(m,n) = fn(m,n) + if (fm(m,n) > 1.0) fm(m,n) = 1.0 + + dumaa_fm_ml = fn(m,n)*( yhi(m,n)/ymean(m,n) ) + if (fm(m,n) > dumaa_fm_ml) fm(m,n) = dumaa_fm_ml + dumaa_fm_ml = 1.0 + (fn(m,n)-1.0)*( ylo(m,n)/ymean(m,n) ) + if (fm(m,n) > dumaa_fm_ml) fm(m,n) = dumaa_fm_ml + + ! Recompute fs (surface fraction) + betayy_ml = ylo(m,n)/yhi(m,n) + dumaa_fs_ml = phiyy_ml**twothird + dumbb_fs_ml = betayy_ml**twothird + + fs(m,n) = ( asub(m,n)*(1.0 - phiyy_ml*dumaa_fs_ml) + 0.625*bsub(m,n)*(1.0 - phiyy_ml*phiyy_ml*dumaa_fs_ml) ) / & + ( asub(m,n)*(1.0 - betayy_ml*dumbb_fs_ml) + 0.625*bsub(m,n)*(1.0 - betayy_ml*betayy_ml*dumbb_fs_ml) ) + + if (fs(m,n) < fn(m,n)) fs(m,n) = fn(m,n) + if (fs(m,n) > fm(m,n)) fs(m,n) = fm(m,n) + + end do + end do + + else + ! ---------------- MODAL branch ---------------- + ! Invert fn = 0.5*(1 - erf(x)) via a few Newton steps on erf(x) = 1 - 2*fn + do n = 1, ntype_aer + do m = 1, nsize_aer(n) + + if ( bin_is_empty(m,n) ) then + fn(m,n) = 0.0 + fs(m,n) = 0.0 + fm(m,n) = 0.0 + cycle + end if + + if (fn(m,n) <= 0.0) then + fn(m,n) = 0.0 + fs(m,n) = 0.0 + fm(m,n) = 0.0 + cycle + else if (fn(m,n) >= 1.0) then + fn(m,n) = 1.0 + fs(m,n) = 1.0 + fm(m,n) = 1.0 + cycle + end if + + rhs_ml = 1.0 - 2.0*fn(m,n) + + ! Winitzki-style initial guess for erfinv(rhs_ml) + if (rhs_ml >= 0.0) then + x = 1.0 + else + x = -1.0 + end if + ! A ~ 0.147 (constant in the approximation) + x = x * sqrt( sqrt( ( (2.0/(pi*0.147)) + 0.5*log(1.0 - rhs_ml*rhs_ml) )**2 - & + log(1.0 - rhs_ml*rhs_ml)/0.147 ) - & + ( (2.0/(pi*0.147)) + 0.5*log(1.0 - rhs_ml*rhs_ml) ) ) + + do it_ml = 1, 4 + fval_ml = erf(x) - rhs_ml + dfval_ml = 2.0/sqrt(pi) * exp( -x*x ) + x = x - fval_ml/dfval_ml + end do + + sigma_ml = alnsign(m,n) + fs(m,n) = 0.5 * ( 1.0 - erf( x - sqrt(2.0) * sigma_ml ) ) + fm(m,n) = 0.5 * ( 1.0 - erf( x - 1.5*sqrt(2.0) * sigma_ml ) ) + + if (fs(m,n) < fn(m,n)) fs(m,n) = fn(m,n) + if (fm(m,n) < fs(m,n)) fm(m,n) = fs(m,n) + if (fm(m,n) > 1.0) fm(m,n) = 1.0 + + end do + end do + end if + + ! ===== Refresh fluxes with the new fractions ===== + if ( (top .and. wbar < 0.0) .or. ( (.not. top) .and. wbar > 0.0 ) ) then + do n = 1, ntype_aer + do m = 1, nsize_aer(n) + fluxn(m,n) = fn(m,n) * wbar + fluxs(m,n) = fs(m,n) * wbar + fluxm(m,n) = fm(m,n) * wbar + end do + end do + flux_fullact = wbar + else + do n = 1, ntype_aer + do m = 1, nsize_aer(n) + fluxn(m,n) = 0.0 + fluxs(m,n) = 0.0 + fluxm(m,n) = 0.0 + end do + end do + flux_fullact = 0.0 + end if + end if + +! Liran add ML end ==================================== + + return + end subroutine activate_ml !---------------------------------------------------------------------- !----------------------------------------------------------------------