diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index d2ecef895..84732e401 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -581,7 +581,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input if (Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then do k=1,LMK do i=1,IM - qvs = Statein%qgrs(i,k2,1) + qvs = Statein%qgrs(i,k,1) qv_mp (i,k) = qvs/(1.-qvs) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) @@ -594,7 +594,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input elseif (Model%imp_physics == Model%imp_physics_thompson) then do k=1,LMK do i=1,IM - qvs = Statein%qgrs(i,k2,1) + qvs = Statein%qgrs(i,k,1) qv_mp (i,k) = qvs/(1.-qvs) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) @@ -701,76 +701,60 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input enddo endif elseif (Model%imp_physics == Model%imp_physics_thompson) then ! Thompson MP - if(Model%kdt == 1 ) then - do k=1,lm - k1 = k + kd - do i=1,im - effrl(i,k1) = Tbd%phy_f3d(i,k,Model%nleffr) - effri(i,k1) = Tbd%phy_f3d(i,k,Model%nieffr) - effrr(i,k1) = 1000. ! rrain_def=1000. - effrs(i,k1) = Tbd%phy_f3d(i,k,Model%nseffr) - enddo + ! + ! 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 (Model%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) + 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) + endif + end do + end do + ! Call Thompson's subroutine to compute effective radii + do i=1,im + ! Initialize to default in units m as in module_mp_thompson.F90 + re_cloud(i,:) = 2.49E-6 + re_ice(i,:) = 4.99E-6 + re_snow(i,:) = 9.99E-6 + call calc_effectRad (tlyr(i,:), plyr(i,:), 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 + ! Scale Thompson's effective radii from meter to micron and apply bounds + do k=1,lm + do i=1,im + re_cloud(i,k) = MAX(2.49, MIN(re_cloud(i,k)*1.e6, 50.)) + re_ice(i,k) = MAX(4.99, MIN(re_ice(i,k)*1.e6, 125.)) + !tgs: clduni has different limits for ice radii: 10.0-150.0 + ! it will raise the low limit from 5 to 10, but the + ! high limit will remain 125. + re_snow(i,k) = MAX(9.99, MIN(re_snow(i,k)*1.e6, 999.)) + end do + end do + do k=1,lm + k1 = k + kd + do i=1,im + effrl(i,k1) = re_cloud (i,k) + effri(i,k1) = re_ice (i,k) + effrr(i,k1) = 1000. ! rrain_def=1000. + effrs(i,k1) = re_snow(i,k) enddo - else ! kdt>1 - if(Model%do_mynnedmf .or. & - Model%imfdeepcnv == Model%imfdeepcnv_gf ) then - !tgs - take into account sub-grid clouds from GF or MYNN PBL - - ! Compute effective radii for QC and QI with sub-grid clouds - do k=1,lm - do i=1,im - ! make NC consistent with sub-grid clouds - if (Model%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) - 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) - endif - end do - end do - ! Call Thompson's subroutine to compute effective radii - do i=1,im - ! Initialize to default in units m as in module_mp_thompson.F90 - re_cloud(i,:) = 2.49E-6 - re_ice(i,:) = 4.99E-6 - re_snow(i,:) = 9.99E-6 - call calc_effectRad (tlyr(i,:), plyr(i,:), 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 - do k=1,lm - do i=1,im - re_cloud(i,k) = MAX(2.49, MIN(re_cloud(i,k)*1.e6, 50.)) - re_ice(i,k) = MAX(4.99, MIN(re_ice(i,k)*1.e6, 125.)) - !tgs: clduni has different limits for ice radii: 10.0-150.0 - ! it will raise the low limit from 5 to 10, but the - ! high limit will remain 125. - re_snow(i,k) = MAX(9.99, MIN(re_snow(i,k)*1.e6, 999.)) - end do - end do - - do k=1,lm - k1 = k + kd - do i=1,im - effrl(i,k1) = re_cloud (i,k) ! Tbd%phy_f3d(i,k,Model%nleffr) - effri(i,k1) = re_ice (i,k) ! Tbd%phy_f3d(i,k,Model%nieffr) - effrr(i,k1) = 1000. ! rrain_def=1000. - effrs(i,k1) = Tbd%phy_f3d(i,k,Model%nseffr) - enddo - enddo - else ! not MYNN or not GF - do k=1,lm - k1 = k + kd - do i=1,im - effrl(i,k1) = Tbd%phy_f3d(i,k,Model%nleffr) - effri(i,k1) = Tbd%phy_f3d(i,k,Model%nieffr) - effrr(i,k1) = 1000. ! rrain_def=1000. - effrs(i,k1) = Tbd%phy_f3d(i,k,Model%nseffr) - enddo - enddo - endif ! MYNN PBL or GF conv - endif ! kdt - else ! neither of the other two cases + enddo + ! Update global arrays + do k=1,lm + k1 = k + kd + do i=1,im + Tbd%phy_f3d(i,k,Model%nleffr) = effrl(i,k1) + Tbd%phy_f3d(i,k,Model%nieffr) = effri(i,k1) + Tbd%phy_f3d(i,k,Model%nseffr) = effrs(i,k1) + enddo + enddo + else ! all other cases cldcov = 0.0 endif @@ -936,14 +920,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input else ! kdt > 1 - do k=1,lm - k1 = k + kd - do i=1,im - Tbd%phy_f3d(i,k,Model%nleffr) = effrl(i,k1) - Tbd%phy_f3d(i,k,Model%nieffr) = effri(i,k1) - Tbd%phy_f3d(i,k,Model%nseffr) = effrs(i,k1) - enddo - enddo ! --- call progcld6 to get Xu-Randall total cloud cover (clouds(:,1:LMK,1)) ! tgs: a short subroutine could be made of progcld5 to ! compute only total cloud fraction. diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 3d22cf33b..466bcbb19 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -728,7 +728,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k))/(1.0_kind_phys-spechum(i,k)) !> - Convert number concentration from moist to dry nc_mp(i,k) = gq0(i,k,ntlnc)/(1.0_kind_phys-spechum(i,k)) - nc_mp(i,k) = nc_mp(i,k) + max(0.0, make_DropletNumber(qc_mp(i,k) * rho_dryair(i,k), nwfa(i,k)) * (1.0/rho_dryair(i,k))) + nc_mp(i,k) = max(0.0, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_dryair(i,k), nwfa(i,k)) * (1.0/rho_dryair(i,k))) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k)/(1.0_kind_phys+qv_mp(i,k)) endif @@ -737,7 +737,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))/(1.0_kind_phys-spechum(i,k)) !> - Convert number concentration from moist to dry ni_mp(i,k) = gq0(i,k,ntinc)/(1.0_kind_phys-spechum(i,k)) - ni_mp(i,k) = ni_mp(i,k) + max(0.0, make_IceNumber(qi_mp(i,k) * rho_dryair(i,k), save_tcp(i,k)) * (1.0/rho_dryair(i,k))) + ni_mp(i,k) = max(0.0, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_dryair(i,k), save_tcp(i,k)) * (1.0/rho_dryair(i,k))) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k)/(1.0_kind_phys+qv_mp(i,k)) endif diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 824c4f63c..3f2ee144e 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -67,9 +67,9 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & real(kind_phys), intent(in ) :: phil(:,:) real(kind_phys), intent(in ) :: area(:) ! Cloud effective radii - real(kind_phys), optional, intent(inout) :: re_cloud(:,:) - real(kind_phys), optional, intent(inout) :: re_ice(:,:) - real(kind_phys), optional, intent(inout) :: re_snow(:,:) + real(kind_phys), optional, intent( out) :: re_cloud(:,:) + real(kind_phys), optional, intent( out) :: re_ice(:,:) + real(kind_phys), optional, intent( out) :: re_snow(:,:) ! MPI information integer, intent(in ) :: mpicomm integer, intent(in ) :: mpirank @@ -319,31 +319,40 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & end if ! Calculate initial cloud effective radii if requested - do i = 1, ncol - do k = 1, nlev - re_cloud(i,k) = 2.49E-6 - re_ice(i,k) = 4.99E-6 - re_snow(i,k) = 9.99E-6 + if (present(re_cloud) .and. present(re_ice) .and. present(re_snow)) then + do i = 1, ncol + do k = 1, nlev + re_cloud(i,k) = 2.49E-6 + re_ice(i,k) = 4.99E-6 + re_snow(i,k) = 9.99E-6 + end do + end do + do i = 1, ncol + call calc_effectRad (tgrs(i,:), prsl(i,:), 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, nlev) end do - end do - do i = 1, ncol - call calc_effectRad (tgrs(i,:), prsl(i,:), 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, nlev) - end do - do i = 1, ncol - do k = 1, nlev - re_cloud(i,k) = MAX(2.49E-6, MIN(re_cloud(i,k), 50.E-6)) - re_ice(i,k) = MAX(4.99E-6, MIN(re_ice(i,k), 125.E-6)) - re_snow(i,k) = MAX(9.99E-6, MIN(re_snow(i,k), 999.E-6)) + do i = 1, ncol + do k = 1, nlev + re_cloud(i,k) = MAX(2.49E-6, MIN(re_cloud(i,k), 50.E-6)) + re_ice(i,k) = MAX(4.99E-6, MIN(re_ice(i,k), 125.E-6)) + re_snow(i,k) = MAX(9.99E-6, MIN(re_snow(i,k), 999.E-6)) + end do end do - end do - ! Convert to micron: required for bit-for-bit identical restarts; - ! otherwise entering mp_thompson_init and converting mu to m and - ! back (without updating re_*) introduces b4b differences. - re_cloud = 1.0E6*re_cloud - re_ice = 1.0E6*re_ice - re_snow = 1.0E6*re_snow + !! Convert to micron: required for bit-for-bit identical restarts; + !! otherwise entering mp_thompson_init and converting mu to m and + !! back (without updating re_*) introduces b4b differences. + !! If this code is used, change units in metadata from m to um! + !re_cloud = 1.0E6*re_cloud + !re_ice = 1.0E6*re_ice + !re_snow = 1.0E6*re_snow + else if (present(re_cloud) .or. present(re_ice) .or. present(re_snow)) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_init:', & + ' all or none of the following optional', & + ' arguments are required: re_cloud, re_ice, re_snow' + errflg = 1 + return + end if !> - Convert number concentrations from dry to moist ni = ni_mp/(1.0_kind_phys+qv_mp) diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 81b2241e1..5bbd85732 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -223,7 +223,7 @@ [re_cloud] standard_name = effective_radius_of_stratiform_cloud_liquid_water_particle_in_um long_name = eff. radius of cloud liquid water particle in micrometer - units = um + units = m dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys @@ -232,7 +232,7 @@ [re_ice] standard_name = effective_radius_of_stratiform_cloud_ice_particle_in_um long_name = eff. radius of cloud ice water particle in micrometer - units = um + units = m dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys @@ -241,7 +241,7 @@ [re_snow] standard_name = effective_radius_of_stratiform_cloud_snow_particle_in_um long_name = effective radius of cloud snow particle in micrometer - units = um + units = m dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys