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

Modify the Thompson scheme to improve radiative fluxes and cloud cover for HR1 #19

Merged
merged 3 commits into from
Nov 8, 2022
Merged

Modify the Thompson scheme to improve radiative fluxes and cloud cover for HR1 #19

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e14df8c
Modify the Thompson and RRTMG(P) schemes to improve radiative fluxes …
ChunxiZhang-NOAA Nov 1, 2022
729b6e6
Remove the lines that are commented out in GFS_rrtmg_pre.F90 and modu…
ChunxiZhang-NOAA Nov 4, 2022
e947b08
Make changes to the Thompson scheme to address reviewer comments
ChunxiZhang-NOAA Nov 7, 2022
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
13 changes: 10 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 @@ -73,7 +73,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
use surface_perturbation, only: cdfnor,ppfbet

! For Thompson MP
use module_mp_thompson, only: calc_effectRad, Nt_c, &
use module_mp_thompson, only: calc_effectRad, &
Nt_c_l, Nt_c_o, &
re_qc_min, re_qc_max, &
re_qi_min, re_qi_max, &
re_qs_min, re_qs_max
Expand Down Expand Up @@ -245,6 +246,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
real (kind=kind_phys) :: alpha0,beta0,m,s,cldtmp,tmp_wt,cdfz
real (kind=kind_phys) :: max_relh
integer :: iflag
integer :: islmsk

integer :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
Expand Down Expand Up @@ -748,7 +750,11 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
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)
nc_mp (i,k) = nt_c*orho(i,k)
if(nint(slmsk(i)) == 1) then
nc_mp (i,k) = Nt_c_l*orho(i,k)
else
nc_mp (i,k) = Nt_c_o*orho(i,k)
endif
ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs)
enddo
enddo
Expand Down Expand Up @@ -877,13 +883,14 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
end do
!> - Call Thompson's subroutine calc_effectRad() to compute effective radii
do i=1,im
islmsk = nint(slmsk(i))
! Effective radii [m] are now intent(out), bounds applied in calc_effectRad
!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,:)*100., qv_mp(i,:), qc_mp(i,:), &
nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), &
effrl(i,:), effri(i,:), effrs(i,:), 1, lm )
effrl(i,:), effri(i,:), effrs(i,:), islmsk, 1, lm )
! Scale Thompson's effective radii from meter to micron
do k=1,lm
effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max))*1.e6
Expand Down
20 changes: 14 additions & 6 deletions physics/GFS_rrtmgp_cloud_mp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ module GFS_rrtmgp_cloud_mp
use rrtmgp_lw_cloud_optics, only: &
radliq_lwr => radliq_lwrLW, radliq_upr => radliq_uprLW,&
radice_lwr => radice_lwrLW, radice_upr => radice_uprLW
use module_mp_thompson, only: calc_effectRad, Nt_c, re_qc_min, re_qc_max, re_qi_min, &
re_qi_max, re_qs_min, re_qs_max
use module_mp_thompson, only: calc_effectRad, Nt_c_l, Nt_c_o, re_qc_min, re_qc_max, &
re_qi_min, re_qi_max, re_qs_min, re_qs_max
use module_mp_thompson_make_number_concentrations, only: make_IceNumber, &
make_DropletNumber, make_RainNumber

Expand Down Expand Up @@ -254,7 +254,7 @@ subroutine GFS_rrtmgp_cloud_mp_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldic
! Update particle size using modified mixing-ratios from Thompson.
call cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, &
i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol,&
mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow)
mraerosol, lsmask, effrin_cldliq, effrin_cldice, effrin_cldsnow)
cld_reliq = effrin_cldliq
cld_reice = effrin_cldice
cld_resnow = effrin_cldsnow
Expand Down Expand Up @@ -820,7 +820,7 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha)
!! \section cmp_reff_Thompson_gen General Algorithm
subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, &
i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol, &
mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow)
mraerosol, lsmask, effrin_cldliq, effrin_cldice, effrin_cldsnow)
implicit none

! Inputs
Expand All @@ -830,6 +830,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
real(kind_phys), intent(in) :: con_eps,con_rd
real(kind_phys), dimension(:,:),intent(in) :: q_lay, p_lay, t_lay
real(kind_phys), dimension(:,:,:),intent(in) :: tracer
real(kind_phys), dimension(:), intent(in) :: lsmask

! Outputs
real(kind_phys), dimension(:,:), intent(inout) :: effrin_cldliq, effrin_cldice, &
Expand All @@ -840,6 +841,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
real(kind_phys) :: rho, orho
real(kind_phys),dimension(nCol,nLev) :: qv_mp, qc_mp, qi_mp, qs_mp, ni_mp, nc_mp, &
nwfa, re_cloud, re_ice, re_snow
integer :: ilsmask

! Prepare cloud mixing-ratios and number concentrations for calc_effectRa
do iLay = 1, nLev
Expand All @@ -863,7 +865,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
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
if (nint(lsmask(iCol)) == 1) then !land
nc_mp(iCol,iLay) = nt_c_l*orho
else
nc_mp(iCol,iLay) = nt_c_o*orho
endif
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, t_lay(iCol,iLay)) * orho
Expand All @@ -873,9 +879,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice

