diff --git a/.gitmodules b/.gitmodules index b2d51bdfe..9f5dcf111 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,3 +2,7 @@ path = physics/Radiation/RRTMGP/rte-rrtmgp url = https://github.com/earth-system-radiation/rte-rrtmgp branch = main +[submodule "physics/MP/TEMPO"] + path = physics/MP/TEMPO/tempo + url = https://github.com/dustinswales/TEMPO.git + branch = feature/add_ccpp_type \ No newline at end of file diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 index eadc041fd..db36e6383 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 @@ -26,6 +26,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& ntdu1, ntdu2, ntdu3, ntdu4, ntdu5, ntss1, ntss2, & ntss3, ntss4, ntss5, ntsu, ntbcb, ntbcl, ntocb, ntocl, ntchm, & imp_physics,imp_physics_nssl, nssl_ccn_on, nssl_invertccn, & + imp_physics_tempo, & imp_physics_thompson, imp_physics_gfdl, imp_physics_zhao_carr, & imp_physics_zhao_carr_pdf, imp_physics_mg, imp_physics_wsm6, & imp_physics_fer_hires, iovr, iovr_rand, iovr_maxrand, iovr_max, & @@ -45,7 +46,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& gasvmr_ccl4, gasvmr_cfc113, aerodp,ext550, clouds6, clouds7, clouds8, & clouds9, cldsa, cldfra, cldfra2d, lwp_ex,iwp_ex, lwp_fc,iwp_fc, & faersw1, faersw2, faersw3, faerlw1, faerlw2, faerlw3, alpha, rrfs_sd, & - aero_dir_fdb, fdb_coef, spp_wts_rad, spp_rad, ico2, ozphys, & + aero_dir_fdb, fdb_coef, spp_wts_rad, spp_rad, ico2, ozphys, tempo_cfg, & errmsg, errflg) use machine, only: kind_phys @@ -72,15 +73,36 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& use surface_perturbation, only: cdfnor,ppfbet ! For Thompson MP - use module_mp_thompson, only: calc_effectRad, & - Nt_c_l, Nt_c_o, & - re_qc_min, re_qc_max, & - re_qi_min, re_qi_max, & - re_qs_min, re_qs_max + use module_mp_thompson, only: calc_effectRad_thmpsn => calc_effectRad, & + Nt_c_l_thmpsn => Nt_c_l, & + Nt_c_o_thmpsn => Nt_c_o, & + re_qc_min_thmpsn => re_qc_min, & + re_qc_max_thmpsn => re_qc_max, & + re_qi_min_thmpsn => re_qi_min, & + re_qi_max_thmpsn => re_qi_max, & + re_qs_min_thmpsn => re_qs_min, & + re_qs_max_thmpsn => re_qs_max use module_mp_thompson_make_number_concentrations, only: & - make_IceNumber, & - make_DropletNumber, & - make_RainNumber + make_IceNumber_thmpsn => make_IceNumber, & + make_DropletNumber_thmpsn => make_DropletNumber, & + make_RainNumber_thmpsn => make_RainNumber + + ! For TEMPO MP + ! DJS to Anders: Change from Thompson -> TEMPO in submodule, propogate change here. + use module_mp_thompson_params, only: Nt_c_l_tempo => Nt_c_l, & + Nt_c_o_tempo => Nt_c_o, & + re_qc_min_tempo => re_qc_min, & + re_qc_max_tempo => re_qc_max, & + re_qi_min_tempo => re_qi_min, & + re_qi_max_tempo => re_qi_max, & + re_qs_min_tempo => re_qs_min, & + re_qs_max_tempo => re_qs_max + use module_mp_thompson_utils, only: calc_effectRad_tempo => calc_effectRad, & + make_IceNumber_tempo => make_IceNumber, & + make_DropletNumber_tempo => make_DropletNumber, & + make_RainNumber_tempo => make_RainNumber + + ! For NRL Ozone use module_ozphys, only: ty_ozphys implicit none @@ -98,6 +120,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& lndp_type, & kdt, imp_physics, & imp_physics_thompson, & + imp_physics_tempo, & imp_physics_gfdl, & imp_physics_zhao_carr, & imp_physics_zhao_carr_pdf, & @@ -233,7 +256,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& qv_mp, qc_mp, qi_mp, qs_mp, & nc_mp, ni_mp, nwfa real (kind=kind_phys), dimension(lm) :: cldfra1d, qv1d, & - & qc1d, qi1d, qs1d, dz1d, p1d, t1d + qc1d, qi1d, qs1d, dz1d, p1d, t1d + ! For TEMPO MP + type(ty_tempo_cfg), intent(in) :: tempo_cfg ! for F-A MP real(kind=kind_phys), dimension(im,lm+LTP+1) :: tem2db, hz @@ -276,7 +301,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& LP1 = LM + 1 ! num of in/out levels - if (imp_physics == imp_physics_thompson) then + if (imp_physics == imp_physics_thompson .or. imp_physics == imp_physics_tempo) then max_relh = 1.5 else max_relh = 1.1 @@ -737,7 +762,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& enddo enddo ! for Thompson MP - prepare variables for calc_effr - if_thompson: if (imp_physics == imp_physics_thompson .and. (ltaerosol .or. mraerosol)) then + if_thompson: if ((imp_physics == imp_physics_thompson .or. imp_physics == imp_physics_tempo) .and. (ltaerosol .or. mraerosol)) then do k=1,LMK do i=1,IM qvs = qlyr(i,k) @@ -752,7 +777,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& nwfa (i,k) = tracer1(i,k,ntwa) enddo enddo - elseif (imp_physics == imp_physics_thompson) then + elseif (imp_physics == imp_physics_thompson .or. imp_physics == imp_physics_tempo) then do k=1,LMK do i=1,IM qvs = qlyr(i,k) @@ -763,9 +788,19 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) if(nint(slmsk(i)) == 1) then - nc_mp (i,k) = Nt_c_l*orho(i,k) + if (imp_physics == imp_physics_thompson) then + nc_mp (i,k) = Nt_c_l_thmpsn*orho(i,k) + endif + if (imp_physics == imp_physics_tempo) then + nc_mp (i,k) = Nt_c_l_tempo*orho(i,k) + endif else - nc_mp (i,k) = Nt_c_o*orho(i,k) + if (imp_physics == imp_physics_thompson) then + nc_mp (i,k) = Nt_c_o_thmpsn*orho(i,k) + endif + if (imp_physics == imp_physics_tempo) then + nc_mp (i,k) = Nt_c_o_tempo*orho(i,k) + endif endif ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs) enddo @@ -878,18 +913,28 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& ! not used yet -- effr_in should always be true for now endif - elseif (imp_physics == imp_physics_thompson) then ! Thompson MP + elseif (imp_physics == imp_physics_thompson .or. imp_physics == imp_physics_tempo) then ! Thompson MP ! ! Compute effective radii for QC, QI, QS with (GF, MYNN) or without (all others) sub-grid clouds ! ! Update number concentration, consistent with sub-grid clouds (GF, MYNN) or without (all others) do k=1,lm do i=1,im - if ((ltaerosol .or. mraerosol) .and. qc_mp(i,k)>1.e-12 .and. nc_mp(i,k)<100.) then - nc_mp(i,k) = make_DropletNumber(qc_mp(i,k)*rho(i,k), nwfa(i,k)*rho(i,k)) * orho(i,k) + if ((ltaerosol .or. mraerosol) .and. qc_mp(i,k)>1.e-12 .and. nc_mp(i,k)<100.) then + if (imp_physics == imp_physics_thompson) then + nc_mp(i,k) = make_DropletNumber_thmpsn(qc_mp(i,k)*rho(i,k), nwfa(i,k)*rho(i,k)) * orho(i,k) + endif + if (imp_physics == imp_physics_tempo) then + nc_mp(i,k) = make_DropletNumber_tempo(qc_mp(i,k)*rho(i,k), nwfa(i,k)*rho(i,k)) * orho(i,k) + endif endif if (qi_mp(i,k)>1.e-12 .and. ni_mp(i,k)<100.) then - ni_mp(i,k) = make_IceNumber(qi_mp(i,k)*rho(i,k), tlyr(i,k)) * orho(i,k) + if (imp_physics == imp_physics_thompson) then + ni_mp(i,k) = make_IceNumber_thmpsn(qi_mp(i,k)*rho(i,k), tlyr(i,k)) * orho(i,k) + endif + if (imp_physics == imp_physics_tempo) then + ni_mp(i,k) = make_IceNumber_tempo(qi_mp(i,k)*rho(i,k), tlyr(i,k)) * orho(i,k) + endif endif end do end do @@ -900,18 +945,38 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& !tgs: progclduni has different limits for ice radii (10.0-150.0) than ! calc_effectRad (4.99-125.0 for WRFv3.8.1; 2.49-125.0 for WRFv4+) ! it will raise the low limit from 5 to 10, but the high limit will remain 125. - call calc_effectRad (tlyr(i,:), plyr(i,:)*100., qv_mp(i,:), qc_mp(i,:), & - nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & - effrl(i,:), effri(i,:), effrs(i,:), islmsk, 1, lm ) - ! Scale Thompson's effective radii from meter to micron - do k=1,lm - effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max))*1.e6 - effri(i,k) = MAX(re_qi_min, MIN(effri(i,k), re_qi_max))*1.e6 - effrs(i,k) = MAX(re_qs_min, MIN(effrs(i,k), re_qs_max))*1.e6 - end do - effrl(i,lmk) = re_qc_min*1.e6 - effri(i,lmk) = re_qi_min*1.e6 - effrs(i,lmk) = re_qs_min*1.e6 + if (imp_physics == imp_physics_thompson) then + call calc_effectRad_thmpsn (tlyr(i,:), plyr(i,:)*100., qv_mp(i,:), qc_mp(i,:), & + nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & + effrl(i,:), effri(i,:), effrs(i,:), islmsk, 1, lm ) + endif + if (imp_physics == imp_physics_tempo) then + call calc_effectRad_tempo (t1d=tlyr(i,:), p1d=plyr(i,:)*100., qv1d=qv_mp(i,:), qc1d=qc_mp(i,:), & + nc1d=nc_mp(i,:), qi1d=qi_mp(i,:), ni1d=ni_mp(i,:), qs1d=qs_mp(i,:), & + re_qc1d=effrl(i,:), re_qi1d=effri(i,:), re_qs1d=effrs(i,:), kts=1, kte=lm, & + lsml=islmsk, configs=tempo_cfg) + end if + ! Scale Thompson/TEMPO's effective radii from meter to micron + if (imp_physics == imp_physics_thompson) then + do k=1,lm + effrl(i,k) = MAX(re_qc_min_thmpsn, MIN(effrl(i,k), re_qc_max_thmpsn))*1.e6 + effri(i,k) = MAX(re_qi_min_thmpsn, MIN(effri(i,k), re_qi_max_thmpsn))*1.e6 + effrs(i,k) = MAX(re_qs_min_thmpsn, MIN(effrs(i,k), re_qs_max_thmpsn))*1.e6 + end do + effrl(i,lmk) = re_qc_min_thmpsn*1.e6 + effri(i,lmk) = re_qi_min_thmpsn*1.e6 + effrs(i,lmk) = re_qs_min_thmpsn*1.e6 + end if + if (imp_physics == imp_physics_tempo) then + do k=1,lm + effrl(i,k) = MAX(re_qc_min_tempo, MIN(effrl(i,k), re_qc_max_tempo))*1.e6 + effri(i,k) = MAX(re_qi_min_tempo, MIN(effri(i,k), re_qi_max_tempo))*1.e6 + effrs(i,k) = MAX(re_qs_min_tempo, MIN(effrs(i,k), re_qs_max_tempo))*1.e6 + end do + effrl(i,lmk) = re_qc_min_tempo*1.e6 + effri(i,lmk) = re_qi_min_tempo*1.e6 + effrs(i,lmk) = re_qs_min_tempo*1.e6 + end if end do effrr(:,:) = 1000. ! rrain_def=1000. ! Update global arrays diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta index a7df9be2a..38af1361e 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta @@ -4,6 +4,7 @@ relative_path = ../../ dependencies = tools/funcphys.f90,hooks/machine.F dependencies = MP/Thompson/module_mp_thompson.F90,MP/Thompson/module_mp_thompson_make_number_concentrations.F90 + dependencies = MP/TEMPO/tempo/module_mp_thompson_params.F90 dependencies = Radiation/RRTMG/radcons.f90,Radiation/radiation_aerosols.f dependencies = Radiation/radiation_astronomy.f,Radiation/radiation_clouds.f,Radiation/radiation_gases.f dependencies = Radiation/RRTMG/radlw_param.f,Radiation/RRTMG/radsw_param.f,Radiation/radiation_cloud_overlap.F90 @@ -266,6 +267,13 @@ dimensions = () type = ty_ozphys intent = in +[tempo_cfg] + standard_name = configuration_for_TEMPO_microphysics + long_name = configuration information for TEMPO microphysics + units = mixed + dimensions = () + type = ty_tempo_cfg + intent = in [iaermdl] standard_name = control_for_aerosol_radiation_scheme long_name = control of aerosol scheme in radiation @@ -462,6 +470,13 @@ dimensions = () type = integer intent = in +[imp_physics_tempo] + standard_name = identifier_for_tempo_microphysics_scheme + long_name = choice of TEMPO microphysics scheme + units = flag + dimensions = () + type = integer + intent = in [imp_physics_thompson] standard_name = identifier_for_thompson_microphysics_scheme long_name = choice of Thompson microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.meta index 7f7ad7532..95f5c1c92 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.meta @@ -6,7 +6,7 @@ dependencies = Radiation/radiation_aerosols.f dependencies = Radiation/radiation_astronomy.f,Radiation/radiation_clouds.f,Radiation/radiation_gases.f dependencies = Radiation/RRTMG/radlw_main.F90,Radiation/RRTMG/radlw_param.f,Radiation/RRTMG/radsw_main.F90,Radiation/RRTMG/radsw_param.f - dependencies = MP/Thompson/module_mp_thompson.F90,photochem/module_ozphys.F90 + dependencies = photochem/module_ozphys.F90 ######################################################################## [ccpp-arg-table] diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 index f61ea76c6..e18b8d5ac 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 @@ -12,10 +12,11 @@ module GFS_rrtmgp_cloud_mp use rrtmgp_lw_cloud_optics, only: & radliq_lwr => radliq_lwrLW, radliq_upr => radliq_uprLW,& radice_lwr => radice_lwrLW, radice_upr => radice_uprLW - use module_mp_thompson, only: calc_effectRad, Nt_c_l, Nt_c_o, re_qc_min, re_qc_max, & - re_qi_min, re_qi_max, re_qs_min, re_qs_max - use module_mp_thompson_make_number_concentrations, only: make_IceNumber, & + use module_mp_thompson_utils, only: calc_effectRad, make_IceNumber, & make_DropletNumber, make_RainNumber + use module_mp_thompson_params, only: Nt_c_l, Nt_c_o, re_qc_min, re_qc_max, & + re_qi_min, re_qi_max, re_qs_min, re_qs_max, configs + real (kind_phys), parameter :: & cld_limit_lower = 0.001, & @@ -901,10 +902,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice ! Compute effective radii for liquid/ice/snow. do iCol=1,nCol ilsmask = nint(lsmask(iCol)) - call calc_effectRad (t_lay(iCol,:), p_lay(iCol,:), qv_mp(iCol,:), qc_mp(iCol,:), & - nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:), & - re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), ilsmask, & - 1, nLev ) + call calc_effectRad (t1d=t_lay(iCol,:), p1d=p_lay(iCol,:), qv1d=qv_mp(iCol,:), qc1d=qc_mp(iCol,:), & + nc1d=nc_mp(iCol,:), qi1d=qi_mp(iCol,:), ni1d=ni_mp(iCol,:), qs1d=qs_mp(iCol,:), & + re_qc1d=re_cloud(iCol,:), re_qi1d=re_ice(iCol,:), re_qs1d=re_snow(iCol,:), kts=1, kte=nLev, & + lsml=ilsmask, configs=configs) + do iLay = 1, nLev re_cloud(iCol,iLay) = MAX(re_qc_min, MIN(re_cloud(iCol,iLay), re_qc_max)) re_ice(iCol,iLay) = MAX(re_qi_min, MIN(re_ice(iCol,iLay), re_qi_max)) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.meta index a0fb78333..7570490a4 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.meta @@ -4,7 +4,7 @@ relative_path = ../../ dependencies = hooks/machine.F dependencies = Radiation/radiation_tools.F90,Radiation/radiation_clouds.f,Radiation/RRTMGP/rrtmgp_lw_cloud_optics.F90 - dependencies = MP/Thompson/module_mp_thompson_make_number_concentrations.F90,MP/Thompson/module_mp_thompson.F90 + dependencies = MP/TEMPO/tempo/module_mp_thompson_utils.F90,MP/TEMPO/tempo/module_mp_thompson_params.F90 ######################################################################## [ccpp-arg-table] diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta index fecb716ed..46802bea6 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta @@ -2,7 +2,7 @@ name = GFS_rrtmgp_setup type = scheme relative_path = ../../ - dependencies = hooks/machine.F,MP/Thompson/module_mp_thompson.F90 + dependencies = hooks/machine.F dependencies = Radiation/radiation_aerosols.f,photochem/module_ozphys.F90 dependencies = Radiation/radiation_gases.f,Radiation/radiation_astronomy.f diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 index d1aadb731..fed7b6c4a 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 @@ -10,14 +10,18 @@ module GFS_suite_interstitial_4 !! subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntrac, ntcw, ntiw, ntclamt, & ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, ntccn, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_nssl, nssl_invertccn, nssl_ccn_on, & + imp_physics_tempo, imp_physics_nssl, nssl_invertccn, nssl_ccn_on, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, dtidx, dtend,& index_of_process_conv_trans, gq0, clw, prsl, save_tcp, con_rd, con_eps, nssl_cccn, nwfa, spechum, ldiag3d, & qdiag3d, save_lnc, save_inc, ntk, ntke, otsptflag, errmsg, errflg) use machine, only: kind_phys - use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber - + use module_mp_thompson_make_number_concentrations, only: & + make_IceNumber_thmpsn => make_IceNumber, & + make_DropletNumber_thmpsn => make_DropletNumber + use module_mp_thompson_utils, only: & + make_IceNumber_tempo => make_IceNumber, & + make_DropletNumber_tempo => make_DropletNumber implicit none ! interface variables @@ -25,7 +29,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr logical, intent(in) :: otsptflag(:)! on/off switch for tracer transport by updraft and integer, intent(in ) :: im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, ntrw, & ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, ntccn, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imp_physics_nssl + imp_physics_tempo, imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imp_physics_nssl logical, intent(in) :: ltaerosol, convert_dry_rho logical, intent(in) :: nssl_ccn_on, nssl_invertccn @@ -225,7 +229,12 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + if (imp_physics == imp_physics_thompson) then + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber_thmpsn(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + endif + if (imp_physics == imp_physics_tempo) then + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber_tempo(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + endif !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -233,8 +242,13 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr !> - Convert moist mixing ratio to dry mixing ratio qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry - ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) + if (imp_physics == imp_physics_thompson) then + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber_thmpsn(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + endif + if (imp_physics == imp_physics_tempo) then + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber_tempo(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + endif !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif @@ -250,13 +264,23 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr !> - Update cloud water mixing ratio qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) !> - Update cloud water number concentration - gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + if (imp_physics == imp_physics_thompson) then + gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber_thmpsn(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + endif + if (imp_physics == imp_physics_tempo) then + gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber_tempo(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + endif endif if (ntinc>0) then !> - Update cloud ice mixing ratio qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) !> - Update cloud ice number concentration - gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + if (imp_physics == imp_physics_thompson) then + gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber_thmpsn(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + endif + if (imp_physics == imp_physics_tempo) then + gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber_tempo(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + endif endif enddo enddo diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta index c6a330ccb..7a248c564 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta @@ -4,6 +4,7 @@ type = scheme relative_path = ../../ dependencies = hooks/machine.F + dependencies = MP/TEMPO/tempo/module_mp_thompson_utils.F90 dependencies = MP/Thompson/module_mp_thompson_make_number_concentrations.F90 ######################################################################## @@ -157,6 +158,13 @@ dimensions = () type = integer intent = in +[imp_physics_tempo] + standard_name = identifier_for_tempo_microphysics_scheme + long_name = choice of TEMPO microphysics scheme + units = flag + dimensions = () + type = integer + intent = in [imp_physics_zhao_carr] standard_name = identifier_for_zhao_carr_microphysics_scheme long_name = choice of Zhao-Carr microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta index 8e25428cc..00fa68a7d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta @@ -4,7 +4,7 @@ relative_path = ../../ dependencies = tools/funcphys.f90,hooks/machine.F dependencies = hooks/physcons.F90,Radiation/RRTMG/radcons.f90 - dependencies = Radiation/radiation_clouds.f,MP/Thompson/module_mp_thompson.F90 + dependencies = Radiation/radiation_clouds.f ######################################################################## [ccpp-arg-table] diff --git a/physics/MP/TEMPO/module_mp_tempo.F90 b/physics/MP/TEMPO/module_mp_tempo.F90 new file mode 100644 index 000000000..f37e316a2 --- /dev/null +++ b/physics/MP/TEMPO/module_mp_tempo.F90 @@ -0,0 +1,1451 @@ +! 3D TEMPO Driver for CCPP +!================================================================================================================= +module module_mp_tempo + + use machine, only: wp => kind_phys, sp => kind_sngl_prec, dp => kind_dbl_prec + use module_mp_thompson_params + use module_mp_thompson_utils, only : create_bins, table_Efrw, table_Efsw, table_dropEvap, & + table_ccnAct, qi_aut_qs, qr_acr_qg_par, qr_acr_qs_par, freezeH2O_par, calc_refl10cm, calc_effectRad + use module_mp_thompson_main, only : mp_thompson_main + use module_mp_radar + + implicit none + +contains + !================================================================================================================= + ! This subroutine handles initialzation of the microphysics scheme including building of lookup tables, + ! allocating arrays for the microphysics scheme, and defining gamma function variables. + subroutine tempo_init(is_aerosol_aware_in, merra2_aerosol_aware_in, is_hail_aware_in, & + mpicomm, mpirank, mpiroot, threads, errmsg, errflg) + + logical, intent(in) :: is_aerosol_aware_in + logical, intent(in) :: merra2_aerosol_aware_in + logical, intent(in) :: is_hail_aware_in + type(MPI_Comm), intent(in) :: mpicomm + integer, intent(in) :: mpirank, mpiroot + integer, intent(in) :: threads + character(len=*), intent(inout) :: errmsg + integer, intent(inout) :: errflg + + integer :: i, j, k, l, m, n + logical :: micro_init + real(wp) :: stime, etime + logical, parameter :: precomputed_tables = .false. + + ! Set module variable is_aerosol_aware/merra2_aerosol_aware + configs%aerosol_aware = is_aerosol_aware_in + merra2_aerosol_aware = merra2_aerosol_aware_in + configs%hail_aware = is_hail_aware_in + if (configs%aerosol_aware .and. merra2_aerosol_aware) then + errmsg = 'Logic error in thompson_init: only one of the two options can be true, ' // & + 'not both: is_aerosol_aware or merra2_aerosol_aware' + errflg = 1 + return + end if + if (mpirank==mpiroot) then + if (configs%aerosol_aware) then + write (*,'(a)') 'Using aerosol-aware version of TEMPO microphysics' + else if(merra2_aerosol_aware) then + write (*,'(a)') 'Using merra2 aerosol-aware version of TEMPO microphysics' + else + write (*,'(a)') 'Using non-aerosol-aware version of TEMPO microphysics' + end if + end if + + micro_init = .false. + + if (configs%hail_aware) then + dimNRHG = NRHG + else + av_g(idx_bg1) = av_g_old + bv_g(idx_bg1) = bv_g_old + dimNRHG = NRHG1 + endif + + ! Allocate space for lookup tables (J. Michalakes 2009Jun08). + if (.not. allocated(tcg_racg)) then + allocate(tcg_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + micro_init = .true. + endif + + ! Rain-graupel (including array above tcg_racg) + if (.not. allocated(tmr_racg)) allocate(tmr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + if (.not. allocated(tcr_gacr)) allocate(tcr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + if (.not. allocated(tnr_racg)) allocate(tnr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + if (.not. allocated(tnr_gacr)) allocate(tnr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + + ! Rain-snow + if (.not. allocated(tcs_racs1)) allocate(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tmr_racs1)) allocate(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tcs_racs2)) allocate(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tmr_racs2)) allocate(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tcr_sacr1)) allocate(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tms_sacr1)) allocate(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tcr_sacr2)) allocate(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tms_sacr2)) allocate(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tnr_racs1)) allocate(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tnr_racs2)) allocate(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tnr_sacr1)) allocate(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.not. allocated(tnr_sacr2)) allocate(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) + + ! Cloud water freezing + if (.not. allocated(tpi_qcfz)) allocate(tpi_qcfz(ntb_c,nbc,ntb_t1,ntb_IN)) + if (.not. allocated(tni_qcfz)) allocate(tni_qcfz(ntb_c,nbc,ntb_t1,ntb_IN)) + + ! Rain freezing + if (.not. allocated(tpi_qrfz)) allocate(tpi_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_IN)) + if (.not. allocated(tpg_qrfz)) allocate(tpg_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_IN)) + if (.not. allocated(tni_qrfz)) allocate(tni_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_IN)) + if (.not. allocated(tnr_qrfz)) allocate(tnr_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_IN)) + + ! Ice growth and conversion to snow + if (.not. allocated(tps_iaus)) allocate(tps_iaus(ntb_i,ntb_i1)) + if (.not. allocated(tni_iaus)) allocate(tni_iaus(ntb_i,ntb_i1)) + if (.not. allocated(tpi_ide)) allocate(tpi_ide(ntb_i,ntb_i1)) + + ! Collision efficiencies + if (.not. allocated(t_efrw)) allocate(t_efrw(nbr,nbc)) + if (.not. allocated(t_efsw)) allocate(t_efsw(nbs,nbc)) + + ! Cloud water + if (.not. allocated(tnr_rev)) allocate(tnr_rev(nbr,ntb_r1,ntb_r)) + if (.not. allocated(tpc_wev)) allocate(tpc_wev(nbc,ntb_c,nbc)) + if (.not. allocated(tnc_wev)) allocate(tnc_wev(nbc,ntb_c,nbc)) + + ! CCN + if (.not. allocated(tnccn_act)) allocate(tnccn_act(ntb_arc,ntb_arw,ntb_art,ntb_arr,ntb_ark)) + + !================================================================================================================= + if_micro_init: if (micro_init) then + + !> - From Martin et al. (1994), assign gamma shape parameter mu for cloud + !! drops according to general dispersion characteristics (disp=~0.25 + !! for maritime and 0.45 for continental) + !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime + !.. to 2 for really dirty air. This not used in 2-moment cloud water + !.. scheme and nu_c used instead and varies from 2 to 15 (integer-only). + mu_c_l = min(15.0_wp, (1000.e6_wp/Nt_c_l + 2.)) + mu_c_o = min(15.0_wp, (1000.e6_wp/Nt_c_o + 2.)) + + !> - Compute Schmidt number to one-third used numerous times + Sc3 = Sc**(1./3.) + + !> - Compute minimum ice diam from mass, min snow/graupel mass from diam + D0i = (xm0i/am_i)**(1./bm_i) + xm0s = am_s * D0s**bm_s + xm0g = am_g(NRHG) * D0g**bm_g + + !> - Compute constants various exponents and gamma() associated with cloud, + !! rain, snow, and graupel + do n = 1, 15 + cce(1,n) = n + 1. + cce(2,n) = bm_r + n + 1. + cce(3,n) = bm_r + n + 4. + cce(4,n) = n + bv_c + 1. + cce(5,n) = bm_r + n + bv_c + 1. + ccg(1,n) = gamma(cce(1,n)) + ccg(2,n) = gamma(cce(2,n)) + ccg(3,n) = gamma(cce(3,n)) + ccg(4,n) = gamma(cce(4,n)) + ccg(5,n) = gamma(cce(5,n)) + ocg1(n) = 1.0 / ccg(1,n) + ocg2(n) = 1.0 / ccg(2,n) + enddo + + cie(1) = mu_i + 1. + cie(2) = bm_i + mu_i + 1. + cie(3) = bm_i + mu_i + bv_i + 1. + cie(4) = mu_i + bv_i + 1. + cie(5) = mu_i + 2. + cie(6) = bm_i*0.5 + mu_i + bv_i + 1. + cie(7) = bm_i*0.5 + mu_i + 1. + cig(1) = gamma(cie(1)) + cig(2) = gamma(cie(2)) + cig(3) = gamma(cie(3)) + cig(4) = gamma(cie(4)) + cig(5) = gamma(cie(5)) + cig(6) = gamma(cie(6)) + cig(7) = gamma(cie(7)) + oig1 = 1.0 / cig(1) + oig2 = 1.0 / cig(2) + obmi = 1.0 / bm_i + + cre(1) = bm_r + 1. + cre(2) = mu_r + 1. + cre(3) = bm_r + mu_r + 1. + cre(4) = bm_r*2. + mu_r + 1. + cre(5) = mu_r + bv_r + 1. + cre(6) = bm_r + mu_r + bv_r + 1. + cre(7) = bm_r*0.5 + mu_r + bv_r + 1. + cre(8) = bm_r + mu_r + bv_r + 3. + cre(9) = mu_r + bv_r + 3. + cre(10) = mu_r + 2. + cre(11) = 0.5*(bv_r + 5. + 2.*mu_r) + cre(12) = bm_r*0.5 + mu_r + 1. + cre(13) = bm_r*2. + mu_r + bv_r + 1. + + do n = 1, 13 + crg(n) = gamma(cre(n)) + enddo + + obmr = 1.0 / bm_r + ore1 = 1.0 / cre(1) + org1 = 1.0 / crg(1) + org2 = 1.0 / crg(2) + org3 = 1.0 / crg(3) + + cse(1) = bm_s + 1. + cse(2) = bm_s + 2. + cse(3) = bm_s*2. + cse(4) = bm_s + bv_s + 1. + cse(5) = bm_s*2. + bv_s + 1. + cse(6) = bm_s*2. + 1. + cse(7) = bm_s + mu_s + 1. + cse(8) = bm_s + mu_s + 2. + cse(9) = bm_s + mu_s + 3. + cse(10) = bm_s + mu_s + bv_s + 1. + cse(11) = bm_s*2. + mu_s + bv_s + 1. + cse(12) = bm_s*2. + mu_s + 1. + cse(13) = bv_s + 2. + cse(14) = bm_s + bv_s + cse(15) = mu_s + 1. + cse(16) = 1.0 + (1.0 + bv_s)/2. + + if (original_thompson) then + cse(17) = cse(16) + mu_s + 1. + cse(18) = bv_s + mu_s + 3. + do n = 1, 18 + csg(n) = gamma(cse(n)) + enddo + else + cse(17) = bm_s + bv_s + 2. + do n = 1, 17 + csg(n) = gamma(cse(n)) + enddo + endif + + oams = 1.0 / am_s + obms = 1.0 / bm_s + ocms = oams**obms + + cge(1,:) = bm_g + 1. + cge(2,:) = mu_g + 1. + cge(3,:) = bm_g + mu_g + 1. + cge(4,:) = bm_g*2. + mu_g + 1. + cge(10,:) = mu_g + 2. + cge(12,:) = bm_g*0.5 + mu_g + 1. + + do m = 1, NRHG + cge(5,m) = bm_g*2. + mu_g + bv_g(m) + 1. + cge(6,m) = bm_g + mu_g + bv_g(m) + 1. + cge(7,m) = bm_g*0.5 + mu_g + bv_g(m) + 1. + cge(8,m) = mu_g + bv_g(m) + 1. ! not used + cge(9,m) = mu_g + bv_g(m) + 3. + cge(11,m) = 0.5*(bv_g(m) + 5. + 2.*mu_g) + enddo + + do m = 1, NRHG + do n = 1, 12 + cgg(n,m) = gamma(cge(n,m)) + enddo + enddo + + oamg = 1.0 / am_g + obmg = 1.0 / bm_g + + do m = 1, NRHG + oamg(m) = 1.0 / am_g(m) + ocmg(m) = oamg(m)**obmg + enddo + + oge1 = 1.0 / cge(1,1) + ogg1 = 1.0 / cgg(1,1) + ogg2 = 1.0 / cgg(2,1) + ogg3 = 1.0 / cgg(3,1) + + !================================================================================================================= + ! Simplify various rate eqns the best we can now. + + ! Rain collecting cloud water and cloud ice + t1_qr_qc = PI * 0.25 * av_r * crg(9) + t1_qr_qi = PI * 0.25 * av_r * crg(9) + t2_qr_qi = PI * 0.25 * am_r*av_r * crg(8) + + ! Graupel collecting cloud water + ! t1_qg_qc = PI*.25*av_g * cgg(9) + + ! Snow collecting cloud water + t1_qs_qc = PI * 0.25 * av_s + + ! Snow collecting cloud ice + t1_qs_qi = PI * 0.25 * av_s + + ! Evaporation of rain; ignore depositional growth of rain. + t1_qr_ev = 0.78 * crg(10) + t2_qr_ev = 0.308 * Sc3 * SQRT(av_r) * crg(11) + + ! Sublimation/depositional growth of snow + t1_qs_sd = 0.86 + t2_qs_sd = 0.28 * Sc3 * SQRT(av_s) + + ! Melting of snow + t1_qs_me = PI * 4. *C_sqrd * olfus * 0.86 + t2_qs_me = PI * 4. *C_sqrd * olfus * 0.28 * Sc3 * SQRT(av_s) + + ! Sublimation/depositional growth of graupel + t1_qg_sd = 0.86 * cgg(10,1) + ! t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11) + + ! Melting of graupel + t1_qg_me = PI * 4. * C_cube * olfus * 0.86 * cgg(10,1) + ! t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11) + + + ! Constants for helping find lookup table indexes. + nic2 = nint(log10(r_c(1))) + nii2 = nint(log10(r_i(1))) + nii3 = nint(log10(Nt_i(1))) + nir2 = nint(log10(r_r(1))) + nir3 = nint(log10(N0r_exp(1))) + nis2 = nint(log10(r_s(1))) + nig2 = nint(log10(r_g(1))) + nig3 = nint(log10(N0g_exp(1))) + niIN2 = nint(log10(Nt_IN(1))) + + ! Create bins of cloud water (from minimum diameter to 100 microns). + Dc(1) = D0c*1.0_dp + dtc(1) = D0c*1.0_dp + do n = 2, nbc + Dc(n) = Dc(n-1) + 1.0e-6_dp + dtc(n) = (Dc(n) - Dc(n-1)) + enddo + + ! Create bins of cloud ice (from min diameter up to 2x min snow size). + call create_bins(numbins=nbi, lowbin=D0i*1.0_dp, highbin=D0s*2.0_dp, & + bins=Di, deltabins=dti) + + ! Create bins of rain (from min diameter up to 5 mm). + call create_bins(numbins=nbr, lowbin=D0r*1.0_dp, highbin=0.005_dp, & + bins=Dr, deltabins=dtr) + + ! Create bins of snow (from min diameter up to 2 cm). + call create_bins(numbins=nbs, lowbin=D0s*1.0_dp, highbin=0.02_dp, & + bins=Ds, deltabins=dts) + + ! Create bins of graupel (from min diameter up to 5 cm). + call create_bins(numbins=nbg, lowbin=D0g*1.0_dp, highbin=0.05_dp, & + bins=Dg, deltabins=dtg) + + ! Create bins of cloud droplet number concentration (1 to 3000 per cc). + call create_bins(numbins=nbc, lowbin=1.0_dp, highbin=3000.0_dp, & + bins=t_Nc) + t_Nc = t_Nc * 1.0e6_dp + nic1 = log(t_Nc(nbc)/t_Nc(1)) + + !================================================================================================================= + ! Create lookup tables for most costly calculations + + ! Assign mpicomm to module variable + mpi_communicator = mpicomm + + ! Standard tables are only written by master MPI task; + ! (physics init cannot be called by multiple threads, + ! hence no need to test for a specific thread number) + if (mpirank==mpiroot) then + thompson_table_writer = .true. + else + thompson_table_writer = .false. + end if + + precomputed_tables_1: if (.not.precomputed_tables) then + + call cpu_time(stime) + + do m = 1, ntb_r + do k = 1, ntb_r1 + do n = 1, dimNRHG + do j = 1, ntb_g + do i = 1, ntb_g1 + tcg_racg(i,j,n,k,m) = 0.0_dp + tmr_racg(i,j,n,k,m) = 0.0_dp + tcr_gacr(i,j,n,k,m) = 0.0_dp + tnr_racg(i,j,n,k,m) = 0.0_dp + tnr_gacr(i,j,n,k,m) = 0.0_dp + enddo + enddo + enddo + enddo + enddo + + do m = 1, ntb_r + do k = 1, ntb_r1 + do j = 1, ntb_t + do i = 1, ntb_s + tcs_racs1(i,j,k,m) = 0.0_dp + tmr_racs1(i,j,k,m) = 0.0_dp + tcs_racs2(i,j,k,m) = 0.0_dp + tmr_racs2(i,j,k,m) = 0.0_dp + tcr_sacr1(i,j,k,m) = 0.0_dp + tms_sacr1(i,j,k,m) = 0.0_dp + tcr_sacr2(i,j,k,m) = 0.0_dp + tms_sacr2(i,j,k,m) = 0.0_dp + tnr_racs1(i,j,k,m) = 0.0_dp + tnr_racs2(i,j,k,m) = 0.0_dp + tnr_sacr1(i,j,k,m) = 0.0_dp + tnr_sacr2(i,j,k,m) = 0.0_dp + enddo + enddo + enddo + enddo + + do m = 1, ntb_IN + do k = 1, ntb_t1 + do j = 1, ntb_r1 + do i = 1, ntb_r + tpi_qrfz(i,j,k,m) = 0.0_dp + tni_qrfz(i,j,k,m) = 0.0_dp + tpg_qrfz(i,j,k,m) = 0.0_dp + tnr_qrfz(i,j,k,m) = 0.0_dp + enddo + enddo + do j = 1, nbc + do i = 1, ntb_c + tpi_qcfz(i,j,k,m) = 0.0_dp + tni_qcfz(i,j,k,m) = 0.0_dp + enddo + enddo + enddo + enddo + + do j = 1, ntb_i1 + do i = 1, ntb_i + tps_iaus(i,j) = 0.0_dp + tni_iaus(i,j) = 0.0_dp + tpi_ide(i,j) = 0.0_dp + enddo + enddo + + do j = 1, nbc + do i = 1, nbr + t_Efrw(i,j) = 0.0 + enddo + do i = 1, nbs + t_Efsw(i,j) = 0.0 + enddo + enddo + + do k = 1, ntb_r + do j = 1, ntb_r1 + do i = 1, nbr + tnr_rev(i,j,k) = 0.0_dp + enddo + enddo + enddo + + do k = 1, nbc + do j = 1, ntb_c + do i = 1, nbc + tpc_wev(i,j,k) = 0.0_dp + tnc_wev(i,j,k) = 0.0_dp + enddo + enddo + enddo + + do m = 1, ntb_ark + do l = 1, ntb_arr + do k = 1, ntb_art + do j = 1, ntb_arw + do i = 1, ntb_arc + tnccn_act(i,j,k,l,m) = 1.0 + enddo + enddo + enddo + enddo + enddo + + if (mpirank==mpiroot) write (*,*)'creating microphysics lookup tables ... ' + if (mpirank==mpiroot) write (*,'(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') & + ' using: mu_c_o=',mu_c_o,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g + + !> - Call table_ccnact() to read a static file containing CCN activation of aerosols. The + !! data were created from a parcel model by Feingold & Heymsfield with + !! further changes by Eidhammer and Kriedenweis + if (mpirank==mpiroot) write(*,*) ' calling table_ccnAct routine' + call table_ccnAct(errmsg, errflg) + if (.not. errflg==0) return + + !> - Call table_efrw() and table_efsw() to creat collision efficiency table + !! between rain/snow and cloud water + if (mpirank==mpiroot) write(*,*) ' creating qc collision eff tables' + call table_Efrw + call table_Efsw + + !> - Call table_dropevap() to creat rain drop evaporation table + if (mpirank==mpiroot) write(*,*) ' creating rain evap table' + call table_dropEvap + + !> - Call qi_aut_qs() to create conversion of some ice mass into snow category + if (mpirank==mpiroot) write(*,*) ' creating ice converting to snow table' + call qi_aut_qs + + call cpu_time(etime) + if (mpirank==mpiroot) print '("Calculating Thompson tables part 1 took ",f10.3," seconds.")', etime-stime + + end if precomputed_tables_1 + + !> - Call radar_init() to initialize various constants for computing radar reflectivity + call cpu_time(stime) + xam_r = am_r + xbm_r = bm_r + xmu_r = mu_r + xam_s = am_s + xbm_s = bm_s + xmu_s = mu_s + xam_g = am_g(idx_bg1) + xbm_g = bm_g + xmu_g = mu_g + call radar_init + call cpu_time(etime) + if (mpirank==mpiroot) print '("Calling radar_init took ",f10.3," seconds.")', etime-stime + + if_not_iiwarm: if (.not. iiwarm) then + + precomputed_tables_2: if (.not.precomputed_tables) then + + call cpu_time(stime) + + !> - Call qr_acr_qg() to create rain collecting graupel & graupel collecting rain table + if (mpirank==mpiroot) write(*,*) ' creating rain collecting graupel table' + call cpu_time(stime) + call qr_acr_qg_par(dimNRHG) + call cpu_time(etime) + if (mpirank==mpiroot) print '("Computing rain collecting graupel table took ",f10.3," seconds.")', etime-stime + + !> - Call qr_acr_qs() to create rain collecting snow & snow collecting rain table + if (mpirank==mpiroot) write (*,*) ' creating rain collecting snow table' + call cpu_time(stime) + call qr_acr_qs_par + call cpu_time(etime) + if (mpirank==mpiroot) print '("Computing rain collecting snow table took ",f10.3," seconds.")', etime-stime + + !> - Call freezeh2o() to create cloud water and rain freezing (Bigg, 1953) table + if (mpirank==mpiroot) write(*,*) ' creating freezing of water drops table' + call cpu_time(stime) + call freezeH2O_par(threads) + call cpu_time(etime) + if (mpirank==mpiroot) print '("Computing freezing of water drops table took ",f10.3," seconds.")', etime-stime + + call cpu_time(etime) + if (mpirank==mpiroot) print '("Calculating Thompson tables part 2 took ",f10.3," seconds.")', etime-stime + + end if precomputed_tables_2 + + endif if_not_iiwarm + + if (mpirank==mpiroot) write(*,*) ' ... DONE microphysical lookup tables' + + endif if_micro_init + + end subroutine tempo_init + + !================================================================================================================= + ! This is a wrapper routine designed to transfer values from 3D to 1D. + ! Required microphysics variables are qv, qc, qr, nr, qi, ni, qs, qg + ! Optional microphysics variables are aerosol aware (nc, nwfa, nifa, nwfa2d, nifa2d), and hail aware (ng, qg) + + subroutine mp_gt_driver(qv, qc, qr, qi, qs, qg, qb, ni, nr, nc, ng, & + nwfa, nifa, nwfa2d, nifa2d, & + tt, th, pii, & + p, w, dz, dt_in, dt_inner, & + sedi_semi, decfl, lsm, & + RAINNC, RAINNCV, & + SNOWNC, SNOWNCV, & + ICENC, ICENCV, & + GRAUPELNC, GRAUPELNCV, SR, & + refl_10cm, diagflag, do_radar_ref, & + max_hail_diam_sfc, & + vt_dbz_wt, first_time_step, & + re_cloud, re_ice, re_snow, & + has_reqc, has_reqi, has_reqs, & + aero_ind_fdb, rand_perturb_on, & + kme_stoch, & + rand_pert, spp_prt_list, spp_var_list, & + spp_stddev_cutoff, n_var_spp, & + ids,ide, jds,jde, kds,kde, & ! domain dims + ims,ime, jms,jme, kms,kme, & ! memory dims + its,ite, jts,jte, kts,kte, & ! tile dims + fullradar_diag, istep, nsteps, & + errmsg, errflg, & + ! Extended diagnostics, array pointers + ! only associated if ext_diag flag is .true. + ext_diag, & + !vts1, txri, txrc, & + prw_vcdc, & + prw_vcde, tpri_inu, tpri_ide_d, & + tpri_ide_s, tprs_ide, tprs_sde_d, & + tprs_sde_s, tprg_gde_d, & + tprg_gde_s, tpri_iha, tpri_wfz, & + tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, & + tprg_rcs, tprs_rcs, & + tprr_rci, tprg_rcg, & + tprw_vcd_c, tprw_vcd_e, tprr_sml, & + tprr_gml, tprr_rcg, & + tprr_rcs, tprv_rev, tten3, qvten3, & + qrten3, qsten3, qgten3, qiten3, niten3, & + nrten3, ncten3, qcten3, & + pfils, pflls) + + !..Subroutine arguments + integer, intent(in):: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + real(wp), dimension(ims:ime, kms:kme, jms:jme), intent(inout):: & + qv, qc, qr, qi, qs, qg, ni, nr + real(wp), dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: & + tt, th + real(wp), dimension(ims:ime, kms:kme, jms:jme), optional, intent(in):: & + pii + real(wp), dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: & + nc, nwfa, nifa, qb, ng + real(wp), dimension(ims:ime, jms:jme), optional, intent(in):: nwfa2d, nifa2d + integer, dimension(ims:ime, jms:jme), intent(in):: lsm + real(wp), dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: & + re_cloud, re_ice, re_snow + real(wp), dimension(ims:ime, kms:kme, jms:jme), intent(inout):: pfils, pflls + integer, intent(in) :: rand_perturb_on, kme_stoch, n_var_spp + real(wp), dimension(:,:), optional, intent(in) :: rand_pert + real(wp), dimension(:), optional, intent(in) :: spp_prt_list + real(wp), dimension(:), intent(in) :: spp_stddev_cutoff + character(len=10), optional, dimension(:), intent(in) :: spp_var_list + integer, intent(in):: has_reqc, has_reqi, has_reqs + + real(wp), dimension(ims:ime, kms:kme, jms:jme), intent(in):: & + p, w, dz + real(wp), dimension(ims:ime, jms:jme), intent(inout):: & + RAINNC, RAINNCV, SR + real(wp), dimension(ims:ime, jms:jme), optional, intent(inout):: & + SNOWNC, SNOWNCV, & + ICENC, ICENCV, & + GRAUPELNC, GRAUPELNCV + real(wp), dimension(ims:ime, kms:kme, jms:jme), intent(inout):: & + refl_10cm + real(wp), dimension(ims:ime, jms:jme), intent(inout):: & + max_hail_diam_sfc + real(wp), dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: & + vt_dbz_wt + logical, intent(in) :: first_time_step + real(wp), intent(in):: dt_in, dt_inner + logical, intent(in) :: sedi_semi + integer, intent(in) :: decfl + ! To support subcycling: current step and maximum number of steps + integer, intent (in) :: istep, nsteps + logical, intent (in) :: fullradar_diag + ! Extended diagnostics, array pointers only associated if ext_diag flag is .true. + logical, intent (in) :: ext_diag + logical, optional, intent(in):: aero_ind_fdb + real(wp), optional, dimension(:,:,:), intent(inout):: & + !vts1, txri, txrc, & + prw_vcdc, & + prw_vcde, tpri_inu, tpri_ide_d, & + tpri_ide_s, tprs_ide, & + tprs_sde_d, tprs_sde_s, tprg_gde_d, & + tprg_gde_s, tpri_iha, tpri_wfz, & + tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, & + tprg_rcs, tprs_rcs, & + tprr_rci, tprg_rcg, & + tprw_vcd_c, tprw_vcd_e, tprr_sml, & + tprr_gml, tprr_rcg, & + tprr_rcs, tprv_rev, tten3, qvten3, & + qrten3, qsten3, qgten3, qiten3, niten3, & + nrten3, ncten3, qcten3 + + !..Local variables + real(wp), dimension(kts:kte):: & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & + ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, & + t1d, p1d, w1d, dz1d, rho, dBZ, pfil1, pfll1 + !..Extended diagnostics, single column arrays + real(wp), dimension(:), allocatable:: & + !vtsk1, txri1, txrc1, & + prw_vcdc1, & + prw_vcde1, tpri_inu1, tpri_ide1_d, & + tpri_ide1_s, tprs_ide1, & + tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, & + tprg_gde1_s, tpri_iha1, tpri_wfz1, & + tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,& + tprg_rcs1, tprs_rcs1, & + tprr_rci1, tprg_rcg1, & + tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, & + tprr_gml1, tprr_rcg1, & + tprr_rcs1, tprv_rev1, tten1, qvten1, & + qrten1, qsten1, qgten1, qiten1, niten1, & + nrten1, ncten1, qcten1 + + real(wp), dimension(kts:kte):: re_qc1d, re_qi1d, re_qs1d + + real(wp), dimension(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic + real(wp) :: dt, pptrain, pptsnow, pptgraul, pptice + real(wp) :: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max + real(wp) :: ygra1, zans1 + real(dp) :: lamg, lam_exp, lamr, N0_min, N0_exp + integer:: lsml + real(wp) :: rand1, rand2, rand3, rand_pert_max + integer:: i, j, k, m + integer:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr + integer:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr + integer:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr + integer:: i_start, j_start, i_end, j_end + logical, optional, intent(in) :: diagflag + integer, optional, intent(in) :: do_radar_ref + logical :: melti = .false. + integer :: ndt, it + + ! CCPP error handling + character(len=*), optional, intent( out) :: errmsg + integer, optional, intent( out) :: errflg + + ! CCPP + if (present(errmsg)) errmsg = '' + if (present(errflg)) errflg = 0 + + ! No need to test for every subcycling step + test_only_once: if (first_time_step .and. istep==1) then + ! Activate this code when removing the guard above + + if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & + (.not.present(tt) .and. .not.(present(th) .and. present(pii))) ) then + if (present(errmsg) .and. present(errflg)) then + write(errmsg, '(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' + errflg = 1 + return + else + write(*,'(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' + stop + end if + end if + + if (configs%aerosol_aware .and. (.not.present(nc) .or. & + .not.present(nwfa) .or. & + .not.present(nifa) .or. & + .not.present(nwfa2d) .or. & + .not.present(nifa2d) )) then + if (present(errmsg) .and. present(errflg)) then + write(errmsg, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', & + ' and nifa2d for aerosol-aware version of Thompson microphysics' + errflg = 1 + return + else + write(*, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', & + ' and nifa2d for aerosol-aware version of Thompson microphysics' + stop + end if + else if (merra2_aerosol_aware .and. (.not.present(nc) .or. & + .not.present(nwfa) .or. & + .not.present(nifa) )) then + if (present(errmsg) .and. present(errflg)) then + write(errmsg, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & + ' for merra2 aerosol-aware version of Thompson microphysics' + errflg = 1 + return + else + write(*, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & + ' for merra2 aerosol-aware version of Thompson microphysics' + stop + end if + else if (.not.configs%aerosol_aware .and. .not.merra2_aerosol_aware .and. & + (present(nwfa) .or. present(nifa) .or. present(nwfa2d) .or. present(nifa2d))) then + write(*,*) 'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware/merra2_aerosol_aware are FALSE' + end if + end if test_only_once + + ! These must be alwyas allocated + !allocate (vtsk1(kts:kte)) + !allocate (txri1(kts:kte)) + !allocate (txrc1(kts:kte)) + allocate_extended_diagnostics: if (ext_diag) then + allocate (prw_vcdc1(kts:kte)) + allocate (prw_vcde1(kts:kte)) + allocate (tpri_inu1(kts:kte)) + allocate (tpri_ide1_d(kts:kte)) + allocate (tpri_ide1_s(kts:kte)) + allocate (tprs_ide1(kts:kte)) + allocate (tprs_sde1_d(kts:kte)) + allocate (tprs_sde1_s(kts:kte)) + allocate (tprg_gde1_d(kts:kte)) + allocate (tprg_gde1_s(kts:kte)) + allocate (tpri_iha1(kts:kte)) + allocate (tpri_wfz1(kts:kte)) + allocate (tpri_rfz1(kts:kte)) + allocate (tprg_rfz1(kts:kte)) + allocate (tprs_scw1(kts:kte)) + allocate (tprg_scw1(kts:kte)) + allocate (tprg_rcs1(kts:kte)) + allocate (tprs_rcs1(kts:kte)) + allocate (tprr_rci1(kts:kte)) + allocate (tprg_rcg1(kts:kte)) + allocate (tprw_vcd1_c(kts:kte)) + allocate (tprw_vcd1_e(kts:kte)) + allocate (tprr_sml1(kts:kte)) + allocate (tprr_gml1(kts:kte)) + allocate (tprr_rcg1(kts:kte)) + allocate (tprr_rcs1(kts:kte)) + allocate (tprv_rev1(kts:kte)) + allocate (tten1(kts:kte)) + allocate (qvten1(kts:kte)) + allocate (qrten1(kts:kte)) + allocate (qsten1(kts:kte)) + allocate (qgten1(kts:kte)) + allocate (qiten1(kts:kte)) + allocate (niten1(kts:kte)) + allocate (nrten1(kts:kte)) + allocate (ncten1(kts:kte)) + allocate (qcten1(kts:kte)) + end if allocate_extended_diagnostics + + !+---+ + i_start = its + j_start = jts + i_end = ite + j_end = jte + + !..For idealized testing by developer. + ! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. & + ! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then + ! i_start = its + 2 + ! i_end = ite + ! j_start = jts + ! j_end = jte + ! endif + + ! dt = dt_in + RAINNC(:,:) = 0.0 + SNOWNC(:,:) = 0.0 + ICENC(:,:) = 0.0 + GRAUPELNC(:,:) = 0.0 + pcp_ra(:,:) = 0.0 + pcp_sn(:,:) = 0.0 + pcp_gr(:,:) = 0.0 + pcp_ic(:,:) = 0.0 + pfils(:,:,:) = 0.0 + pflls(:,:,:) = 0.0 + rand_pert_max = 0.0 + ndt = max(nint(dt_in/dt_inner),1) + dt = dt_in/ndt + if(dt_in .le. dt_inner) dt= dt_in + + !Get the Thompson MP SPP magnitude and standard deviation cutoff, + !then compute rand_pert_max + + if (rand_perturb_on .ne. 0) then + do k =1,n_var_spp + select case (spp_var_list(k)) + case('mp') + rand_pert_max = spp_prt_list(k)*spp_stddev_cutoff(k) + end select + enddo + endif + + do it = 1, ndt + + qc_max = 0. + qr_max = 0. + qs_max = 0. + qi_max = 0. + qg_max = 0 + ni_max = 0. + nr_max = 0. + imax_qc = 0 + imax_qr = 0 + imax_qi = 0 + imax_qs = 0 + imax_qg = 0 + imax_ni = 0 + imax_nr = 0 + jmax_qc = 0 + jmax_qr = 0 + jmax_qi = 0 + jmax_qs = 0 + jmax_qg = 0 + jmax_ni = 0 + jmax_nr = 0 + kmax_qc = 0 + kmax_qr = 0 + kmax_qi = 0 + kmax_qs = 0 + kmax_qg = 0 + kmax_ni = 0 + kmax_nr = 0 + + j_loop: do j = j_start, j_end + i_loop: do i = i_start, i_end + + !+---+-----------------------------------------------------------------+ + !..Introduce stochastic parameter perturbations by creating as many scalar rand1, rand2, ... + !.. variables as needed to perturb different pieces of microphysics. gthompsn 21Mar2018 + ! Setting spp_mp_opt to 1 gives graupel Y-intercept pertubations (2^0) + ! 2 gives cloud water distribution gamma shape parameter perturbations (2^1) + ! 4 gives CCN & IN activation perturbations (2^2) + ! 3 gives both 1+2 + ! 5 gives both 1+4 + ! 6 gives both 2+4 + ! 7 gives all 1+2+4 + ! For now (22Mar2018), standard deviation should be up to 0.75 and cut-off at 3.0 + ! stddev in order to constrain the various perturbations from being too extreme. + !+---+-----------------------------------------------------------------+ + rand1 = 0.0 + rand2 = 0.0 + rand3 = 0.0 + if (rand_perturb_on .ne. 0) then + if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1) + m = RSHIFT(ABS(rand_perturb_on),1) + if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1)*2. + m = RSHIFT(ABS(rand_perturb_on),2) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+rand_pert_max) + m = RSHIFT(ABS(rand_perturb_on),3) + endif + !+---+-----------------------------------------------------------------+ + + pptrain = 0. + pptsnow = 0. + pptgraul = 0. + pptice = 0. + RAINNCV(i,j) = 0. + IF ( PRESENT (snowncv) ) THEN + SNOWNCV(i,j) = 0. + ENDIF + IF ( PRESENT (icencv) ) THEN + ICENCV(i,j) = 0. + ENDIF + IF ( PRESENT (graupelncv) ) THEN + GRAUPELNCV(i,j) = 0. + ENDIF + SR(i,j) = 0. + + do k = kts, kte + if (present(tt)) then + t1d(k) = tt(i,k,j) + else + t1d(k) = th(i,k,j)*pii(i,k,j) + end if + p1d(k) = p(i,k,j) + w1d(k) = w(i,k,j) + dz1d(k) = dz(i,k,j) + qv1d(k) = qv(i,k,j) + qc1d(k) = qc(i,k,j) + qi1d(k) = qi(i,k,j) + qr1d(k) = qr(i,k,j) + qs1d(k) = qs(i,k,j) + qg1d(k) = qg(i,k,j) + ni1d(k) = ni(i,k,j) + nr1d(k) = nr(i,k,j) + rho(k) = RoverRv * p1d(k) / (R * t1d(k) * (qv1d(k)+RoverRv)) + + ! These arrays are always allocated and must be initialized + !vtsk1(k) = 0. + !txrc1(k) = 0. + !txri1(k) = 0. + initialize_extended_diagnostics: if (ext_diag) then + prw_vcdc1(k) = 0. + prw_vcde1(k) = 0. + tpri_inu1(k) = 0. + tpri_ide1_d(k) = 0. + tpri_ide1_s(k) = 0. + tprs_ide1(k) = 0. + tprs_sde1_d(k) = 0. + tprs_sde1_s(k) = 0. + tprg_gde1_d(k) = 0. + tprg_gde1_s(k) = 0. + tpri_iha1(k) = 0. + tpri_wfz1(k) = 0. + tpri_rfz1(k) = 0. + tprg_rfz1(k) = 0. + tprs_scw1(k) = 0. + tprg_scw1(k) = 0. + tprg_rcs1(k) = 0. + tprs_rcs1(k) = 0. + tprr_rci1(k) = 0. + tprg_rcg1(k) = 0. + tprw_vcd1_c(k) = 0. + tprw_vcd1_e(k) = 0. + tprr_sml1(k) = 0. + tprr_gml1(k) = 0. + tprr_rcg1(k) = 0. + tprr_rcs1(k) = 0. + tprv_rev1(k) = 0. + tten1(k) = 0. + qvten1(k) = 0. + qrten1(k) = 0. + qsten1(k) = 0. + qgten1(k) = 0. + qiten1(k) = 0. + niten1(k) = 0. + nrten1(k) = 0. + ncten1(k) = 0. + qcten1(k) = 0. + endif initialize_extended_diagnostics + enddo + lsml = lsm(i,j) + if (configs%aerosol_aware .or. merra2_aerosol_aware) then + do k = kts, kte + nc1d(k) = nc(i,k,j) + nwfa1d(k) = nwfa(i,k,j) + nifa1d(k) = nifa(i,k,j) + enddo + else + do k = kts, kte + if(lsml == 1) then + nc1d(k) = Nt_c_l / rho(k) + else + nc1d(k) = Nt_c_o / rho(k) + endif + nwfa1d(k) = nwfa_default + nifa1d(k) = nifa_default + enddo + endif + + ! ng and qb are optional hail-aware variables + if ((present(ng)) .and. (present(qb))) then + configs%hail_aware = .true. + do k = kts, kte + ng1d(k) = ng(i,k,j) + qb1d(k) = qb(i,k,j) + enddo + else + do k = kte, kts, -1 + ! This is the one-moment graupel formulation + if (qg1d(k) > R1) then + ygra1 = log10(max(1.e-9, qg1d(k)*rho(k))) + zans1 = 3.4 + 2.0/7.0*(ygra1+8.0) + ! zans1 = max(2.0, min(zans1, 6.0)) + N0_exp = max(gonv_min, min(10.0**(zans1), gonv_max)) + lam_exp = (n0_exp*am_g(idx_bg1)*cgg(1,1) / (rho(k)*qg1d(k)))**oge1 + lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg + ng1d(k) = cgg(2,1) * ogg3*rho(k) * qg1d(k) * lamg**bm_g / am_g(idx_bg1) + ng1d(k) = max(R2, (ng1d(k)/rho(k))) + qb1d(k) = qg1d(k) / rho_g(idx_bg1) + else + ng1d(k) = 0 + qb1d(k) = 0 + endif + enddo + endif + + !> - Call mp_thompson() + call mp_thompson_main(qv1d=qv1d, qc1d=qc1d, qi1d=qi1d, qr1d=qr1d, qs1d=qs1d, qg1d=qg1d, qb1d=qb1d, & + ni1d=ni1d, nr1d=nr1d, nc1d=nc1d, ng1d=ng1d, nwfa1d=nwfa1d, nifa1d=nifa1d, t1d=t1d, p1d=p1d, & + w1d=w1d, dzq=dz1d, pptrain=pptrain, pptsnow=pptsnow, pptgraul=pptgraul, pptice=pptice, & + rand1=rand1, rand2=rand3, rand3=rand3, & + ext_diag=ext_diag, sedi_semi=sedi_semi, decfl=decfl, & + prw_vcdc1=prw_vcdc1, & + prw_vcde1=prw_vcde1, & + tpri_inu1=tpri_inu1, tpri_ide1_d=tpri_ide1_d, tpri_ide1_s=tpri_ide1_s, tprs_ide1=tprs_ide1, & + tprs_sde1_d=tprs_sde1_d, tprs_sde1_s=tprs_sde1_s, & + tprg_gde1_d=tprg_gde1_d, tprg_gde1_s=tprg_gde1_s, tpri_iha1=tpri_iha1, tpri_wfz1=tpri_wfz1, & + tpri_rfz1=tpri_rfz1, tprg_rfz1=tprg_rfz1, tprs_scw1=tprs_scw1, tprg_scw1=tprg_scw1, & + tprg_rcs1=tprg_rcs1, tprs_rcs1=tprs_rcs1, tprr_rci1=tprr_rci1, & + tprg_rcg1=tprg_rcg1, tprw_vcd1_c=tprw_vcd1_c, & + tprw_vcd1_e=tprw_vcd1_e, tprr_sml1=tprr_sml1, tprr_gml1=tprr_gml1, tprr_rcg1=tprr_rcg1, & + tprr_rcs1=tprr_rcs1, tprv_rev1=tprv_rev1, & + tten1=tten1, qvten1=qvten1, qrten1=qrten1, qsten1=qsten1, & + qgten1=qgten1, qiten1=qiten1, niten1=niten1, nrten1=nrten1, ncten1=ncten1, qcten1=qcten1, & + pfil1=pfil1, pfll1=pfll1, lsml=lsml, & + kts=kts, kte=kte, dt=dt, ii=i, jj=j, configs=configs) + + + pcp_ra(i,j) = pcp_ra(i,j) + pptrain + pcp_sn(i,j) = pcp_sn(i,j) + pptsnow + pcp_gr(i,j) = pcp_gr(i,j) + pptgraul + pcp_ic(i,j) = pcp_ic(i,j) + pptice + RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice + RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice + IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN + ! Add ice to snow if separate ice not present + IF ( .NOT.PRESENT(icencv) .OR. .NOT.PRESENT(icenc) ) THEN + SNOWNCV(i,j) = pptsnow + pptice + SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice + ELSE + SNOWNCV(i,j) = pptsnow + SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + ENDIF + ENDIF + ! Use separate ice if present (as in FV3) + IF ( PRESENT(icencv) .AND. PRESENT(icenc) ) THEN + ICENCV(i,j) = pptice + ICENC(i,j) = ICENC(i,j) + pptice + ENDIF + IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN + GRAUPELNCV(i,j) = pptgraul + GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul + ENDIF + SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12) + + + + !..Reset lowest model level to initial state aerosols (fake sfc source). + !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol + !.. number tendency (number per kg per second). + if (configs%aerosol_aware) then + if ( present (aero_ind_fdb) ) then + if ( .not. aero_ind_fdb) then + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt + endif + else + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt + end if + + do k = kts, kte + nc(i,k,j) = nc1d(k) + nwfa(i,k,j) = nwfa1d(k) + nifa(i,k,j) = nifa1d(k) + enddo + endif + + if (merra2_aerosol_aware) then + do k = kts, kte + nc(i,k,j) = nc1d(k) + nwfa(i,k,j) = nwfa1d(k) + nifa(i,k,j) = nifa1d(k) + enddo + endif + + if ((present(ng)) .and. (present(qb))) then + do k = kts, kte + ng(i,k,j) = ng1d(k) + qb(i,k,j) = qb1d(k) + enddo + endif + + do k = kts, kte + qv(i,k,j) = qv1d(k) + qc(i,k,j) = qc1d(k) + qi(i,k,j) = qi1d(k) + qr(i,k,j) = qr1d(k) + qs(i,k,j) = qs1d(k) + qg(i,k,j) = qg1d(k) + ni(i,k,j) = ni1d(k) + nr(i,k,j) = nr1d(k) + pfils(i,k,j) = pfils(i,k,j) + pfil1(k) + pflls(i,k,j) = pflls(i,k,j) + pfll1(k) + if (present(tt)) then + tt(i,k,j) = t1d(k) + else + th(i,k,j) = t1d(k)/pii(i,k,j) + endif + + if (qc1d(k) .gt. qc_max) then + imax_qc = i + jmax_qc = j + kmax_qc = k + qc_max = qc1d(k) + elseif (qc1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qc ', qc1d(k), & + ' at i,j,k=', i,j,k + endif + if (qr1d(k) .gt. qr_max) then + imax_qr = i + jmax_qr = j + kmax_qr = k + qr_max = qr1d(k) + elseif (qr1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qr ', qr1d(k), & + ' at i,j,k=', i,j,k + endif + if (nr1d(k) .gt. nr_max) then + imax_nr = i + jmax_nr = j + kmax_nr = k + nr_max = nr1d(k) + elseif (nr1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative nr ', nr1d(k), & + ' at i,j,k=', i,j,k + endif + if (qs1d(k) .gt. qs_max) then + imax_qs = i + jmax_qs = j + kmax_qs = k + qs_max = qs1d(k) + elseif (qs1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qs ', qs1d(k), & + ' at i,j,k=', i,j,k + endif + if (qi1d(k) .gt. qi_max) then + imax_qi = i + jmax_qi = j + kmax_qi = k + qi_max = qi1d(k) + elseif (qi1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qi ', qi1d(k), & + ' at i,j,k=', i,j,k + endif + if (qg1d(k) .gt. qg_max) then + imax_qg = i + jmax_qg = j + kmax_qg = k + qg_max = qg1d(k) + elseif (qg1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qg ', qg1d(k), & + ' at i,j,k=', i,j,k + endif + if (ni1d(k) .gt. ni_max) then + imax_ni = i + jmax_ni = j + kmax_ni = k + ni_max = ni1d(k) + elseif (ni1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative ni ', ni1d(k), & + ' at i,j,k=', i,j,k + endif + if (qv1d(k) .lt. 0.0) then + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qv ', qv1d(k), & + ' at i,j,k=', i,j,k + if (k.lt.kte-2 .and. k.gt.kts+1) then + write(*,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) + qv(i,k,j) = max(1.e-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j))) + else + qv(i,k,j) = 1.e-7 + endif + endif + enddo + + assign_extended_diagnostics: if (ext_diag) then + do k=kts,kte + !vts1(i,k,j) = vtsk1(k) + !txri(i,k,j) = txri(i,k,j) + txri1(k) + !txrc(i,k,j) = txrc(i,k,j) + txrc1(k) + prw_vcdc(i,k,j) = prw_vcdc(i,k,j) + prw_vcdc1(k) + prw_vcde(i,k,j) = prw_vcde(i,k,j) + prw_vcde1(k) + tpri_inu(i,k,j) = tpri_inu(i,k,j) + tpri_inu1(k) + tpri_ide_d(i,k,j) = tpri_ide_d(i,k,j) + tpri_ide1_d(k) + tpri_ide_s(i,k,j) = tpri_ide_s(i,k,j) + tpri_ide1_s(k) + tprs_ide(i,k,j) = tprs_ide(i,k,j) + tprs_ide1(k) + tprs_sde_s(i,k,j) = tprs_sde_s(i,k,j) + tprs_sde1_s(k) + tprs_sde_d(i,k,j) = tprs_sde_d(i,k,j) + tprs_sde1_d(k) + tprg_gde_d(i,k,j) = tprg_gde_d(i,k,j) + tprg_gde1_d(k) + tprg_gde_s(i,k,j) = tprg_gde_s(i,k,j) + tprg_gde1_s(k) + tpri_iha(i,k,j) = tpri_iha(i,k,j) + tpri_iha1(k) + tpri_wfz(i,k,j) = tpri_wfz(i,k,j) + tpri_wfz1(k) + tpri_rfz(i,k,j) = tpri_rfz(i,k,j) + tpri_rfz1(k) + tprg_rfz(i,k,j) = tprg_rfz(i,k,j) + tprg_rfz1(k) + tprs_scw(i,k,j) = tprs_scw(i,k,j) + tprs_scw1(k) + tprg_scw(i,k,j) = tprg_scw(i,k,j) + tprg_scw1(k) + tprg_rcs(i,k,j) = tprg_rcs(i,k,j) + tprg_rcs1(k) + tprs_rcs(i,k,j) = tprs_rcs(i,k,j) + tprs_rcs1(k) + tprr_rci(i,k,j) = tprr_rci(i,k,j) + tprr_rci1(k) + tprg_rcg(i,k,j) = tprg_rcg(i,k,j) + tprg_rcg1(k) + tprw_vcd_c(i,k,j) = tprw_vcd_c(i,k,j) + tprw_vcd1_c(k) + tprw_vcd_e(i,k,j) = tprw_vcd_e(i,k,j) + tprw_vcd1_e(k) + tprr_sml(i,k,j) = tprr_sml(i,k,j) + tprr_sml1(k) + tprr_gml(i,k,j) = tprr_gml(i,k,j) + tprr_gml1(k) + tprr_rcg(i,k,j) = tprr_rcg(i,k,j) + tprr_rcg1(k) + tprr_rcs(i,k,j) = tprr_rcs(i,k,j) + tprr_rcs1(k) + tprv_rev(i,k,j) = tprv_rev(i,k,j) + tprv_rev1(k) + tten3(i,k,j) = tten3(i,k,j) + tten1(k) + qvten3(i,k,j) = qvten3(i,k,j) + qvten1(k) + qrten3(i,k,j) = qrten3(i,k,j) + qrten1(k) + qsten3(i,k,j) = qsten3(i,k,j) + qsten1(k) + qgten3(i,k,j) = qgten3(i,k,j) + qgten1(k) + qiten3(i,k,j) = qiten3(i,k,j) + qiten1(k) + niten3(i,k,j) = niten3(i,k,j) + niten1(k) + nrten3(i,k,j) = nrten3(i,k,j) + nrten1(k) + ncten3(i,k,j) = ncten3(i,k,j) + ncten1(k) + qcten3(i,k,j) = qcten3(i,k,j) + qcten1(k) + + enddo + endif assign_extended_diagnostics + + if (ndt>1 .and. it==ndt) then + + SR(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j))/(RAINNC(i,j)+1.e-12) + RAINNCV(i,j) = RAINNC(i,j) + IF ( PRESENT (snowncv) ) THEN + SNOWNCV(i,j) = SNOWNC(i,j) + ENDIF + IF ( PRESENT (icencv) ) THEN + ICENCV(i,j) = ICENC(i,j) + ENDIF + IF ( PRESENT (graupelncv) ) THEN + GRAUPELNCV(i,j) = GRAUPELNC(i,j) + ENDIF + endif + + ! Diagnostic calculations only for last step + ! if Thompson MP is called multiple times + last_step_only: IF ((ndt>1 .and. it==ndt) .or. & + (nsteps>1 .and. istep==nsteps) .or. & + (nsteps==1 .and. ndt==1)) THEN + +!! max_hail_diam_sfc(i,j) = hail_mass_99th_percentile(kts, kte, qg1d, t1d, p1d, qv1d) + + !> - Call calc_refl10cm() + + diagflag_present: IF ( PRESENT (diagflag) ) THEN + if (diagflag .and. do_radar_ref == 1) then + ! + ! Only set melti to true at the output times + if (fullradar_diag) then + melti=.true. + else + melti=.false. + endif + ! + if (present(vt_dbz_wt)) then + call calc_refl10cm (qv1d=qv1d, qc1d=qc1d, qr1d=qr1d, nr1d=nr1d, qs1d=qs1d, qg1d=qg1d, & + ng1d=ng1d, qb1d=qb1d, t1d=t1d, p1d=p1d, dBZ=dBZ, rand1=rand1, kts=kts, kte=kte, ii=i, jj=j, & + melti=melti, vt_dBZ=vt_dbz_wt(i,:,j), & + first_time_step=first_time_step, configs=configs) + else + call calc_refl10cm (qv1d=qv1d, qc1d=qc1d, qr1d=qr1d, nr1d=nr1d, qs1d=qs1d, qg1d=qg1d, & + ng1d=ng1d, qb1d=qb1d, t1d=t1d, p1d=p1d, dBZ=dBZ, rand1=rand1, kts=kts, kte=kte, ii=i, jj=j, & + melti=melti, configs=configs) + end if + do k = kts, kte + refl_10cm(i,k,j) = max(-35., dBZ(k)) + enddo + endif + ENDIF diagflag_present + + IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN + do k = kts, kte + re_qc1d(k) = re_qc_min + re_qi1d(k) = re_qi_min + re_qs1d(k) = re_qs_min + enddo + !> - Call calc_effectrad() + call calc_effectRad (t1d=t1d, p1d=p1d, qv1d=qv1d, qc1d=qc1d, & + nc1d=nc1d, qi1d=qi1d, ni1d=ni1d, qs1d=qs1d, & + re_qc1d=re_qc1d, re_qi1d=re_qi1d, re_qs1d=re_qs1d, & + kts=kts, kte=kte, lsml=lsml, configs=configs) + do k = kts, kte + re_cloud(i,k,j) = max(re_qc_min, min(re_qc1d(k), re_qc_max)) + re_ice(i,k,j) = max(re_qi_min, min(re_qi1d(k), re_qi_max)) + re_snow(i,k,j) = max(re_qs_min, min(re_qs1d(k), re_qs_max)) + enddo + ENDIF + ENDIF last_step_only + + enddo i_loop + enddo j_loop + + ! DEBUG - GT + ! write(*,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', & + ! 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', & + ! 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', & + ! 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', & + ! 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', & + ! 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', & + ! 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & + ! 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' + ! END DEBUG - GT + enddo ! end of nt loop + + do j = j_start, j_end + do k = kts, kte + do i = i_start, i_end + pfils(i,k,j) = pfils(i,k,j)/dt_in + pflls(i,k,j) = pflls(i,k,j)/dt_in + enddo + enddo + enddo + + ! These are always allocated + !deallocate (vtsk1) + !deallocate (txri1) + !deallocate (txrc1) + deallocate_extended_diagnostics: if (ext_diag) then + deallocate (prw_vcdc1) + deallocate (prw_vcde1) + deallocate (tpri_inu1) + deallocate (tpri_ide1_d) + deallocate (tpri_ide1_s) + deallocate (tprs_ide1) + deallocate (tprs_sde1_d) + deallocate (tprs_sde1_s) + deallocate (tprg_gde1_d) + deallocate (tprg_gde1_s) + deallocate (tpri_iha1) + deallocate (tpri_wfz1) + deallocate (tpri_rfz1) + deallocate (tprg_rfz1) + deallocate (tprs_scw1) + deallocate (tprg_scw1) + deallocate (tprg_rcs1) + deallocate (tprs_rcs1) + deallocate (tprr_rci1) + deallocate (tprg_rcg1) + deallocate (tprw_vcd1_c) + deallocate (tprw_vcd1_e) + deallocate (tprr_sml1) + deallocate (tprr_gml1) + deallocate (tprr_rcg1) + deallocate (tprr_rcs1) + deallocate (tprv_rev1) + deallocate (tten1) + deallocate (qvten1) + deallocate (qrten1) + deallocate (qsten1) + deallocate (qgten1) + deallocate (qiten1) + deallocate (niten1) + deallocate (nrten1) + deallocate (ncten1) + deallocate (qcten1) + end if deallocate_extended_diagnostics + + END SUBROUTINE mp_gt_driver + !> @} + + !>\ingroup aathompson + SUBROUTINE tempo_finalize() + + IMPLICIT NONE + + if (ALLOCATED(tcg_racg)) DEALLOCATE(tcg_racg) + if (ALLOCATED(tmr_racg)) DEALLOCATE(tmr_racg) + if (ALLOCATED(tcr_gacr)) DEALLOCATE(tcr_gacr) + if (ALLOCATED(tnr_racg)) DEALLOCATE(tnr_racg) + if (ALLOCATED(tnr_gacr)) DEALLOCATE(tnr_gacr) + + if (ALLOCATED(tcs_racs1)) DEALLOCATE(tcs_racs1) + if (ALLOCATED(tmr_racs1)) DEALLOCATE(tmr_racs1) + if (ALLOCATED(tcs_racs2)) DEALLOCATE(tcs_racs2) + if (ALLOCATED(tmr_racs2)) DEALLOCATE(tmr_racs2) + if (ALLOCATED(tcr_sacr1)) DEALLOCATE(tcr_sacr1) + if (ALLOCATED(tms_sacr1)) DEALLOCATE(tms_sacr1) + if (ALLOCATED(tcr_sacr2)) DEALLOCATE(tcr_sacr2) + if (ALLOCATED(tms_sacr2)) DEALLOCATE(tms_sacr2) + if (ALLOCATED(tnr_racs1)) DEALLOCATE(tnr_racs1) + if (ALLOCATED(tnr_racs2)) DEALLOCATE(tnr_racs2) + if (ALLOCATED(tnr_sacr1)) DEALLOCATE(tnr_sacr1) + if (ALLOCATED(tnr_sacr2)) DEALLOCATE(tnr_sacr2) + + if (ALLOCATED(tpi_qcfz)) DEALLOCATE(tpi_qcfz) + if (ALLOCATED(tni_qcfz)) DEALLOCATE(tni_qcfz) + + if (ALLOCATED(tpi_qrfz)) DEALLOCATE(tpi_qrfz) + if (ALLOCATED(tpg_qrfz)) DEALLOCATE(tpg_qrfz) + if (ALLOCATED(tni_qrfz)) DEALLOCATE(tni_qrfz) + if (ALLOCATED(tnr_qrfz)) DEALLOCATE(tnr_qrfz) + + if (ALLOCATED(tps_iaus)) DEALLOCATE(tps_iaus) + if (ALLOCATED(tni_iaus)) DEALLOCATE(tni_iaus) + if (ALLOCATED(tpi_ide)) DEALLOCATE(tpi_ide) + + if (ALLOCATED(t_Efrw)) DEALLOCATE(t_Efrw) + if (ALLOCATED(t_Efsw)) DEALLOCATE(t_Efsw) + + if (ALLOCATED(tnr_rev)) DEALLOCATE(tnr_rev) + if (ALLOCATED(tpc_wev)) DEALLOCATE(tpc_wev) + if (ALLOCATED(tnc_wev)) DEALLOCATE(tnc_wev) + + if (ALLOCATED(tnccn_act)) DEALLOCATE(tnccn_act) + + END SUBROUTINE tempo_finalize + +end module module_mp_tempo + !+---+-----------------------------------------------------------------+ + !ctrlL + !+---+-----------------------------------------------------------------+ + !+---+-----------------------------------------------------------------+ + diff --git a/physics/MP/TEMPO/mp_tempo.F90 b/physics/MP/TEMPO/mp_tempo.F90 new file mode 100644 index 000000000..24aecfe58 --- /dev/null +++ b/physics/MP/TEMPO/mp_tempo.F90 @@ -0,0 +1,1007 @@ +!>\file mp_tempo.F90 +!! This file contains aerosol-aware Tempo MP scheme. + + +!>\defgroup aatempo Aerosol-Aware Tempo MP Module +!! This module contains the aerosol-aware Tempo microphysics scheme. +module mp_tempo + + use mpi_f08 + use machine, only : kind_phys + + use module_mp_thompson_params + use module_mp_thompson_utils + use module_mp_thompson_main + use module_mp_tempo + + implicit none + + public :: mp_tempo_init, mp_tempo_run, mp_tempo_finalize + + private + + logical :: is_initialized = .False. + + integer, parameter :: ext_ndiag3d = 37 + + contains + +!> This subroutine is a wrapper around the actual tempo_init(). +!! \section arg_table_mp_tempo_init Argument Table +!! \htmlinclude mp_tempo_init.html +!! + subroutine mp_tempo_init(ncol, nlev, con_g, con_rd, con_eps, & + restart, imp_physics, & + imp_physics_thompson, convert_dry_rho, & + spechum, qc, qr, qi, qs, qg, ni, nr, & + chw, vh, & + is_aerosol_aware, merra2_aerosol_aware, & + is_hail_aware, & + nc, nwfa2d, nifa2d, & + nwfa, nifa, tgrs, prsl, phil, area, & + aerfld, mpicomm, mpirank, mpiroot, & + threads, diag3d, & + errmsg, errflg) + + implicit none + + ! Interface variables + integer, intent(in ) :: ncol + integer, intent(in ) :: nlev + real(kind_phys), intent(in ) :: con_g, con_rd, con_eps + logical, intent(in ) :: restart + integer, intent(in ) :: imp_physics + integer, intent(in ) :: imp_physics_thompson + ! Hydrometeors + logical, intent(in ) :: convert_dry_rho + real(kind_phys), intent(inout) :: spechum(:,:) + real(kind_phys), intent(inout) :: qc(:,:) + real(kind_phys), intent(inout) :: qr(:,:) + real(kind_phys), intent(inout) :: qi(:,:) + real(kind_phys), intent(inout) :: qs(:,:) + real(kind_phys), intent(inout) :: qg(:,:) + real(kind_phys), intent(inout) :: ni(:,:) + real(kind_phys), intent(inout) :: nr(:,:) + real(kind_phys), intent(inout), optional :: chw(:,:) + real(kind_phys), intent(inout), optional :: vh(:,:) + + ! Aerosols + logical, intent(in ) :: is_aerosol_aware + logical, intent(in ) :: merra2_aerosol_aware + logical, intent(in ) :: is_hail_aware + real(kind_phys), intent(inout), optional :: nc(:,:) + real(kind_phys), intent(inout), optional :: nwfa(:,:) + real(kind_phys), intent(inout), optional :: nifa(:,:) + real(kind_phys), intent(inout), optional :: nwfa2d(:) + real(kind_phys), intent(inout), optional :: nifa2d(:) + real(kind_phys), intent(in) :: aerfld(:,:,:) + ! State variables + real(kind_phys), intent(in ) :: tgrs(:,:) + real(kind_phys), intent(in ) :: prsl(:,:) + real(kind_phys), intent(in ) :: phil(:,:) + real(kind_phys), intent(in ) :: area(:) + ! MPI information + type(MPI_Comm), intent(in ) :: mpicomm + integer, intent(in ) :: mpirank + integer, intent(in ) :: mpiroot + ! Threading/blocking information + integer, intent(in ) :: threads + ! Extended diagnostics + real(kind_phys), intent(in ),optional :: diag3d(:,:,:) + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! + real(kind_phys) :: qv(1:ncol,1:nlev) ! kg kg-1 (water vapor mixing ratio) + real(kind_phys) :: hgt(1:ncol,1:nlev) ! m + real(kind_phys) :: rho(1:ncol,1:nlev) ! kg m-3 + real(kind_phys) :: orho(1:ncol,1:nlev) ! m3 kg-1 + real(kind_phys) :: nc_local(1:ncol,1:nlev) ! needed because nc is only allocated if is_aerosol_aware is true + ! + real (kind=kind_phys) :: h_01, z1, niIN3, niCCN3 + integer :: i, k + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + if (is_initialized) return + + ! Consistency checks + if (imp_physics/=imp_physics_thompson) then + write(errmsg,'(*(a))') "Logic error: namelist choice of microphysics is different from Thompson MP" + errflg = 1 + return + end if + + if (present(diag3d)) then + if (size(diag3d,dim=3) /= ext_ndiag3d) then + write(errmsg,'(*(a))') "Logic error: number of diagnostic 3d arrays from model does not match requirements" + errflg = 1 + return + end if + end if + + if (is_aerosol_aware .and. merra2_aerosol_aware) then + write(errmsg,'(*(a))') "Logic error: Only one Thompson aerosol option can be true, either is_aerosol_aware or merra2_aerosol_aware)" + errflg = 1 + return + end if + + ! Call Thompson init + call tempo_init(is_aerosol_aware_in=is_aerosol_aware, & + merra2_aerosol_aware_in=merra2_aerosol_aware, & + is_hail_aware_in=is_hail_aware, & + mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads, errmsg=errmsg, errflg=errflg) + if (errflg /= 0) return + + ! For restart runs, the init is done here + if (restart) then + is_initialized = .true. + return + end if + + ! Geopotential height in m2 s-2 to height in m + hgt = phil/con_g + + ! Ensure non-negative mass mixing ratios of all water variables + where(spechum<0) spechum = 1.0E-10 ! COMMENT, gthompsn, spechum should *never* be identically zero. + where(qc<0) qc = 0.0 + where(qr<0) qr = 0.0 + where(qi<0) qi = 0.0 + where(qs<0) qs = 0.0 + where(qg<0) qg = 0.0 + + !> - Convert specific humidity to water vapor mixing ratio. + !> - Also, hydrometeor variables are mass or number mixing ratio + !> - either kg of species per kg of dry air, or per kg of (dry + vapor). + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if + + + qv = spechum/(1.0_kind_phys-spechum) + + if (convert_dry_rho) then + qc = qc/(1.0_kind_phys-spechum) + qr = qr/(1.0_kind_phys-spechum) + qi = qi/(1.0_kind_phys-spechum) + qs = qs/(1.0_kind_phys-spechum) + qg = qg/(1.0_kind_phys-spechum) + + ni = ni/(1.0_kind_phys-spechum) + nr = nr/(1.0_kind_phys-spechum) + if (is_hail_aware) then + chw = chw/(1.0_kind_phys-spechum) + vh = vh/(1.0_kind_phys-spechum) + endif + if (is_aerosol_aware .or. merra2_aerosol_aware) then + nc = nc/(1.0_kind_phys-spechum) + nwfa = nwfa/(1.0_kind_phys-spechum) + nifa = nifa/(1.0_kind_phys-spechum) + end if + end if + + ! Density of moist air in kg m-3 and inverse density of air + rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) + orho = 1.0/rho + + ! Ensure we have 1st guess ice number where mass non-zero but no number. + where(qi .LE. 0.0) ni=0.0 + where(qi .GT. 0 .and. ni .LE. 0.0) ni = make_IceNumber(qi*rho, tgrs) * orho + where(qi .EQ. 0.0 .and. ni .GT. 0.0) ni=0.0 + + ! Ensure we have 1st guess rain number where mass non-zero but no number. + where(qr .LE. 0.0) nr=0.0 + where(qr .GT. 0 .and. nr .LE. 0.0) nr = make_RainNumber(qr*rho, tgrs) * orho + where(qr .EQ. 0.0 .and. nr .GT. 0.0) nr=0.0 + + if (is_hail_aware) then + where(qg .LE. 0.0) chw=0.0 + where(qg .LE. 0.0) vh=0.0 + endif + + !..Check for existing aerosol data, both CCN and IN aerosols. If missing + !.. fill in just a basic vertical profile, somewhat boundary-layer following. + if (is_aerosol_aware) then + + ! Potential cloud condensation nuclei (CCN) + if (MAXVAL(nwfa) .lt. eps) then + if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial CCN aerosols.' + do i = 1, ncol + if (hgt(i,1).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1).ge.2500.0) then + h_01 = 0.01 + else + h_01 = 0.8*cos(hgt(i,1)*0.001 - 1.0) + endif + niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 + nwfa(i,1) = naCCN1+naCCN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niCCN3) + z1 = hgt(i,2)-hgt(i,1) + nwfa2d(i) = nwfa(i,1) * 0.000196 * (50./z1) + do k = 2, nlev + nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) + enddo + enddo + else + if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosols are present.' + if (MAXVAL(nwfa2d) .lt. eps) then + !+---+-----------------------------------------------------------------+ + !..Scale the lowest level aerosol data into an emissions rate. This is + !.. very far from ideal, but need higher emissions where larger amount + !.. of (climo) existing and lesser emissions where there exists fewer to + !.. begin as a first-order simplistic approach. Later, proper connection to + !.. emission inventory would be better, but, for now, scale like this: + !.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per second per grid box unit + !.. that was tested as ~(20kmx20kmx50m = 2.E10 m**-3) + !+---+-----------------------------------------------------------------+ + if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' + do i = 1, ncol + z1 = hgt(i,2)-hgt(i,1) + nwfa2d(i) = nwfa(i,1) * 0.000196 * (50./z1) + enddo + else + if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosol surface emission rates are present.' + endif + endif + + ! Potential ice nuclei (IN) + if (MAXVAL(nifa) .lt. eps) then + if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial IN aerosols.' + do i = 1, ncol + if (hgt(i,1).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1).ge.2500.0) then + h_01 = 0.01 + else + h_01 = 0.8*cos(hgt(i,1)*0.001 - 1.0) + endif + niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01 + nifa(i,1) = naIN1+naIN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niIN3) + nifa2d(i) = 0. + do k = 2, nlev + nifa(i,k) = naIN1+naIN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niIN3) + enddo + enddo + else + if (mpirank==mpiroot) write(*,*) ' Apparently initial IN aerosols are present.' + if (MAXVAL(nifa2d) .lt. eps) then + if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial IN aerosol surface emission rates, set to zero.' + ! calculate IN surface flux here, right now just set to zero + nifa2d = 0. + else + if (mpirank==mpiroot) write(*,*) ' Apparently initial IN aerosol surface emission rates are present.' + endif + endif + + ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. + where(qc .LE. 0.0) nc=0.0 + where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho + where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc = 0.0 + + ! Ensure non-negative aerosol number concentrations. + where(nwfa .LE. 0.0) nwfa = 1.1E6 + where(nifa .LE. 0.0) nifa = naIN1*0.01 + + ! Copy to local array for calculating cloud effective radii below + nc_local = nc + + else if (merra2_aerosol_aware) then + + ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. + where(qc .LE. 0.0) nc=0.0 + where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho + where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc = 0.0 + + else + + ! Constant droplet concentration for single moment cloud water as in + ! module_mp_thompson.F90, only needed for effective radii calculation + nc_local = Nt_c_l/rho + + end if + + if (convert_dry_rho) then + !qc = qc/(1.0_kind_phys+qv) + !qr = qr/(1.0_kind_phys+qv) + !qi = qi/(1.0_kind_phys+qv) + !qs = qs/(1.0_kind_phys+qv) + !qg = qg/(1.0_kind_phys+qv) + + ni = ni/(1.0_kind_phys+qv) + nr = nr/(1.0_kind_phys+qv) + if (is_hail_aware) then + chw = chw/(1.0_kind_phys+qv) + vh = vh/(1.0_kind_phys+qv) + endif + if (is_aerosol_aware .or. merra2_aerosol_aware) then + nc = nc/(1.0_kind_phys+qv) + nwfa = nwfa/(1.0_kind_phys+qv) + nifa = nifa/(1.0_kind_phys+qv) + end if + end if + + is_initialized = .true. + + end subroutine mp_tempo_init + + +!> \section arg_table_mp_tempo_run Argument Table +!! \htmlinclude mp_tempo_run.html +!! +!>\ingroup aatempo +!>\section gen_tempo_hrrr Tempo MP General Algorithm +!>@{ + subroutine mp_tempo_run(ncol, nlev, con_g, con_rd, & + con_eps, convert_dry_rho, & + spechum, qc, qr, qi, qs, qg, ni, nr, & + chw, vh, & + is_aerosol_aware, is_hail_aware, & + merra2_aerosol_aware, nc, nwfa, nifa,& + nwfa2d, nifa2d, aero_ind_fdb, & + tgrs, prsl, phii, omega, & + sedi_semi, decfl, islmsk, dtp, & + dt_inner, & + first_time_step, istep, nsteps, & + prcp, rain, graupel, ice, snow, sr, & + refl_10cm, fullradar_diag, & + max_hail_diam_sfc, & + do_radar_ref, aerfld, & + mpicomm, mpirank, mpiroot, blkno, & + diag3d, reset_diag3d, & + spp_wts_mp, spp_mp, n_var_spp, & + spp_prt_list, spp_var_list, & + spp_stddev_cutoff, & + cplchm, pfi_lsan, pfl_lsan, & + errmsg, errflg) + + implicit none + + ! Interface variables + + ! Dimensions and constants + integer, intent(in ) :: ncol + integer, intent(in ) :: nlev + real(kind_phys), intent(in ) :: con_g + real(kind_phys), intent(in ) :: con_rd + real(kind_phys), intent(in ) :: con_eps + ! Hydrometeors + logical, intent(in ) :: convert_dry_rho + real(kind_phys), intent(inout) :: spechum(:,:) + real(kind_phys), intent(inout) :: qc(:,:) + real(kind_phys), intent(inout) :: qr(:,:) + real(kind_phys), intent(inout) :: qi(:,:) + real(kind_phys), intent(inout) :: qs(:,:) + real(kind_phys), intent(inout) :: qg(:,:) + real(kind_phys), intent(inout) :: ni(:,:) + real(kind_phys), intent(inout) :: nr(:,:) + real(kind_phys), optional, intent(inout) :: chw(:,:) + real(kind_phys), optional, intent(inout) :: vh(:,:) + ! Aerosols + logical, intent(in) :: is_aerosol_aware, fullradar_diag + logical, intent(in) :: merra2_aerosol_aware, is_hail_aware + real(kind_phys), optional, intent(inout) :: nc(:,:) + real(kind_phys), optional, intent(inout) :: nwfa(:,:) + real(kind_phys), optional, intent(inout) :: nifa(:,:) + real(kind_phys), optional, intent(in ) :: nwfa2d(:) + real(kind_phys), optional, intent(in ) :: nifa2d(:) + real(kind_phys), intent(in) :: aerfld(:,:,:) + logical, optional, intent(in ) :: aero_ind_fdb + ! State variables and timestep information + real(kind_phys), intent(inout) :: tgrs(:,:) + real(kind_phys), intent(in ) :: prsl(:,:) + real(kind_phys), intent(in ) :: phii(:,:) + real(kind_phys), intent(in ) :: omega(:,:) + integer, intent(in ) :: islmsk(:) + real(kind_phys), intent(in ) :: dtp + logical, intent(in ) :: first_time_step + integer, intent(in ) :: istep, nsteps + real, intent(in ) :: dt_inner + ! Precip/rain/snow/graupel fall amounts and fraction of frozen precip + real(kind_phys), intent(inout) :: prcp(:) + real(kind_phys), optional, intent(inout) :: rain(:) + real(kind_phys), optional, intent(inout) :: graupel(:) + real(kind_phys), optional, intent(inout) :: ice(:) + real(kind_phys), optional, intent(inout) :: snow(:) + real(kind_phys), intent( out) :: sr(:) + ! Radar reflectivity + real(kind_phys), intent(inout) :: refl_10cm(:,:) + real(kind_phys), intent(inout) :: max_hail_diam_sfc(:) + logical, intent(in ) :: do_radar_ref + logical, intent(in) :: sedi_semi + integer, intent(in) :: decfl + ! MPI and block information + integer, intent(in) :: blkno + type(MPI_Comm), intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + ! Extended diagnostic output + real(kind_phys), target, intent(inout), optional :: diag3d(:,:,:) + logical, intent(in) :: reset_diag3d + + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! SPP + integer, intent(in) :: spp_mp + integer, intent(in) :: n_var_spp + real(kind_phys), optional,intent(in) :: spp_wts_mp(:,:) + real(kind_phys), optional,intent(in) :: spp_prt_list(:) + character(len=10),optional,intent(in) :: spp_var_list(:) + real(kind_phys), optional,intent(in) :: spp_stddev_cutoff(:) + + logical, intent (in) :: cplchm + ! ice and liquid water 3d precipitation fluxes - only allocated if cplchm is .true. + real(kind=kind_phys), optional, intent(inout), dimension(:,:) :: pfi_lsan + real(kind=kind_phys), optional, intent(inout), dimension(:,:) :: pfl_lsan + + ! Local variables + + ! Reduced time step if subcycling is used + real(kind_phys) :: dtstep + integer :: ndt + ! Air density + real(kind_phys) :: rho(1:ncol,1:nlev) !< kg m-3 + ! Water vapor mixing ratio (instead of specific humidity) + real(kind_phys) :: qv(1:ncol,1:nlev) !< kg kg-1 + ! Vertical velocity and level width + real(kind_phys) :: w(1:ncol,1:nlev) !< m s-1 + real(kind_phys) :: dz(1:ncol,1:nlev) !< m + ! Rain/snow/graupel fall amounts + real(kind_phys) :: rain_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: graupel_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: ice_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: snow_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: delta_rain_mp(1:ncol) ! mm + real(kind_phys) :: delta_graupel_mp(1:ncol) ! mm + real(kind_phys) :: delta_ice_mp(1:ncol) ! mm + real(kind_phys) :: delta_snow_mp(1:ncol) ! mm + + real(kind_phys) :: pfils(1:ncol,1:nlev,1) + real(kind_phys) :: pflls(1:ncol,1:nlev,1) + ! Radar reflectivity + logical :: diagflag ! must be true if do_radar_ref is true, not used otherwise + integer :: do_radar_ref_mp ! integer instead of logical do_radar_ref + ! Effective cloud radii - turned off in CCPP (taken care off in radiation) + logical, parameter :: do_effective_radii = .false. + integer, parameter :: has_reqc = 0 + integer, parameter :: has_reqi = 0 + integer, parameter :: has_reqs = 0 + integer, parameter :: kme_stoch = 1 + integer :: spp_mp_opt + ! Dimensions used in mp_gt_driver + integer :: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + ! Pointer arrays for extended diagnostics + !real(kind_phys), dimension(:,:,:), pointer :: vts1 => null() + !real(kind_phys), dimension(:,:,:), pointer :: txri => null() + !real(kind_phys), dimension(:,:,:), pointer :: txrc => null() + real(kind_phys), dimension(:,:,:), pointer :: prw_vcdc => null() + real(kind_phys), dimension(:,:,:), pointer :: prw_vcde => null() + real(kind_phys), dimension(:,:,:), pointer :: tpri_inu => null() + real(kind_phys), dimension(:,:,:), pointer :: tpri_ide_d => null() + real(kind_phys), dimension(:,:,:), pointer :: tpri_ide_s => null() + real(kind_phys), dimension(:,:,:), pointer :: tprs_ide => null() + real(kind_phys), dimension(:,:,:), pointer :: tprs_sde_d => null() + real(kind_phys), dimension(:,:,:), pointer :: tprs_sde_s => null() + real(kind_phys), dimension(:,:,:), pointer :: tprg_gde_d => null() + real(kind_phys), dimension(:,:,:), pointer :: tprg_gde_s => null() + real(kind_phys), dimension(:,:,:), pointer :: tpri_iha => null() + real(kind_phys), dimension(:,:,:), pointer :: tpri_wfz => null() + real(kind_phys), dimension(:,:,:), pointer :: tpri_rfz => null() + real(kind_phys), dimension(:,:,:), pointer :: tprg_rfz => null() + real(kind_phys), dimension(:,:,:), pointer :: tprs_scw => null() + real(kind_phys), dimension(:,:,:), pointer :: tprg_scw => null() + real(kind_phys), dimension(:,:,:), pointer :: tprg_rcs => null() + real(kind_phys), dimension(:,:,:), pointer :: tprs_rcs => null() + real(kind_phys), dimension(:,:,:), pointer :: tprr_rci => null() + real(kind_phys), dimension(:,:,:), pointer :: tprg_rcg => null() + real(kind_phys), dimension(:,:,:), pointer :: tprw_vcd_c => null() + real(kind_phys), dimension(:,:,:), pointer :: tprw_vcd_e => null() + real(kind_phys), dimension(:,:,:), pointer :: tprr_sml => null() + real(kind_phys), dimension(:,:,:), pointer :: tprr_gml => null() + real(kind_phys), dimension(:,:,:), pointer :: tprr_rcg => null() + real(kind_phys), dimension(:,:,:), pointer :: tprr_rcs => null() + real(kind_phys), dimension(:,:,:), pointer :: tprv_rev => null() + real(kind_phys), dimension(:,:,:), pointer :: tten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: qvten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: qrten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: qsten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: qgten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: qiten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: niten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: nrten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: ncten3 => null() + real(kind_phys), dimension(:,:,:), pointer :: qcten3 => null() + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + if (first_time_step .and. istep==1 .and. blkno==1) then + write(*,*) 'Calling TEMPO microphysics run' + ! Check initialization state + if (.not.is_initialized) then + write(errmsg, fmt='((a))') 'mp_tempo_run called before mp_tempo_init' + errflg = 1 + return + end if + ! Check forr optional arguments of aerosol-aware microphysics + if (is_aerosol_aware .and. .not. (present(nc) .and. & + present(nwfa) .and. & + present(nifa) .and. & + present(nwfa2d) .and. & + present(nifa2d) )) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_tempo_run:', & + ' aerosol-aware microphysics require all of the', & + ' following optional arguments:', & + ' nc, nwfa, nifa, nwfa2d, nifa2d' + errflg = 1 + return + else if (merra2_aerosol_aware .and. .not. (present(nc) .and. & + present(nwfa) .and. & + present(nifa) )) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_tempo_run:', & + ' merra2 aerosol-aware microphysics require the', & + ' following optional arguments: nc, nwfa, nifa' + errflg = 1 + return + end if + ! Consistency cheecks - subcycling and inner loop at the same time are not supported + if (nsteps>1 .and. dt_inner < dtp) then + write(errmsg,'(*(a))') "Logic error: Subcycling and inner loop cannot be used at the same time" + errflg = 1 + return + else if (mpirank==mpiroot .and. nsteps>1) then + write(*,'(a,i0,a,a,f6.2,a)') 'Tempo MP is using ', nsteps, ' substep(s) per time step with an ', & + 'effective time step of ', dtp/real(nsteps, kind=kind_phys), ' seconds' + else if (mpirank==mpiroot .and. dt_inner < dtp) then + ndt = max(nint(dtp/dt_inner),1) + write(*,'(a,i0,a,a,f6.2,a)') 'Tempo MP is using ', ndt, ' inner loops per time step with an ', & + 'effective time step of ', dtp/real(ndt, kind=kind_phys), ' seconds' + end if + end if + + ! Set stochastic physics selection to apply all perturbations + if ( spp_mp==7 ) then + spp_mp_opt=7 + else + spp_mp_opt=0 + endif + + ! Set reduced time step if subcycling is used + if (nsteps>1) then + dtstep = dtp/real(nsteps, kind=kind_phys) + else + dtstep = dtp + end if + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if + + !> - Convert specific humidity to water vapor mixing ratio. + !> - Also, hydrometeor variables are mass or number mixing ratio + !> - either kg of species per kg of dry air, or per kg of (dry + vapor). + + ! DH* - do this only if istep == 1? Would be ok if it was + ! guaranteed that nothing else in the same subcycle group + ! was using these arrays, but it is somewhat dangerous. + qv = spechum/(1.0_kind_phys-spechum) + + if (convert_dry_rho) then + qc = qc/(1.0_kind_phys-spechum) + qr = qr/(1.0_kind_phys-spechum) + qi = qi/(1.0_kind_phys-spechum) + qs = qs/(1.0_kind_phys-spechum) + qg = qg/(1.0_kind_phys-spechum) + + ni = ni/(1.0_kind_phys-spechum) + nr = nr/(1.0_kind_phys-spechum) + if (is_hail_aware) then + chw = chw/(1.0_kind_phys-spechum) + vh = vh/(1.0_kind_phys-spechum) + endif + if (is_aerosol_aware .or. merra2_aerosol_aware) then + nc = nc/(1.0_kind_phys-spechum) + nwfa = nwfa/(1.0_kind_phys-spechum) + nifa = nifa/(1.0_kind_phys-spechum) + end if + end if + ! *DH + + !> - Density of air in kg m-3 + rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) + + !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 + w = -omega/(rho*con_g) + + !> - Layer width in m from geopotential in m2 s-2 + dz = (phii(:,2:nlev+1) - phii(:,1:nlev)) / con_g + + ! Accumulated values inside Tempo scheme, not used; + ! only use delta and add to inout variables (different units) + rain_mp = 0 + graupel_mp = 0 + ice_mp = 0 + snow_mp = 0 + delta_rain_mp = 0 + delta_graupel_mp = 0 + delta_ice_mp = 0 + delta_snow_mp = 0 + + ! Flags for calculating radar reflectivity; diagflag is redundant + if (do_radar_ref) then + diagflag = .true. + do_radar_ref_mp = 1 + else + diagflag = .false. + do_radar_ref_mp = 0 + end if + + ! Set internal dimensions + ids = 1 + ims = 1 + its = 1 + ide = ncol + ime = ncol + ite = ncol + jds = 1 + jms = 1 + jts = 1 + jde = 1 + jme = 1 + jte = 1 + kds = 1 + kms = 1 + kts = 1 + kde = nlev + kme = nlev + kte = nlev + if(cplchm) then + pfi_lsan = 0.0 + pfl_lsan = 0.0 + end if + + ! Set pointers for extended diagnostics + set_extended_diagnostic_pointers: if (present(diag3d)) then + if (reset_diag3d) then + diag3d = 0.0 + end if + !vts1 => diag3d(:,:,X:X) + !txri => diag3d(:,:,X:X) + !txrc => diag3d(:,:,X:X) + prw_vcdc => diag3d(:,:,1:1) + prw_vcde => diag3d(:,:,2:2) + tpri_inu => diag3d(:,:,3:3) + tpri_ide_d => diag3d(:,:,4:4) + tpri_ide_s => diag3d(:,:,5:5) + tprs_ide => diag3d(:,:,6:6) + tprs_sde_d => diag3d(:,:,7:7) + tprs_sde_s => diag3d(:,:,8:8) + tprg_gde_d => diag3d(:,:,9:9) + tprg_gde_s => diag3d(:,:,10:10) + tpri_iha => diag3d(:,:,11:11) + tpri_wfz => diag3d(:,:,12:12) + tpri_rfz => diag3d(:,:,13:13) + tprg_rfz => diag3d(:,:,14:14) + tprs_scw => diag3d(:,:,15:15) + tprg_scw => diag3d(:,:,16:16) + tprg_rcs => diag3d(:,:,17:17) + tprs_rcs => diag3d(:,:,18:18) + tprr_rci => diag3d(:,:,19:19) + tprg_rcg => diag3d(:,:,20:20) + tprw_vcd_c => diag3d(:,:,21:21) + tprw_vcd_e => diag3d(:,:,22:22) + tprr_sml => diag3d(:,:,23:23) + tprr_gml => diag3d(:,:,24:24) + tprr_rcg => diag3d(:,:,25:25) + tprr_rcs => diag3d(:,:,26:26) + tprv_rev => diag3d(:,:,27:27) + tten3 => diag3d(:,:,28:28) + qvten3 => diag3d(:,:,29:29) + qrten3 => diag3d(:,:,30:30) + qsten3 => diag3d(:,:,31:31) + qgten3 => diag3d(:,:,32:32) + qiten3 => diag3d(:,:,33:33) + niten3 => diag3d(:,:,34:34) + nrten3 => diag3d(:,:,35:35) + ncten3 => diag3d(:,:,36:36) + qcten3 => diag3d(:,:,37:37) + end if set_extended_diagnostic_pointers + !> - Call mp_gt_driver() with or without aerosols, with or without effective radii, ... + if (is_aerosol_aware .or. merra2_aerosol_aware) then + if (is_hail_aware) then + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, qb=vh, ni=ni, nr=nr, & + nc=nc, ng=chw, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & + sedi_semi=sedi_semi, decfl=decfl, lsm=islmsk, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + aero_ind_fdb=aero_ind_fdb, rand_perturb_on=spp_mp_opt, & + kme_stoch=kme_stoch, & + rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & + spp_prt_list=spp_prt_list, n_var_spp=n_var_spp, & + spp_stddev_cutoff=spp_stddev_cutoff, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + fullradar_diag=fullradar_diag, istep=istep, nsteps=nsteps, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & + ! Extended diagnostics + ext_diag=present(diag3d), & + ! vts1=vts1, txri=txri, txrc=txrc, & + prw_vcdc=prw_vcdc, & + prw_vcde=prw_vcde, tpri_inu=tpri_inu, tpri_ide_d=tpri_ide_d, & + tpri_ide_s=tpri_ide_s, tprs_ide=tprs_ide, & + tprs_sde_d=tprs_sde_d, & + tprs_sde_s=tprs_sde_s, tprg_gde_d=tprg_gde_d, & + tprg_gde_s=tprg_gde_s, tpri_iha=tpri_iha, & + tpri_wfz=tpri_wfz, tpri_rfz=tpri_rfz, tprg_rfz=tprg_rfz, & + tprs_scw=tprs_scw, tprg_scw=tprg_scw, tprg_rcs=tprg_rcs, & + tprs_rcs=tprs_rcs, & + tprr_rci=tprr_rci, tprg_rcg=tprg_rcg, tprw_vcd_c=tprw_vcd_c, & + tprw_vcd_e=tprw_vcd_e, tprr_sml=tprr_sml, tprr_gml=tprr_gml, & + tprr_rcg=tprr_rcg, tprr_rcs=tprr_rcs, & + tprv_rev=tprv_rev, tten3=tten3, & + qvten3=qvten3, qrten3=qrten3, qsten3=qsten3, qgten3=qgten3, & + qiten3=qiten3, niten3=niten3, nrten3=nrten3, ncten3=ncten3, & + qcten3=qcten3, pfils=pfils, pflls=pflls) + else + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & + nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & + sedi_semi=sedi_semi, decfl=decfl, lsm=islmsk, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + aero_ind_fdb=aero_ind_fdb, rand_perturb_on=spp_mp_opt, & + kme_stoch=kme_stoch, & + rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & + spp_prt_list=spp_prt_list, n_var_spp=n_var_spp, & + spp_stddev_cutoff=spp_stddev_cutoff, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + fullradar_diag=fullradar_diag, istep=istep, nsteps=nsteps, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & + ! Extended diagnostics + ext_diag=present(diag3d), & + ! vts1=vts1, txri=txri, txrc=txrc, & + prw_vcdc=prw_vcdc, & + prw_vcde=prw_vcde, tpri_inu=tpri_inu, tpri_ide_d=tpri_ide_d, & + tpri_ide_s=tpri_ide_s, tprs_ide=tprs_ide, & + tprs_sde_d=tprs_sde_d, & + tprs_sde_s=tprs_sde_s, tprg_gde_d=tprg_gde_d, & + tprg_gde_s=tprg_gde_s, tpri_iha=tpri_iha, & + tpri_wfz=tpri_wfz, tpri_rfz=tpri_rfz, tprg_rfz=tprg_rfz, & + tprs_scw=tprs_scw, tprg_scw=tprg_scw, tprg_rcs=tprg_rcs, & + tprs_rcs=tprs_rcs, & + tprr_rci=tprr_rci, tprg_rcg=tprg_rcg, tprw_vcd_c=tprw_vcd_c, & + tprw_vcd_e=tprw_vcd_e, tprr_sml=tprr_sml, tprr_gml=tprr_gml, & + tprr_rcg=tprr_rcg, tprr_rcs=tprr_rcs, & + tprv_rev=tprv_rev, tten3=tten3, & + qvten3=qvten3, qrten3=qrten3, qsten3=qsten3, qgten3=qgten3, & + qiten3=qiten3, niten3=niten3, nrten3=nrten3, ncten3=ncten3, & + qcten3=qcten3, pfils=pfils, pflls=pflls) + endif + else + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & + sedi_semi=sedi_semi, decfl=decfl, lsm=islmsk, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & + rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & + spp_prt_list=spp_prt_list, n_var_spp=n_var_spp, & + spp_stddev_cutoff=spp_stddev_cutoff, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + fullradar_diag=fullradar_diag, istep=istep, nsteps=nsteps, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & + ! Extended diagnostics + ext_diag=present(diag3d), & + ! vts1=vts1, txri=txri, txrc=txrc, & + prw_vcdc=prw_vcdc, & + prw_vcde=prw_vcde, tpri_inu=tpri_inu, tpri_ide_d=tpri_ide_d, & + tpri_ide_s=tpri_ide_s, tprs_ide=tprs_ide, & + tprs_sde_d=tprs_sde_d, & + tprs_sde_s=tprs_sde_s, tprg_gde_d=tprg_gde_d, & + tprg_gde_s=tprg_gde_s, tpri_iha=tpri_iha, & + tpri_wfz=tpri_wfz, tpri_rfz=tpri_rfz, tprg_rfz=tprg_rfz, & + tprs_scw=tprs_scw, tprg_scw=tprg_scw, tprg_rcs=tprg_rcs, & + tprs_rcs=tprs_rcs, & + tprr_rci=tprr_rci, tprg_rcg=tprg_rcg, tprw_vcd_c=tprw_vcd_c, & + tprw_vcd_e=tprw_vcd_e, tprr_sml=tprr_sml, tprr_gml=tprr_gml, & + tprr_rcg=tprr_rcg, tprr_rcs=tprr_rcs, & + tprv_rev=tprv_rev, tten3=tten3, & + qvten3=qvten3, qrten3=qrten3, qsten3=qsten3, qgten3=qgten3, & + qiten3=qiten3, niten3=niten3, nrten3=nrten3, ncten3=ncten3, & + qcten3=qcten3, pfils=pfils, pflls=pflls) + end if + if (errflg/=0) return + + ! DH* - do this only if istep == nsteps? Would be ok if it was + ! guaranteed that nothing else in the same subcycle group + ! was using these arrays, but it is somewhat dangerous. + + !> - Convert water vapor mixing ratio back to specific humidity + spechum = qv/(1.0_kind_phys+qv) + + if (convert_dry_rho) then + qc = qc/(1.0_kind_phys+qv) + qr = qr/(1.0_kind_phys+qv) + qi = qi/(1.0_kind_phys+qv) + qs = qs/(1.0_kind_phys+qv) + qg = qg/(1.0_kind_phys+qv) + + ni = ni/(1.0_kind_phys+qv) + nr = nr/(1.0_kind_phys+qv) + if (is_hail_aware) then + chw = chw/(1.0_kind_phys+qv) + vh = vh/(1.0_kind_phys+qv) + endif + if (is_aerosol_aware .or. merra2_aerosol_aware) then + nc = nc/(1.0_kind_phys+qv) + nwfa = nwfa/(1.0_kind_phys+qv) + nifa = nifa/(1.0_kind_phys+qv) + end if + end if + ! *DH + + !> - Convert rainfall deltas from mm to m (on physics timestep); add to inout variables + ! "rain" in Tempo MP refers to precipitation (total of liquid rainfall+snow+graupel+ice) + prcp = prcp + max(0.0, delta_rain_mp/1000.0_kind_phys) + graupel = graupel + max(0.0, delta_graupel_mp/1000.0_kind_phys) + ice = ice + max(0.0, delta_ice_mp/1000.0_kind_phys) + snow = snow + max(0.0, delta_snow_mp/1000.0_kind_phys) + rain = rain + max(0.0, (delta_rain_mp - (delta_graupel_mp + delta_ice_mp + delta_snow_mp))/1000.0_kind_phys) + + ! Recompute sr at last subcycling step + if (nsteps>1 .and. istep == nsteps) then + ! Unlike inside mp_gt_driver, rain does not contain frozen precip + sr = (snow + graupel + ice)/(rain + snow + graupel + ice +1.e-12) + end if + + ! output instantaneous ice/snow and rain water 3d precipitation fluxes + if(cplchm) then + pfi_lsan(:,:) = pfils(:,:,1) + pfl_lsan(:,:) = pflls(:,:,1) + end if + + unset_extended_diagnostic_pointers: if (present(diag3d)) then + !vts1 => null() + !txri => null() + !txrc => null() + prw_vcdc => null() + prw_vcde => null() + tpri_inu => null() + tpri_ide_d => null() + tpri_ide_s => null() + tprs_ide => null() + tprs_sde_d => null() + tprs_sde_s => null() + tprg_gde_d => null() + tprg_gde_s => null() + tpri_iha => null() + tpri_wfz => null() + tpri_rfz => null() + tprg_rfz => null() + tprs_scw => null() + tprg_scw => null() + tprg_rcs => null() + tprs_rcs => null() + tprr_rci => null() + tprg_rcg => null() + tprw_vcd_c => null() + tprw_vcd_e => null() + tprr_sml => null() + tprr_gml => null() + tprr_rcg => null() + tprr_rcs => null() + tprv_rev => null() + tten3 => null() + qvten3 => null() + qrten3 => null() + qsten3 => null() + qgten3 => null() + qiten3 => null() + niten3 => null() + nrten3 => null() + ncten3 => null() + qcten3 => null() + end if unset_extended_diagnostic_pointers + + end subroutine mp_tempo_run +!>@} + +!> \section arg_table_mp_tempo_finalize Argument Table +!! \htmlinclude mp_tempo_finalize.html +!! + subroutine mp_tempo_finalize(errmsg, errflg) + + implicit none + + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + if (.not.is_initialized) return + + call tempo_finalize() + + is_initialized = .false. + + end subroutine mp_tempo_finalize + + subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + ! To calculate nifa and nwfa from bins of aerosols. + ! In GOCART and MERRA2, aerosols are given as mixing ratio (kg/kg). To + ! convert from kg/kg to #/kg, the "unit mass" (mass of one particle) + ! within the mass bins is calculated. A lognormal size distribution + ! within aerosol bins is used to find the size based upon the median + ! mass. NIFA is mainly summarized over five dust bins and NWFA over the + ! other 10 bins. The parameters besides each bins are carefully tuned + ! for a good performance of the scheme. + ! + ! The fields for the last index of the aerfld array + ! are specified as below. + ! 1: dust bin 1, 0.1 to 1.0 micrometers + ! 2: dust bin 2, 1.0 to 1.8 micrometers + ! 3: dust bin 3, 1.8 to 3.0 micrometers + ! 4: dust bin 4, 3.0 to 6.0 micrometers + ! 5: dust bin 5, 6.0 to 10.0 micrometers + ! 6: sea salt bin 1, 0.03 to 0.1 micrometers + ! 7: sea salt bin 2, 0.1 to 0.5 micrometers + ! 8: sea salt bin 3, 0.5 to 1.5 micrometers + ! 9: sea salt bin 4, 1.5 to 5.0 micrometers + ! 10: sea salt bin 5, 5.0 to 10.0 micrometers + ! 11: Sulfate, 0.35 (mean) micrometers + ! 15: water-friendly organic carbon, 0.35 (mean) micrometers + ! + ! Bin densities are as follows: + ! 1: dust bin 1: 2500 kg/m2 + ! 2-5: dust bin 2-5: 2650 kg/m2 + ! 6-10: sea salt bins 6-10: 2200 kg/m2 + ! 11: sulfate: 1700 kg/m2 + ! 15: organic carbon: 1800 kg/m2 + + implicit none + integer, intent(in)::ncol, nlev + real (kind=kind_phys), dimension(:,:,:), intent(in) :: aerfld + real (kind=kind_phys), dimension(:,:), intent(out ):: nifa, nwfa + + nifa=(aerfld(:,:,1)/4.0737762+aerfld(:,:,2)/30.459203+aerfld(:,:,3)/153.45048+ & + aerfld(:,:,4)/1011.5142+ aerfld(:,:,5)/5683.3501)*1.e15 + + nwfa=((aerfld(:,:,6)/0.0045435214+aerfld(:,:,7)/0.2907854+aerfld(:,:,8)/12.91224+ & + aerfld(:,:,9)/206.2216+ aerfld(:,:,10)/4326.23)*9.+aerfld(:,:,11)/0.3053104*5+ & + aerfld(:,:,15)/0.3232698*8)*1.e15 + end subroutine get_niwfa + +end module mp_tempo diff --git a/physics/MP/TEMPO/mp_tempo.meta b/physics/MP/TEMPO/mp_tempo.meta new file mode 100644 index 000000000..c24dceec2 --- /dev/null +++ b/physics/MP/TEMPO/mp_tempo.meta @@ -0,0 +1,858 @@ +[ccpp-table-properties] + name = mp_tempo + type = scheme + dependencies = ../../hooks/machine.F + dependencies = ../module_mp_radar.F90 + dependencies = tempo/module_mp_thompson_params.F90 + dependencies = tempo/module_mp_thompson_utils.F90 + dependencies = tempo/module_mp_thompson_main.F90 + dependencies = module_mp_tempo.F90 + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_init + type = scheme +[ncol] + standard_name = horizontal_dimension + long_name = horizontal dimension + units = count + dimensions = () + type = integer + intent = in +[nlev] + standard_name = vertical_layer_dimension + long_name = number of vertical levels + units = count + dimensions = () + type = integer + intent = in +[con_g] + standard_name = gravitational_acceleration + long_name = gravitational acceleration + units = m s-2 + dimensions = () + type = real + kind = kind_phys + intent = in +[con_rd] + standard_name = gas_constant_of_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in +[imp_physics] + standard_name = control_for_microphysics_scheme + long_name = choice of microphysics scheme + units = flag + dimensions = () + type = integer + intent = in +[imp_physics_thompson] + standard_name = identifier_for_thompson_microphysics_scheme + long_name = choice of Thompson microphysics scheme + units = flag + dimensions = () + type = integer + intent = in +[convert_dry_rho] + standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air + long_name = flag for converting hydrometeors from moist to dry air + units = flag + dimensions = () + type = logical + intent = in +[spechum] + standard_name = specific_humidity + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qc] + standard_name = cloud_liquid_water_mixing_ratio + long_name = cloud water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qr] + standard_name = rain_mixing_ratio + long_name = rain water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qi] + standard_name = cloud_ice_mixing_ratio + long_name = ice water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qs] + standard_name = snow_mixing_ratio + long_name = snow water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qg] + standard_name = graupel_mixing_ratio + long_name = graupel mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[chw] + standard_name = mass_number_concentration_of_graupel_in_air + long_name = graupel number concentration + units = kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[vh] + standard_name = graupel_volume + long_name = graupel particle volume + units = m3 kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[ni] + standard_name = mass_number_concentration_of_cloud_ice_water_crystals_in_air + long_name = ice number concentration + units = kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[nr] + standard_name = mass_number_concentration_of_rain_water_in_air + long_name = rain number concentration + units = kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[is_aerosol_aware] + standard_name = flag_for_aerosol_physics + long_name = flag for aerosol-aware physics + units = flag + dimensions = () + type = logical + intent = in +[is_hail_aware] + standard_name = flag_for_hail_physics + long_name = flag for hail-aware physics + units = flag + dimensions = () + type = logical + intent = in +[merra2_aerosol_aware] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in +[nc] + standard_name = mass_number_concentration_of_cloud_liquid_water_particles_in_air + long_name = cloud droplet number concentration + units = kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nwfa2d] + standard_name = tendency_of_hygroscopic_aerosols_at_surface_adjacent_layer + long_name = instantaneous fake water-friendly surface aerosol source + units = kg-1 s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nifa2d] + standard_name = tendency_of_nonhygroscopic_ice_nucleating_aerosols_at_surface_adjacent_layer + long_name = instantaneous fake ice-friendly surface aerosol source + units = kg-1 s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nwfa] + standard_name = mass_number_concentration_of_hygroscopic_aerosols + long_name = number concentration of water-friendly aerosols + units = kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nifa] + standard_name = mass_number_concentration_of_nonhygroscopic_ice_nucleating_aerosols + long_name = number concentration of ice-friendly aerosols + units = kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[tgrs] + standard_name = air_temperature + long_name = model layer mean temperature + units = K + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[prsl] + standard_name = air_pressure + long_name = mean layer pressure + units = Pa + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[phil] + standard_name = geopotential + long_name = geopotential at model layer centers + units = m2 s-2 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[area] + standard_name = cell_area + long_name = area of the grid cell + units = m2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in +[aerfld] + standard_name = mass_mixing_ratio_of_aerosol_from_gocart_or_merra2 + long_name = mass mixing ratio of aerosol from gocart or merra2 + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension,number_of_aerosol_tracers_MG) + type = real + kind = kind_phys + intent = in +[mpicomm] + standard_name = mpi_communicator + long_name = MPI communicator + units = index + dimensions = () + type = MPI_Comm + intent = in +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in +[threads] + standard_name = number_of_openmp_threads + long_name = number of OpenMP threads available to scheme + units = count + dimensions = () + type = integer + intent = in +[diag3d] + standard_name = extended_diagnostics_output_from_thompson_microphysics + long_name = set of 3d arrays for extended diagnostics output from thompson microphysics + units = none + dimensions = (horizontal_dimension,vertical_layer_dimension,number_of_3d_diagnostic_output_arrays_from_thompson_microphysics) + type = real + kind = kind_phys + intent = in + optional = True +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_run + type = scheme +[ncol] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[nlev] + standard_name = vertical_layer_dimension + long_name = number of vertical levels + units = count + dimensions = () + type = integer + intent = in +[con_g] + standard_name = gravitational_acceleration + long_name = gravitational acceleration + units = m s-2 + dimensions = () + type = real + kind = kind_phys + intent = in +[con_rd] + standard_name = gas_constant_of_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[convert_dry_rho] + standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air + long_name = flag for converting hydrometeors from moist to dry air + units = flag + dimensions = () + type = logical + intent = in +[spechum] + standard_name = specific_humidity_of_new_state + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qc] + standard_name = cloud_liquid_water_mixing_ratio_of_new_state + long_name = cloud water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qr] + standard_name = rain_mixing_ratio_of_new_state + long_name = rain water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qi] + standard_name = cloud_ice_mixing_ratio_of_new_state + long_name = ice water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qs] + standard_name = snow_mixing_ratio_of_new_state + long_name = snow water mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[qg] + standard_name = graupel_mixing_ratio_of_new_state + long_name = graupel mixing ratio wrt dry+vapor (no condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[chw] + standard_name = mass_number_concentration_of_graupel_of_new_state + long_name = graupel number concentration + units = kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[vh] + standard_name = graupel_volume_of_new_state + long_name = graupel particle volume + units = m3 kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[ni] + standard_name = mass_number_concentration_of_cloud_ice_water_crystals_in_air_of_new_state + long_name = ice number concentration + units = kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[nr] + standard_name = mass_number_concentration_of_rain_of_new_state + long_name = rain number concentration + units = kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[is_aerosol_aware] + standard_name = flag_for_aerosol_physics + long_name = flag for aerosol-aware physics + units = flag + dimensions = () + type = logical + intent = in +[is_hail_aware] + standard_name = flag_for_hail_physics + long_name = flag for hail-aware physics + units = flag + dimensions = () + type = logical + intent = in +[merra2_aerosol_aware] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in +[nc] + standard_name = mass_number_concentration_of_cloud_liquid_water_particles_in_air_of_new_state + long_name = cloud droplet number concentration + units = kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nwfa] + standard_name = mass_number_concentration_of_hygroscopic_aerosols_of_new_state + long_name = number concentration of water-friendly aerosols + units = kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nifa] + standard_name = mass_number_concentration_of_nonhygroscopic_ice_nucleating_aerosols_of_new_state + long_name = number concentration of ice-friendly aerosols + units = kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[nwfa2d] + standard_name = tendency_of_hygroscopic_aerosols_at_surface_adjacent_layer + long_name = instantaneous fake water-friendly surface aerosol source + units = kg-1 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[nifa2d] + standard_name = tendency_of_nonhygroscopic_ice_nucleating_aerosols_at_surface_adjacent_layer + long_name = instantaneous fake ice-friendly surface aerosol source + units = kg-1 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[aero_ind_fdb] + standard_name = do_smoke_aerosol_indirect_feedback + long_name = flag for wfa ifa emission indirect feedback + units = flag + dimensions = () + type = logical + intent = in +[tgrs] + standard_name = air_temperature_of_new_state + long_name = model layer mean temperature + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[prsl] + standard_name = air_pressure + long_name = mean layer pressure + units = Pa + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[phii] + standard_name = geopotential_at_interface + long_name = geopotential at model layer interfaces + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_interface_dimension) + type = real + kind = kind_phys + intent = in +[omega] + standard_name = lagrangian_tendency_of_air_pressure + long_name = layer mean vertical velocity + units = Pa s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[sedi_semi] + standard_name = flag_for_semi_Lagrangian_sedi_rain + long_name = flag for semi Lagrangian sedi of rain + units = flag + dimensions = () + type = logical + intent = in +[decfl] + standard_name = deformed_CFL_factor + long_name = deformed CFL factor + units = count + dimensions = () + type = integer + intent = in +[islmsk] + standard_name = sea_land_ice_mask + long_name = sea/land/ice mask (=0/1/2) + units = flag + dimensions = (horizontal_loop_extent) + type = integer + intent = in +[dtp] + standard_name = timestep_for_physics + long_name = physics timestep + units = s + dimensions = () + type = real + kind = kind_phys + intent = in +[dt_inner] + standard_name = time_step_for_inner_loop + long_name = time step for inner loop + units = s + dimensions = () + type = real + kind = kind_phys + intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[istep] + standard_name = ccpp_loop_counter + long_name = loop counter for subcycling loops in CCPP + units = index + dimensions = () + type = integer + intent = in +[nsteps] + standard_name = ccpp_loop_extent + long_name = loop extent for subcycling loops in CCPP + units = count + dimensions = () + type = integer + intent = in +[prcp] + standard_name = lwe_thickness_of_explicit_precipitation_amount + long_name = explicit precipitation (rain, ice, snow, graupel) on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[rain] + standard_name = lwe_thickness_of_explicit_rain_amount + long_name = explicit rain fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[graupel] + standard_name = lwe_thickness_of_graupel_amount + long_name = graupel fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[ice] + standard_name = lwe_thickness_of_ice_amount + long_name = ice fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[snow] + standard_name = lwe_thickness_of_snow_amount + long_name = snow fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[sr] + standard_name = ratio_of_snowfall_to_rainfall + long_name = ratio of snowfall to large-scale rainfall + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[refl_10cm] + standard_name = radar_reflectivity_10cm + long_name = instantaneous refl_10cm + units = dBZ + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[max_hail_diam_sfc] + standard_name = max_hail_diameter_sfc + long_name = instantaneous maximum hail diameter at lowest model level + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[fullradar_diag] + standard_name = do_full_radar_reflectivity + long_name = flag for computing full radar reflectivity + units = flag + dimensions = () + type = logical + intent = in +[do_radar_ref] + standard_name = flag_for_radar_reflectivity + long_name = flag for radar reflectivity + units = flag + dimensions = () + type = logical + intent = in +[aerfld] + standard_name = mass_mixing_ratio_of_aerosol_from_gocart_or_merra2 + long_name = mass mixing ratio of aerosol from gocart or merra2 + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_aerosol_tracers_MG) + type = real + kind = kind_phys + intent = in +[mpicomm] + standard_name = mpi_communicator + long_name = MPI communicator + units = index + dimensions = () + type = MPI_Comm + intent = in +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in +[blkno] + standard_name = ccpp_block_number + long_name = number of block for explicit data blocking in CCPP + units = index + dimensions = () + type = integer + intent = in +[diag3d] + standard_name = extended_diagnostics_output_from_thompson_microphysics + long_name = set of 3d arrays for extended diagnostics output from thompson microphysics + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_3d_diagnostic_output_arrays_from_thompson_microphysics) + type = real + kind = kind_phys + intent = inout + optional = True +[reset_diag3d] + standard_name = flag_reset_extended_diagnostics_output_arrays_from_thompson_microphysics + long_name = flag for resetting extended diagnostics output arrays from thompson microphysics + units = flag + dimensions = () + type = logical + intent = in +[spp_wts_mp] + standard_name = spp_weights_for_microphysics_scheme + long_name = spp weights for microphysics scheme + units = 1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + intent = in + optional = True +[spp_mp] + standard_name = control_for_microphysics_spp_perturbations + long_name = control for microphysics spp perturbations + units = count + dimensions = () + type = integer + intent = in +[n_var_spp] + standard_name = number_of_perturbed_spp_schemes + long_name = number of perturbed spp schemes + units = count + dimensions = () + type = integer + intent = in +[spp_prt_list] + standard_name = magnitude_of_spp_perturbations + long_name = magnitude of spp perturbations + units = 1 + dimensions = (number_of_perturbed_spp_schemes) + type = real + kind = kind_phys + intent = in + optional = True +[spp_stddev_cutoff] + standard_name = magnitude_of_spp_standard_deviation_cutoff + long_name = magnitude of spp standard deviation cutoff + units = 1 + dimensions = (number_of_perturbed_spp_schemes) + type = real + kind = kind_phys + intent = in + optional = True +[spp_var_list] + standard_name = perturbed_spp_schemes + long_name = perturbed spp schemes + units = none + dimensions = (number_of_perturbed_spp_schemes) + type = character + kind = len=10 + intent = in + optional = True +[cplchm] + standard_name = flag_for_chemistry_coupling + long_name = flag controlling cplchm collection (default off) + units = flag + dimensions = () + type = logical + intent = in +[pfi_lsan] + standard_name = ice_flux_due_to_large_scale_precipitation + long_name = instantaneous 3D flux of ice from nonconvective precipitation + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[pfl_lsan] + standard_name = liquid_flux_due_to_large_scale_precipitation + long_name = instantaneous 3D flux of liquid water from nonconvective precipitation + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_finalize + type = scheme +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/MP/TEMPO/mp_tempo_post.F90 b/physics/MP/TEMPO/mp_tempo_post.F90 new file mode 100644 index 000000000..e8b531f83 --- /dev/null +++ b/physics/MP/TEMPO/mp_tempo_post.F90 @@ -0,0 +1,150 @@ +module mp_tempo_post + + use mpi_f08 + use machine, only : kind_phys + + implicit none + + public :: mp_tempo_post_init, mp_tempo_post_run, mp_tempo_post_finalize + + private + + logical :: is_initialized = .false. + + logical :: apply_limiter + +contains + +!! \section arg_table_mp_tempo_post_init Argument Table +!! \htmlinclude mp_tempo_post_init.html +!! + subroutine mp_tempo_post_init(ttendlim, errmsg, errflg) + + implicit none + + ! Interface variables + real(kind_phys), intent(in) :: ttendlim + + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Local variables + integer :: i + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Check initialization state + if (is_initialized) return + + if (ttendlim < 0) then + apply_limiter = .false. + else + apply_limiter = .true. + end if + + is_initialized = .true. + + end subroutine mp_tempo_post_init + +!> \section arg_table_mp_tempo_post_run Argument Table +!! \htmlinclude mp_tempo_post_run.html +!! + subroutine mp_tempo_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, ttendlim, & + kdt, mpicomm, mpirank, mpiroot, errmsg, errflg) + + implicit none + + ! Interface variables + integer, intent(in) :: ncol + integer, intent(in) :: nlev + real(kind_phys), dimension(:,:), intent(in) :: tgrs_save + real(kind_phys), dimension(:,:), intent(inout) :: tgrs + real(kind_phys), dimension(:,:), intent(in) :: prslk + real(kind_phys), intent(in) :: dtp + real(kind_phys), intent(in) :: ttendlim + integer, intent(in) :: kdt + ! MPI information + type(MPI_Comm), intent(in ) :: mpicomm + integer, intent(in ) :: mpirank + integer, intent(in ) :: mpiroot + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Local variables + real(kind_phys), dimension(1:ncol,1:nlev) :: mp_tend + integer :: i, k +#ifdef DEBUG + integer :: events +#endif + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Check initialization state + if (.not.is_initialized) then + write(errmsg, fmt='((a))') 'mp_tempo_post_run called before mp_tempo_post_init' + errflg = 1 + return + end if + + ! If limiter is deactivated, return immediately + if (.not.apply_limiter) return + + ! mp_tend and ttendlim are expressed in potential temperature + mp_tend = (tgrs - tgrs_save)/prslk + +#ifdef DEBUG + events = 0 +#endif + do k=1,nlev + do i=1,ncol + mp_tend(i,k) = max( -ttendlim*dtp, min( ttendlim*dtp, mp_tend(i,k) ) ) + +#ifdef DEBUG + if (tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) .ne. tgrs(i,k)) then + write(0,'(a,3i6,3e16.7)') "mp_tempo_post_run mp_tend limiter: kdt, i, k, t_old, t_new, t_lim:", & + & kdt, i, k, tgrs_save(i,k), tgrs(i,k), tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) + events = events + 1 + end if +#endif + tgrs(i,k) = tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) + end do + end do + +#ifdef DEBUG + if (events > 0) then + write(0,'(a,i0,a,i0,a,i0)') "mp_tempo_post_run: ttendlim applied ", events, "/", nlev*ncol, & + & " times at timestep ", kdt + end if +#endif + + end subroutine mp_tempo_post_run + +!! \section arg_table_mp_tempo_post_finalize Argument Table +!! \htmlinclude mp_tempo_post_finalize.html +!! + subroutine mp_tempo_post_finalize(errmsg, errflg) + + implicit none + + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! initialize ccpp error handling variables + errmsg = '' + errflg = 0 + + ! Check initialization state + if (.not. is_initialized) return + + is_initialized = .false. + + end subroutine mp_tempo_post_finalize + +end module mp_tempo_post diff --git a/physics/MP/TEMPO/mp_tempo_post.meta b/physics/MP/TEMPO/mp_tempo_post.meta new file mode 100644 index 000000000..e51edbf27 --- /dev/null +++ b/physics/MP/TEMPO/mp_tempo_post.meta @@ -0,0 +1,154 @@ +[ccpp-table-properties] + name = mp_tempo_post + type = scheme + dependencies = ../../hooks/machine.F + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_post_init + type = scheme +[ttendlim] + standard_name = max_tendency_of_air_potential_temperature_due_to_large_scale_precipitation + long_name = temperature tendency limiter per physics time step + units = K s-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_post_run + type = scheme +[ncol] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[nlev] + standard_name = vertical_layer_dimension + long_name = number of vertical levels + units = count + dimensions = () + type = integer + intent = in +[tgrs_save] + standard_name = air_temperature_save + long_name = air temperature before entering a physics scheme + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[tgrs] + standard_name = air_temperature_of_new_state + long_name = model layer mean temperature + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[prslk] + standard_name = dimensionless_exner_function + long_name = dimensionless Exner function at model layer centers + units = none + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[dtp] + standard_name = timestep_for_physics + long_name = physics timestep + units = s + dimensions = () + type = real + kind = kind_phys + intent = in +[ttendlim] + standard_name = max_tendency_of_air_potential_temperature_due_to_large_scale_precipitation + long_name = temperature tendency limiter per physics time step + units = K s-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[kdt] + standard_name = index_of_timestep + long_name = current forecast iteration + units = index + dimensions = () + type = integer + intent = in +[mpicomm] + standard_name = mpi_communicator + long_name = MPI communicator + units = index + dimensions = () + type = MPI_Comm + intent = in +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_post_finalize + type = scheme +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/MP/TEMPO/mp_tempo_pre.F90 b/physics/MP/TEMPO/mp_tempo_pre.F90 new file mode 100644 index 000000000..1e5b7b92d --- /dev/null +++ b/physics/MP/TEMPO/mp_tempo_pre.F90 @@ -0,0 +1,44 @@ +!>\file mp_tempo_pre.F90 +!! + +! CCPP license goes here, as well as further documentation +!>\ingroup aatempo +module mp_tempo_pre + + use machine, only : kind_phys + + implicit none + + public :: mp_tempo_pre_run + + private + + contains + +!> \section arg_table_mp_tempo_pre_run Argument Table +!! \htmlinclude mp_tempo_pre_run.html +!! + subroutine mp_tempo_pre_run(ncol, nlev, tgrs, tgrs_save, errmsg, errflg) + + implicit none + + ! Interface variables + integer, intent(in ) :: ncol + integer, intent(in ) :: nlev + real(kind_phys), intent(in ) :: tgrs(:,:) + real(kind_phys), intent( out) :: tgrs_save(:,:) + + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Save current air temperature for tendency limiters in mp_tempo_post + tgrs_save = tgrs + + end subroutine mp_tempo_pre_run + +end module mp_tempo_pre diff --git a/physics/MP/TEMPO/mp_tempo_pre.meta b/physics/MP/TEMPO/mp_tempo_pre.meta new file mode 100644 index 000000000..2c6b44c34 --- /dev/null +++ b/physics/MP/TEMPO/mp_tempo_pre.meta @@ -0,0 +1,54 @@ +[ccpp-table-properties] + name = mp_tempo_pre + type = scheme + dependencies = ../../hooks/machine.F + +######################################################################## +[ccpp-arg-table] + name = mp_tempo_pre_run + type = scheme +[ncol] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[nlev] + standard_name = vertical_layer_dimension + long_name = number of vertical levels + units = count + dimensions = () + type = integer + intent = in +[tgrs] + standard_name = air_temperature_of_new_state + long_name = model layer mean temperature + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[tgrs_save] + standard_name = air_temperature_save + long_name = air temperature before entering a physics scheme + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/MP/TEMPO/tempo b/physics/MP/TEMPO/tempo new file mode 160000 index 000000000..2d82021f1 --- /dev/null +++ b/physics/MP/TEMPO/tempo @@ -0,0 +1 @@ +Subproject commit 2d82021f15306512a58efa867a2f076511bccb2a diff --git a/physics/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index 979405cdb..5af54c378 100644 --- a/physics/Radiation/radiation_clouds.f +++ b/physics/Radiation/radiation_clouds.f @@ -3315,7 +3315,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & & modify_qvapor, max_relh, & & kts,kte, debug_flag) ! - USE module_mp_thompson , ONLY : rsif, rslf + USE module_mp_thompson_utils , ONLY : rsif, rslf IMPLICIT NONE ! INTEGER, INTENT(IN):: kts, kte