Skip to content

Commit

Permalink
Merge pull request #79 from tanyasmirnova/thomp_nc_density
Browse files Browse the repository at this point in the history 
Remove inconsistencies in computation of air density with Thompson MP
  • Loading branch information
@DomHeinzeller
DomHeinzeller authored Feb 26, 2021
2 parents b2c7bd5 + f06e592 commit c82e501
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 86 deletions.
8 changes: 5 additions & 3 deletions physics/GFS_rrtmg_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) = prsl(i,k2)/(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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
140 changes: 68 additions & 72 deletions physics/GFS_rrtmgp_thompsonmp_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -75,123 +75,119 @@ 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)
q_lay, & ! water-vapor mixing ratio (kg/kg)
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
cld_iwp, & ! Cloud ice water path
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
cld_condensate(1:nCol,1:nLev,1) = tracer(1:nCol,1:nLev,i_cldliq) ! -liquid water
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
rho = p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev))
orho = 1./rho

! 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 = 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))
nc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq_nc) / (1.-q_lay(iCol,iLay))
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)) * 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
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,:), 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
Expand All @@ -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()
Expand All @@ -241,19 +237,19 @@ 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)
qs_lay, & ! Saturation vapor-pressure (Pa)
relhum, & ! Relative humidity
cld_mr, & ! Total cloud mixing ratio
alpha ! Scheme parameter (default=100)

! Outputs
real(kind_phys) :: cld_frac_XuRandall

Expand All @@ -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
11 changes: 6 additions & 5 deletions physics/GFS_suite_interstitial.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_dryair
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)
Expand Down Expand Up @@ -742,16 +742,17 @@ 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 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_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, 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
Expand All @@ -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_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, 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
Expand Down
Loading

0 comments on commit c82e501

Please sign in to comment.