! Compute effective radii for liquid/ice/snow.
do iCol=1,nCol
ilsmask = nint(lsmask(iCol))
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 )
re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), ilsmask, &
1, nLev )
do iLay = 1, nLev
re_cloud(iCol,iLay) = MAX(re_qc_min, MIN(re_cloud(iCol,iLay), re_qc_max))
re_ice(iCol,iLay) = MAX(re_qi_min, MIN(re_ice(iCol,iLay), re_qi_max))
Expand Down
82 changes: 62 additions & 20 deletions physics/module_mp_thompson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,9 @@ MODULE module_mp_thompson
!.. scheme. In 2-moment cloud water, Nt_c represents a maximum of
!.. droplet concentration and nu_c is also variable depending on local
!.. droplet number concentration.
REAL, PARAMETER :: Nt_c = 100.E6
!REAL, PARAMETER :: Nt_c = 100.E6
ChunxiZhang-NOAA marked this conversation as resolved.
Show resolved Hide resolved
REAL, PARAMETER :: Nt_c_o = 50.E6
ChunxiZhang-NOAA marked this conversation as resolved.
Show resolved Hide resolved
REAL, PARAMETER :: Nt_c_l = 100.E6
REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6

!..Declaration of constants for assumed CCN/IN aerosols when none in
Expand All @@ -108,7 +110,7 @@ MODULE module_mp_thompson
REAL, PARAMETER, PRIVATE:: mu_r = 0.0
REAL, PARAMETER, PRIVATE:: mu_g = 0.0
REAL, PARAMETER, PRIVATE:: mu_i = 0.0
REAL, PRIVATE:: mu_c
REAL, PRIVATE:: mu_c_o, mu_c_l

!..Sum of two gamma distrib for snow (Field et al. 2005).
!.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
Expand Down Expand Up @@ -150,7 +152,6 @@ MODULE module_mp_thompson
REAL, PARAMETER, PRIVATE:: fv_s = 100.0
REAL, PARAMETER, PRIVATE:: av_g = 442.0
REAL, PARAMETER, PRIVATE:: bv_g = 0.89
REAL, PARAMETER, PRIVATE:: av_i = 1493.9
REAL, PARAMETER, PRIVATE:: bv_i = 1.0
REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8
REAL, PARAMETER, PRIVATE:: bv_c = 2.0
Expand Down Expand Up @@ -534,7 +535,8 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
!.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
!.. to 2 for really dirty air. This not used in 2-moment cloud water
!.. scheme and nu_c used instead and varies from 2 to 15 (integer-only).
mu_c = MIN(15., (1000.E6/Nt_c + 2.))
mu_c_l = MIN(15., (1000.E6/Nt_c_l + 2.))
mu_c_o = MIN(15., (1000.E6/Nt_c_o + 2.))

!> - Compute Schmidt number to one-third used numerous times
Sc3 = Sc**(1./3.)
Expand Down Expand Up @@ -889,7 +891,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &

