From 5f7968b2628b8ab8838d9e9d7cdb4a498f95f009 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Tue, 23 Feb 2021 01:51:19 +0000 Subject: [PATCH 1/5] 1. Fixed units inconsistency for water friendly aerosols in the computation of number concentration of liquid water. 2. Replaced dry-air density with the moist-air density for consistency with Thompson MP. --- physics/GFS_rrtmg_pre.F90 | 10 ++++++---- physics/GFS_suite_interstitial.F90 | 10 +++++----- physics/mp_thompson.F90 | 13 +++++++------ 3 files changed, 18 insertions(+), 15 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 109df3b65..ab488c687 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -288,8 +288,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & plyr(i,k1) = prsl(i,k2) * 0.01 ! pa to mb (hpa) tlyr(i,k1) = tgrs(i,k2) prslk1(i,k1) = prslk(i,k2) - rho(i,k1) = plyr(i,k1)/(con_rd*tlyr(i,k1)) - orho(i,k1) = 1.0/rho(i,k1) !> - Compute relative humidity. es = min( prsl(i,k2), fpvs( tgrs(i,k2) ) ) ! fpvs and prsl in pa @@ -636,6 +634,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) + rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) @@ -649,6 +649,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) + rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) @@ -761,7 +763,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do k=1,lm do i=1,im if (ltaerosol .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)) * orho(i,k) + nc_mp(i,k) = make_DropletNumber(qc_mp(i,k)*rho(i,k), nwfa(i,k)*rho(i,k)) * orho(i,k) 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) @@ -774,7 +776,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & !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,:), qv_mp(i,:), qc_mp(i,:), & + 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,:), & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, lm ) end do diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 551f0e600..f9a26bb57 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -717,7 +717,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ! local variables integer :: i,k,n,tracers - real(kind=kind_phys), dimension(im,levs) :: rho_dryair + real(kind=kind_phys), dimension(im,levs) :: rho_air real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio) @@ -767,16 +767,16 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then do k=1,levs do i=1,im - !> - Density of air in kg m-3 - rho_dryair(i,k) = prsl(i,k) / (con_rd*save_tcp(i,k)) !> - Convert specific humidity to dry mixing ratio qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) + !> - Density of air in kg m-3 + rho_air(i,k) = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) if (ntlnc>0) then !> - Convert moist mixing ratio to dry mixing ratio 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_dryair(i,k), nwfa(i,k)) * (one/rho_dryair(i,k))) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_air(i,k), nwfa(i,k)*rho_air(i,k)) * (one/rho_air(i,k))) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -785,7 +785,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to 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_dryair(i,k), save_tcp(i,k)) * (one/rho_dryair(i,k))) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_air(i,k), save_tcp(i,k)) * (one/rho_air(i,k))) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index ec19945b0..82ddc95be 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -159,10 +159,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Geopotential height in m2 s-2 to height in m hgt = phil/con_g - ! Density of air in kg m-3 and inverse density of air - rho = prsl/(con_rd*tgrs) - orho = 1.0/rho - ! Prior to calling the functions: make_DropletNumber, make_IceNumber, make_RainNumber, ! the incoming mixing ratios should be converted to units of mass/num per cubic meter ! rather than per kg of air. So, to pass back to the model state variables, @@ -178,6 +174,10 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & qs_mp = qs/(1.0_kind_phys-spechum) qg_mp = qg/(1.0_kind_phys-spechum) + ! Density of air in kg m-3 and inverse density of air + rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) + orho = 1.0/rho + !> - Convert number concentrations from moist to dry ni_mp = ni/(1.0_kind_phys-spechum) nr_mp = nr/(1.0_kind_phys-spechum) @@ -304,7 +304,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! If qc is in boundary conditions but nc is not, calculate nc from qc, rho and nwfa if (maxval(qc_mp)>0.0 .and. maxval(nc_mp)==0.0) then - nc_mp = make_DropletNumber(qc_mp*rho, nwfa) * orho + nc_mp = make_DropletNumber(qc_mp*rho, nwfa*rho) * orho end if ! If nc is in boundary conditions but qc is not, reset nc to zero @@ -428,6 +428,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! Air density real(kind_phys) :: rho(1:ncol,1:nlev) !< kg m-3 + !rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) ! Hydrometeors real(kind_phys) :: qv_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) real(kind_phys) :: qc_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) @@ -510,7 +511,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if !> - Density of air in kg m-3 - rho = prsl/(con_rd*tgrs) + rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 w = -omega/(rho*con_g) From 9cd399d66063956ae51c52b89e0ae14972fd9040 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Tue, 23 Feb 2021 17:40:22 +0000 Subject: [PATCH 2/5] Replace dry-air density to moist-air density in RRTMGP. Fix gthe units of water-friendly aerosols in the call to make_DropletNumber. --- physics/GFS_rrtmgp_thompsonmp_pre.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index ea27f3d2b..de1fa3547 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -147,11 +147,11 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do ! Cloud particle sizes and number concentrations... ! First, prepare cloud mixing-ratios and number concentrations for Calc_Re - rho = p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)) - orho = 1./rho do iLay = 1, nLev do iCol = 1, nCol qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) + rho(iCol,iLay) = 0.622*p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)*(qv_mp(iCol,iLay)+0.622)) + orho(iCol,iLay) = 1./rho(iCol,iLay) qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay)) @@ -169,7 +169,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do do iLay = 1, nLev do iCol = 1, nCol if (ltaerosol .and. qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then - nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)) * orho(iCol,iLay) + nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)*rho(iCol,iLay)) * orho(iCol,iLay) endif if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho(iCol,iLay), t_lay(iCol,iLay)) * orho(iCol,iLay) From a2d800b1136af477f3982766e2ee8ac94a2b6fd6 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 25 Feb 2021 12:09:14 -0700 Subject: [PATCH 3/5] Update submodule pointer for physics/rte-rrtmgp --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 566bee9cd..33c8a984c 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 566bee9cd6f9977e82d75d9b4964b20b1ff6163d +Subproject commit 33c8a984c17cf41be5d4c2928242e1b4239bfc40 From 5a9d2567d2cee270cf60e9fcdf35eed64105d4ca Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Thu, 25 Feb 2021 19:17:34 +0000 Subject: [PATCH 4/5] Update submodule pointer for physics/rte-rrtmgp --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 566bee9cd..33c8a984c 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 566bee9cd6f9977e82d75d9b4964b20b1ff6163d +Subproject commit 33c8a984c17cf41be5d4c2928242e1b4239bfc40 From 4c74498e664b912a46df829bb828893a05c12ade Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 26 Feb 2021 07:24:07 -0700 Subject: [PATCH 5/5] Bugfixes, performance improvements and formatting updates in physics/GFS_suite_interstitial.F90 and physics/GFS_rrtmgp_thompsonmp_pre.F90 --- physics/GFS_rrtmgp_thompsonmp_pre.F90 | 140 +++++++++++++------------- physics/GFS_suite_interstitial.F90 | 11 +- 2 files changed, 74 insertions(+), 77 deletions(-) diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index de1fa3547..58e1bddea 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -16,24 +16,24 @@ module GFS_rrtmgp_thompsonmp_pre make_RainNumber use rrtmgp_lw_cloud_optics, only: radliq_lwr, radliq_upr, radice_lwr, radice_upr implicit none - + ! Parameters specific to THOMPSON MP scheme. real(kind_phys), parameter :: & rerain_def = 1000.0 ! Default rain radius to 1000 microns - + public GFS_rrtmgp_thompsonmp_pre_init, GFS_rrtmgp_thompsonmp_pre_run, GFS_rrtmgp_thompsonmp_pre_finalize - -contains + +contains ! ###################################################################################### ! ###################################################################################### subroutine GFS_rrtmgp_thompsonmp_pre_init() end subroutine GFS_rrtmgp_thompsonmp_pre_init - + ! ###################################################################################### ! ###################################################################################### !! \section arg_table_GFS_rrtmgp_thompsonmp_pre_run !! \htmlinclude GFS_rrtmgp_thompsonmp_pre_run.html -!! +!! subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, doLWrad, & i_cldliq, i_cldice, i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, i_cldliq_nc, & i_cldice_nc, i_twa, effr_in, p_lev, p_lay, tv_lay, t_lay, effrin_cldliq, & @@ -42,14 +42,14 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do imfdeepcnv_gf, doGP_cldoptics_PADE, doGP_cldoptics_LUT, & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, & cld_rerain, precip_frac, errmsg, errflg) - - ! Inputs + + ! Inputs integer, intent(in) :: & nCol, & ! Number of horizontal grid points nLev, & ! Number of vertical layers ncnd, & ! Number of cloud condensation types. - nTracers, & ! Number of tracers from model. - i_cldliq, & ! Index into tracer array for cloud liquid amount. + nTracers, & ! Number of tracers from model. + i_cldliq, & ! Index into tracer array for cloud liquid amount. i_cldice, & ! cloud ice amount. i_cldrain, & ! cloud rain amount. i_cldsnow, & ! cloud snow amount. @@ -61,9 +61,9 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do imfdeepcnv, & ! Choice of mass-flux deep convection scheme imfdeepcnv_gf ! Flag for Grell-Freitas deep convection scheme logical, intent(in) :: & - doSWrad, & ! Call SW radiation? - doLWrad, & ! Call LW radiation - effr_in, & ! Use cloud effective radii provided by model? + doSWrad, & ! Call SW radiation? + doLWrad, & ! Call LW radiation + effr_in, & ! Use cloud effective radii provided by model? uni_cld, & ! Use provided cloud-fraction? lmfshal, & ! Flag for mass-flux shallow convection scheme used by Xu-Randall lmfdeep2, & ! Flag for some scale-aware mass-flux convection scheme active @@ -75,7 +75,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do con_g, & ! Physical constant: gravitational constant con_rd ! Physical constant: gas-constant for dry air - real(kind_phys), dimension(nCol,nLev), intent(in) :: & + real(kind_phys), dimension(nCol,nLev), intent(in) :: & tv_lay, & ! Virtual temperature (K) t_lay, & ! Temperature (K) qs_lay, & ! Saturation vapor pressure (Pa) @@ -83,13 +83,13 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do relhum, & ! Relative humidity p_lay, & ! Pressure at model-layers (Pa) cld_frac_mg ! Cloud-fraction from MG scheme. WTF????? - real(kind_phys), dimension(nCol,nLev+1), intent(in) :: & + real(kind_phys), dimension(nCol,nLev+1), intent(in) :: & p_lev ! Pressure at model-level interfaces (Pa) real(kind_phys), dimension(nCol, nLev, nTracers),intent(in) :: & - tracer ! Cloud condensate amount in layer by type () - + tracer ! Cloud condensate amount in layer by type () + ! In/Outs - real(kind_phys), dimension(nCol,nLev), intent(inout) :: & + real(kind_phys), dimension(nCol,nLev), intent(inout) :: & cld_frac, & ! Total cloud fraction cld_lwp, & ! Cloud liquid water path cld_reliq, & ! Cloud liquid effective radius @@ -97,31 +97,32 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do cld_reice, & ! Cloud ice effecive radius cld_swp, & ! Cloud snow water path cld_resnow, & ! Cloud snow effective radius - cld_rwp, & ! Cloud rain water path + cld_rwp, & ! Cloud rain water path cld_rerain, & ! Cloud rain effective radius precip_frac, & ! Precipitation fraction effrin_cldliq, & ! Effective radius for liquid cloud-particles (microns) effrin_cldice, & ! Effective radius for ice cloud-particles (microns) - effrin_cldsnow ! Effective radius for snow cloud-particles (microns) - - ! Outputs + effrin_cldsnow ! Effective radius for snow cloud-particles (microns) + + ! Outputs character(len=*), intent(out) :: & errmsg ! Error message - integer, intent(out) :: & + integer, intent(out) :: & errflg ! Error flag - + ! Local variables real(kind_phys) :: alpha0, pfac, tem1, cld_mr real(kind_phys), dimension(nCol, nLev, min(4,ncnd)) :: cld_condensate integer :: iCol,iLay,l - real(kind_phys), dimension(nCol,nLev) :: deltaP, deltaZ, rho, orho, re_cloud, re_ice,& + real(kind_phys) :: rho, orho + real(kind_phys), dimension(nCol,nLev) :: deltaP, deltaZ, re_cloud, re_ice,& re_snow, qv_mp, qc_mp, qi_mp, qs_mp, nc_mp, ni_mp, nwfa logical :: top_at_1 - + ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - + if (.not. (doSWrad .or. doLWrad)) return ! Cloud condensate @@ -129,29 +130,30 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do cld_condensate(1:nCol,1:nLev,2) = tracer(1:nCol,1:nLev,i_cldice) ! -ice water cld_condensate(1:nCol,1:nLev,3) = tracer(1:nCol,1:nLev,i_cldrain) ! -rain water cld_condensate(1:nCol,1:nLev,4) = tracer(1:nCol,1:nLev,i_cldsnow) + &! -snow + grapuel - tracer(1:nCol,1:nLev,i_cldgrpl) - + tracer(1:nCol,1:nLev,i_cldgrpl) + ! Cloud water path (g/m2) - deltaP = abs(p_lev(:,2:nLev+1)-p_lev(:,1:nLev))/100. + deltaP = abs(p_lev(:,2:nLev+1)-p_lev(:,1:nLev))/100. do iLay = 1, nLev do iCol = 1, nCol - ! Compute liquid/ice condensate path from mixing ratios (kg/kg)->(g/m2) + ! Compute liquid/ice condensate path from mixing ratios (kg/kg)->(g/m2) tem1 = (1.0e5/con_g) * deltaP(iCol,iLay) cld_lwp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,1) * tem1) cld_iwp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,2) * tem1) cld_rwp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,3) * tem1) - cld_swp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,4) * tem1) + cld_swp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,4) * tem1) enddo - enddo - + enddo + ! Cloud particle sizes and number concentrations... - - ! First, prepare cloud mixing-ratios and number concentrations for Calc_Re + + ! Prepare cloud mixing-ratios and number concentrations for calc_effectRad, + ! and update number concentrations, consistent with sub-grid clouds do iLay = 1, nLev - do iCol = 1, nCol + do iCol = 1, nCol qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) - rho(iCol,iLay) = 0.622*p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)*(qv_mp(iCol,iLay)+0.622)) - orho(iCol,iLay) = 1./rho(iCol,iLay) + rho = 0.622*p_lay(iCol,iLay)/(con_rd*t_lay(iCol,iLay)*(qv_mp(iCol,iLay)+0.622)) + orho = 1./rho qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay)) @@ -159,24 +161,18 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do ni_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice_nc) / (1.-q_lay(iCol,iLay)) if (ltaerosol) then nwfa(iCol,iLay) = tracer(iCol,iLay,i_twa) + if (qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then + nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho + endif else - nc_mp(iCol,iLay) = nt_c*orho(iCol,iLay) - endif - enddo - enddo - - ! Update number concentration, consistent with sub-grid clouds - do iLay = 1, nLev - do iCol = 1, nCol - if (ltaerosol .and. qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then - nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)*rho(iCol,iLay)) * orho(iCol,iLay) + nc_mp(iCol,iLay) = nt_c*orho endif if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then - ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho(iCol,iLay), t_lay(iCol,iLay)) * orho(iCol,iLay) + ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho, t_lay(iCol,iLay)) * orho endif enddo enddo - + ! Compute effective radii for liquid/ice/snow using subgrid scale clouds ! Call Thompson's subroutine to compute effective radii do iCol=1,nCol @@ -184,14 +180,14 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:), & re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), 1, nLev ) enddo - - ! Scale Thompson's effective radii from meter to micron + + ! Scale Thompson's effective radii from meter to micron effrin_cldliq(1:nCol,1:nLev) = re_cloud(1:nCol,1:nLev)*1.e6 effrin_cldice(1:nCol,1:nLev) = re_ice(1:nCol,1:nLev)*1.e6 effrin_cldsnow(1:nCol,1:nLev) = re_snow(1:nCol,1:nLev)*1.e6 - - ! Bound effective radii for RRTMGP, LUT's for cloud-optics go from - ! 2.5 - 21.5 microns for liquid clouds, + + ! Bound effective radii for RRTMGP, LUT's for cloud-optics go from + ! 2.5 - 21.5 microns for liquid clouds, ! 10 - 180 microns for ice-clouds if (doGP_cldoptics_PADE .or. doGP_cldoptics_LUT) then where(effrin_cldliq .lt. radliq_lwr) effrin_cldliq = radliq_lwr @@ -205,32 +201,32 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do cld_reice(1:nCol,1:nLev) = effrin_cldice(1:nCol,1:nLev) cld_resnow(1:nCol,1:nLev) = effrin_cldsnow(1:nCol,1:nLev) cld_rerain(1:nCol,1:nLev) = rerain_def - + ! Compute cloud-fraction. Else, use value provided - if(.not. do_mynnedmf .or. imfdeepcnv .ne. imfdeepcnv_gf ) then ! MYNN PBL or GF conv + if(.not. do_mynnedmf .or. imfdeepcnv .ne. imfdeepcnv_gf ) then ! MYNN PBL or GF conv ! Cloud-fraction if (uni_cld) then - cld_frac(1:nCol,1:nLev) = cld_frac_mg(1:nCol,1:nLev) + cld_frac(1:nCol,1:nLev) = cld_frac_mg(1:nCol,1:nLev) else if( lmfshal) alpha0 = 100. ! Default (from GATE simulations) if(.not. lmfshal) alpha0 = 2000. - ! Xu-Randall (1996) cloud-fraction + ! Xu-Randall (1996) cloud-fraction do iLay = 1, nLev do iCol = 1, nCol cld_mr = cld_condensate(iCol,iLay,1) + cld_condensate(iCol,iLay,2) + & cld_condensate(iCol,iLay,4) - cld_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & - qs_lay(iCol,iLay), relhum(iCol,iLay), cld_mr, alpha0) + cld_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & + qs_lay(iCol,iLay), relhum(iCol,iLay), cld_mr, alpha0) enddo enddo endif endif - + ! Precipitation fraction (Hack. For now use cloud-fraction) precip_frac(1:nCol,1:nLev) = cld_frac(1:nCol,1:nLev) - + end subroutine GFS_rrtmgp_thompsonmp_pre_run - + ! ###################################################################################### ! ###################################################################################### subroutine GFS_rrtmgp_thompsonmp_pre_finalize() @@ -241,11 +237,11 @@ end subroutine GFS_rrtmgp_thompsonmp_pre_finalize ! Xu-Randall(1996) A Semiempirical Cloudiness Parameterization for Use in Climate Models ! https://doi.org/10.1175/1520-0469(1996)053<3084:ASCPFU>2.0.CO;2 ! - ! cld_frac = {1-exp[-alpha*cld_mr/((1-relhum)*qs_lay)**lambda]}*relhum**P + ! cld_frac = {1-exp[-alpha*cld_mr/((1-relhum)*qs_lay)**lambda]}*relhum**P ! ! ###################################################################################### function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) - + ! Inputs real(kind_phys), intent(in) :: & p_lay, & ! Pressure (Pa) @@ -253,7 +249,7 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) relhum, & ! Relative humidity cld_mr, & ! Total cloud mixing ratio alpha ! Scheme parameter (default=100) - + ! Outputs real(kind_phys) :: cld_frac_XuRandall @@ -262,21 +258,21 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) ! Parameters real(kind_phys) :: & - lambda = 0.50, & ! + lambda = 0.50, & ! P = 0.25 - + clwt = 1.0e-6 * (p_lay*0.001) if (cld_mr > clwt) then onemrh = max(1.e-10, 1.0 - relhum) tem1 = alpha / min(max((onemrh*qs_lay)**lambda,0.0001),1.0) tem2 = max(min(tem1*(cld_mr - clwt), 50.0 ), 0.0 ) tem3 = sqrt(sqrt(relhum)) ! This assumes "p" = 0.25. Identical, but cheaper than relhum**p - ! + ! cld_frac_XuRandall = max( tem3*(1.0-exp(-tem2)), 0.0 ) else cld_frac_XuRandall = 0.0 endif - + return end function end module GFS_rrtmgp_thompsonmp_pre diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 464ede8cc..b51b6f9c8 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -692,7 +692,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ! local variables integer :: i,k,n,tracers - real(kind=kind_phys), dimension(im,levs) :: rho_air + real(kind=kind_phys) :: rho, orho real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio) @@ -744,14 +744,15 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to do i=1,im !> - Convert specific humidity to dry mixing ratio qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) - !> - Density of air in kg m-3 - rho_air(i,k) = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) + !> - Density of air in kg m-3 and inverse density + rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) + orho = one/rho if (ntlnc>0) then !> - Convert moist mixing ratio to dry mixing ratio 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_air(i,k), nwfa(i,k)*rho_air(i,k)) * (one/rho_air(i,k))) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -760,7 +761,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to 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_air(i,k), save_tcp(i,k)) * (one/rho_air(i,k))) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif