Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove inconsistencies in computation of air density with Thompson MP #79

Merged
merged 8 commits into from
Feb 26, 2021
Merged

Remove inconsistencies in computation of air density with Thompson MP #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5f7968b
1. Fixed units inconsistency for water friendly aerosols in the
tanyasmirnova Feb 23, 2021
a0b31a5
Merge branch 'gsl/develop' of https://github.com/NOAA-GSL/ccpp-physic…
tanyasmirnova Feb 23, 2021
9cd399d
Replace dry-air density to moist-air density in RRTMGP.
tanyasmirnova Feb 23, 2021
a2d800b
Update submodule pointer for physics/rte-rrtmgp
climbfuji Feb 25, 2021
5a9d256
Update submodule pointer for physics/rte-rrtmgp
tanyasmirnova Feb 25, 2021
8449be5
Merge pull request #4 from climbfuji/thomp_nc_density
tanyasmirnova Feb 25, 2021
4c74498
Bugfixes, performance improvements and formatting updates in physics/…
climbfuji Feb 26, 2021
f06e592
Merge pull request #5 from climbfuji/bugfixes_tanya_thomp_nc_density
tanyasmirnova Feb 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions physics/GFS_rrtmg_pre.F90
View file
Open in desktop
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)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are inconsistencies with how I am testing the limits of QC and NC in mp_thompson code. And why are we needing to make droplet number in radiation code? Shouldn't this only be done in one place (MP code) and the resulting number just be given to radiation schemes?

I don't like all the duplication of code for fear that something goes awry/different in one than another place. If we need to make another isolation of ensuring species number and mass jointly make sense, then let's do it once and never more than once. But radiation code seems weird to me calling droplet or ice number fixes.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@joeolson42 @tanyasmirnova @hannahcbarnes was the reason for this the logic in the sgscloud_radpre and sgscloud_radpost schemes, i.e. the modification of qc and qi before radiation with boundary layer clouds from MYNN PBL in sgscloud_radpre, and then restoring the old values in ``sgscloud_radpost`?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have to add the call to make_DropletNumber here only for the case when there are subgrid clouds created in PBL and convective schemes. They create sub-grid scale mixing ratios, but no NCs. After that the calc_effectRad is called using consistent mixing rations and nc_mp and ni_mp. Effective radii are needed for the radiation. After call to radiation, the sub-grid clouds are subtracted to restore the resolved mixing ratios from Thompson MP.

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
Open in desktop
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
Open in desktop
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