if (mpirank==mpiroot) write (*,*)'creating microphysics lookup tables ... '
if (mpirank==mpiroot) write (*,'(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
' using: mu_c_o=',mu_c_o,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g

!> - Call table_ccnact() to read a static file containing CCN activation of aerosols. The
!! data were created from a parcel model by Feingold & Heymsfield with
Expand Down Expand Up @@ -982,7 +984,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
nwfa, nifa, nwfa2d, nifa2d, &
tt, th, pii, &
p, w, dz, dt_in, dt_inner, &
sedi_semi, decfl, &
sedi_semi, decfl, lsm, &
RAINNC, RAINNCV, &
SNOWNC, SNOWNCV, &
ICENC, ICENCV, &
Expand Down Expand Up @@ -1037,6 +1039,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
nc, nwfa, nifa
REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d
INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN):: lsm
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
re_cloud, re_ice, re_snow
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: pfils, pflls
Expand Down Expand Up @@ -1117,6 +1120,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
REAL:: dt, pptrain, pptsnow, pptgraul, pptice
REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
INTEGER:: lsml
REAL:: rand1, rand2, rand3, rand_pert_max
INTEGER:: i, j, k, m
INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
Expand Down Expand Up @@ -1419,8 +1423,13 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
nifa1d(k) = nifa(i,k,j)
enddo
else
lsml = lsm(i,j)
do k = kts, kte
nc1d(k) = Nt_c/rho(k)
if(lsml == 1) then
nc1d(k) = Nt_c_l/rho(k)
else
nc1d(k) = Nt_c_o/rho(k)
endif
nwfa1d(k) = 11.1E6
nifa1d(k) = naIN1*0.01
enddo
Expand All @@ -1429,7 +1438,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
!> - Call mp_thompson()
call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, &
pptrain, pptsnow, pptgraul, pptice, &
lsml, pptrain, pptsnow, pptgraul, pptice, &
#if ( WRF_CHEM == 1 )
rainprod1d, evapprod1d, &
#endif
Expand Down Expand Up @@ -1698,7 +1707,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
enddo
!> - Call calc_effectrad()
call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
re_qc1d, re_qi1d, re_qs1d, kts, kte)
re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
do k = kts, kte
re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max))
re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max))
Expand Down Expand Up @@ -1841,7 +1850,7 @@ END SUBROUTINE thompson_finalize
!> @{
subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, &
pptrain, pptsnow, pptgraul, pptice, &
lsml, pptrain, pptsnow, pptgraul, pptice, &
#if ( WRF_CHEM == 1 )
rainprod, evapprod, &
#endif
Expand Down Expand Up @@ -1879,6 +1888,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq
REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
REAL, INTENT(IN):: dt
INTEGER, INTENT(IN):: lsml
REAL, INTENT(IN):: rand1, rand2, rand3
! Extended diagnostics, most arrays only allocated if ext_diag is true
LOGICAL, INTENT(IN) :: ext_diag
Expand Down Expand Up @@ -1982,6 +1992,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
REAL:: Ef_ra, Ef_sa, Ef_ga
REAL:: dtsave, odts, odt, odzq, hgt_agl, SR
REAL:: xslw1, ygra1, zans1, eva_factor
REAL:: av_i
INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
INTEGER, DIMENSION(5):: ksed1
INTEGER:: nir, nis, nig, nii, nic, niin
Expand All @@ -2006,6 +2017,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
odt = 1./dt
odts = 1./dtsave
iexfrq = 1
! Transition value of coefficient matching at crossover from cloud ice to snow
av_i = av_s * D0s ** (bv_s - bv_i)
ChunxiZhang-NOAA marked this conversation as resolved.
Show resolved Hide resolved

!+---+-----------------------------------------------------------------+
!> - Initialize Source/sink terms. First 2 chars: "pr" represents source/sink of
Expand Down Expand Up @@ -2210,7 +2223,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
endif
nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
/ am_r*lamc**bm_r)
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if (lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
else
qc1d(k) = 0.0
nc1d(k) = 0.0
Expand Down Expand Up @@ -2919,13 +2938,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &

!> - Deposition nucleation of dust/mineral from DeMott et al (2010)
!! we may need to relax the temperature and ssati constraints.
if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps &
if ( (ssati(k).ge. 0.15) .or. (ssatw(k).gt. eps &
.and. temp(k).lt.253.15) ) then
if (dustyIce .AND. (is_aerosol_aware .or. merra2_aerosol_aware)) then
xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
xnc = xnc*(1.0 + 50.*rand3)
else
xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
xnc = MIN(1000.E3, TNO*EXP(ATO*(T_0-temp(k))))
endif
xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
Expand Down Expand Up @@ -3273,7 +3292,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
lami = cie(2)/5.E-6
xni = MIN(4999.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
elseif (xDi.gt. 300.E-6) then
elseif (xDi.gt. 300.E-6) then
lami = cie(2)/300.E-6
xni = cig(1)*oig2*xri/am_i*lami**bm_i
niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
Expand Down Expand Up @@ -3389,7 +3408,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if(lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
L_qc(k) = .true.
else
rc(k) = R1
Expand Down Expand Up @@ -3560,7 +3585,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
if (is_aerosol_aware .or. merra2_aerosol_aware) then
xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k)))
else
xnc = Nt_c
if(lsml == 1) then
xnc = Nt_c_l
else
xnc = Nt_c_o
endif
endif
pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho

Expand Down Expand Up @@ -3630,7 +3659,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
if (rc(k).eq.R1) L_qc(k) = .false.
nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if(lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
temp(k) = t1d(k) + DT*tten(k)
rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
Expand Down Expand Up @@ -4235,7 +4270,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
xDi = (bm_i + mu_i + 1.) * ilami
if (xDi.lt. 5.E-6) then
lami = cie(2)/5.E-6
elseif (xDi.gt. 300.E-6) then
elseif (xDi.gt. 300.E-6) then
lami = cie(2)/300.E-6
endif
ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
Expand Down Expand Up @@ -5749,7 +5784,7 @@ END FUNCTION delta_p
!! distribution, not the second part, which is the larger sizes.

subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
& re_qc1d, re_qi1d, re_qs1d, kts, kte)
& re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)

IMPLICIT NONE

Expand All @@ -5766,6 +5801,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
DOUBLE PRECISION:: lamc, lami
LOGICAL:: has_qc, has_qi, has_qs
INTEGER:: inu_c
INTEGER:: lsml
real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
& 504,720,990,1320,1716,2184,2730,3360,4080,4896/)

Expand All @@ -5781,7 +5817,13 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
rc(k) = MAX(R1, qc1d(k)*rho(k))
nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
if( lsml == 1) then
nc(k) = Nt_c_l
else
nc(k) = Nt_c_o
endif
endif
if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true.
ri(k) = MAX(R1, qi1d(k)*rho(k))
ni(k) = MAX(R2, ni1d(k)*rho(k))
Expand Down
Loading