Skip to content

Commit

Permalink
Merge pull request NCAR#30 from climbfuji/thompson_mp_cloud_effective…
Browse files Browse the repository at this point in the history 
…_radii_cleanup

Cleanup of Thompson MP cloud effective radii calculation
  • Loading branch information
@DomHeinzeller
DomHeinzeller authored Jun 3, 2020
2 parents 64e8ff6 + 2aa91ae commit 02bdecf
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 110 deletions.
134 changes: 55 additions & 79 deletions physics/GFS_rrtmg_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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.
Expand Down
4 changes: 2 additions & 2 deletions physics/GFS_suite_interstitial.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
61 changes: 35 additions & 26 deletions physics/mp_thompson.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions physics/mp_thompson.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 02bdecf

Please sign in to comment.