Skip to content

Commit e14df8c

Browse files
ChunxiZhang-NOAAChunxiZhang-NOAA
committed
Modify the Thompson and RRTMG(P) schemes to improve radiative fluxes and cloud cover
1 parent aa7bff5 commit e14df8c

5 files changed

+106
-35
lines changed

‎physics/GFS_rrtmg_pre.F90

+11-3
Original file line numberDiff line numberDiff line change
@@ -73,7 +73,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
7373
use surface_perturbation, only: cdfnor,ppfbet
7474

7575
! For Thompson MP
76-
use module_mp_thompson, only: calc_effectRad, Nt_c, &
76+
use module_mp_thompson, only: calc_effectRad, &
77+
Nt_c_l, Nt_c_o, &
7778
re_qc_min, re_qc_max, &
7879
re_qi_min, re_qi_max, &
7980
re_qs_min, re_qs_max
@@ -245,6 +246,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
245246
real (kind=kind_phys) :: alpha0,beta0,m,s,cldtmp,tmp_wt,cdfz
246247
real (kind=kind_phys) :: max_relh
247248
integer :: iflag
249+
integer :: islmsk
248250

249251
integer :: ids, ide, jds, jde, kds, kde, &
250252
ims, ime, jms, jme, kms, kme, &
@@ -748,7 +750,12 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
748750
qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs)
749751
qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs)
750752
qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs)
751-
nc_mp (i,k) = nt_c*orho(i,k)
753+
!nc_mp (i,k) = nt_c*orho(i,k)
754+
if(slmsk(i) == 0.) then
755+
nc_mp (i,k) = Nt_c_o*orho(i,k)
756+
else
757+
nc_mp (i,k) = Nt_c_l*orho(i,k)
758+
endif
752759
ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs)
753760
enddo
754761
enddo
@@ -877,13 +884,14 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, &
877884
end do
878885
!> - Call Thompson's subroutine calc_effectRad() to compute effective radii
879886
do i=1,im
887+
islmsk = nint(slmsk(i))
880888
! Effective radii [m] are now intent(out), bounds applied in calc_effectRad
881889
!tgs: progclduni has different limits for ice radii (10.0-150.0) than
882890
! calc_effectRad (4.99-125.0 for WRFv3.8.1; 2.49-125.0 for WRFv4+)
883891
! it will raise the low limit from 5 to 10, but the high limit will remain 125.
884892
call calc_effectRad (tlyr(i,:), plyr(i,:)*100., qv_mp(i,:), qc_mp(i,:), &
885893
nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), &
886-
effrl(i,:), effri(i,:), effrs(i,:), 1, lm )
894+
effrl(i,:), effri(i,:), effrs(i,:), islmsk, 1, lm )
887895
! Scale Thompson's effective radii from meter to micron
888896
do k=1,lm
889897
effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max))*1.e6

‎physics/GFS_rrtmgp_cloud_mp.F90

+14-6
Original file line numberDiff line numberDiff line change
@@ -12,8 +12,8 @@ module GFS_rrtmgp_cloud_mp
1212
use rrtmgp_lw_cloud_optics, only: &
1313
radliq_lwr => radliq_lwrLW, radliq_upr => radliq_uprLW,&
1414
radice_lwr => radice_lwrLW, radice_upr => radice_uprLW
15-
use module_mp_thompson, only: calc_effectRad, Nt_c, re_qc_min, re_qc_max, re_qi_min, &
16-
re_qi_max, re_qs_min, re_qs_max
15+
use module_mp_thompson, only: calc_effectRad, Nt_c_l, Nt_c_o, re_qc_min, re_qc_max, &
16+
re_qi_min, re_qi_max, re_qs_min, re_qs_max
1717
use module_mp_thompson_make_number_concentrations, only: make_IceNumber, &
1818
make_DropletNumber, make_RainNumber
1919

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

826826
! Inputs
@@ -830,6 +830,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
830830
real(kind_phys), intent(in) :: con_eps,con_rd
831831
real(kind_phys), dimension(:,:),intent(in) :: q_lay, p_lay, t_lay
832832
real(kind_phys), dimension(:,:,:),intent(in) :: tracer
833+
real(kind_phys), dimension(:), intent(in) :: lsmask
833834

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

844846
! Prepare cloud mixing-ratios and number concentrations for calc_effectRa
845847
do iLay = 1, nLev
@@ -863,7 +865,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
863865
nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho
864866
endif
865867
else
866-
nc_mp(iCol,iLay) = nt_c*orho
868+
if (nint(lsmask(iCol)) == 0) then !land
869+
nc_mp(iCol,iLay) = nt_c_o*orho
870+
else
871+
nc_mp(iCol,iLay) = nt_c_l*orho
872+
endif
867873
endif
868874
if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then
869875
ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho, t_lay(iCol,iLay)) * orho
@@ -873,9 +879,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
873879

874880
! Compute effective radii for liquid/ice/snow.
875881
do iCol=1,nCol
882+
ilsmask = nint(lsmask(iCol))
876883
call calc_effectRad (t_lay(iCol,:), p_lay(iCol,:), qv_mp(iCol,:), qc_mp(iCol,:), &
877884
nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:), &
878-
re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), 1, nLev )
885+
re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), ilsmask, &
886+
1, nLev )
879887
do iLay = 1, nLev
880888
re_cloud(iCol,iLay) = MAX(re_qc_min, MIN(re_cloud(iCol,iLay), re_qc_max))
881889
re_ice(iCol,iLay) = MAX(re_qi_min, MIN(re_ice(iCol,iLay), re_qi_max))

‎physics/module_mp_thompson.F90

+67-21
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,9 @@ MODULE module_mp_thompson
9292
!.. scheme. In 2-moment cloud water, Nt_c represents a maximum of
9393
!.. droplet concentration and nu_c is also variable depending on local
9494
!.. droplet number concentration.
95-
REAL, PARAMETER :: Nt_c = 100.E6
95+
!REAL, PARAMETER :: Nt_c = 100.E6
96+
REAL, PARAMETER :: Nt_c_o = 50.E6
97+
REAL, PARAMETER :: Nt_c_l = 100.E6
9698
REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6
9799

98100
!..Declaration of constants for assumed CCN/IN aerosols when none in
@@ -108,7 +110,8 @@ MODULE module_mp_thompson
108110
REAL, PARAMETER, PRIVATE:: mu_r = 0.0
109111
REAL, PARAMETER, PRIVATE:: mu_g = 0.0
110112
REAL, PARAMETER, PRIVATE:: mu_i = 0.0
111-
REAL, PRIVATE:: mu_c
113+
!REAL, PRIVATE:: mu_c
114+
REAL, PRIVATE:: mu_c_o, mu_c_l
112115

113116
!..Sum of two gamma distrib for snow (Field et al. 2005).
114117
!.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
@@ -150,7 +153,7 @@ MODULE module_mp_thompson
150153
REAL, PARAMETER, PRIVATE:: fv_s = 100.0
151154
REAL, PARAMETER, PRIVATE:: av_g = 442.0
152155
REAL, PARAMETER, PRIVATE:: bv_g = 0.89
153-
REAL, PARAMETER, PRIVATE:: av_i = 1493.9
156+
!REAL, PARAMETER, PRIVATE:: av_i = 1493.9
154157
REAL, PARAMETER, PRIVATE:: bv_i = 1.0
155158
REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8
156159
REAL, PARAMETER, PRIVATE:: bv_c = 2.0
@@ -534,7 +537,9 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
534537
!.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
535538
!.. to 2 for really dirty air. This not used in 2-moment cloud water
536539
!.. scheme and nu_c used instead and varies from 2 to 15 (integer-only).
537-
mu_c = MIN(15., (1000.E6/Nt_c + 2.))
540+
!mu_c = MIN(15., (1000.E6/Nt_c + 2.))
541+
mu_c_l = MIN(15., (1000.E6/Nt_c_l + 2.))
542+
mu_c_o = MIN(15., (1000.E6/Nt_c_o + 2.))
538543

539544
!> - Compute Schmidt number to one-third used numerous times
540545
Sc3 = Sc**(1./3.)
@@ -889,7 +894,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
889894

890895
if (mpirank==mpiroot) write (*,*)'creating microphysics lookup tables ... '
891896
if (mpirank==mpiroot) write (*,'(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
892-
' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
897+
' using: mu_c_o=',mu_c_o,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
893898

894899
!> - Call table_ccnact() to read a static file containing CCN activation of aerosols. The
895900
!! data were created from a parcel model by Feingold & Heymsfield with
@@ -982,7 +987,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
982987
nwfa, nifa, nwfa2d, nifa2d, &
983988
tt, th, pii, &
984989
p, w, dz, dt_in, dt_inner, &
985-
sedi_semi, decfl, &
990+
sedi_semi, decfl, lsm, &
986991
RAINNC, RAINNCV, &
987992
SNOWNC, SNOWNCV, &
988993
ICENC, ICENCV, &
@@ -1037,6 +1042,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
10371042
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
10381043
nc, nwfa, nifa
10391044
REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d
1045+
INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN):: lsm
10401046
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
10411047
re_cloud, re_ice, re_snow
10421048
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: pfils, pflls
@@ -1117,6 +1123,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
11171123
REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
11181124
REAL:: dt, pptrain, pptsnow, pptgraul, pptice
11191125
REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
1126+
INTEGER:: lsml
11201127
REAL:: rand1, rand2, rand3, rand_pert_max
11211128
INTEGER:: i, j, k, m
11221129
INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
@@ -1419,8 +1426,14 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
14191426
nifa1d(k) = nifa(i,k,j)
14201427
enddo
14211428
else
1429+
lsml = lsm(i,j)
14221430
do k = kts, kte
1423-
nc1d(k) = Nt_c/rho(k)
1431+
!nc1d(k) = Nt_c/rho(k)
1432+
if(lsml == 0) then
1433+
nc1d(k) = Nt_c_o/rho(k)
1434+
else
1435+
nc1d(k) = Nt_c_l/rho(k)
1436+
endif
14241437
nwfa1d(k) = 11.1E6
14251438
nifa1d(k) = naIN1*0.01
14261439
enddo
@@ -1429,7 +1442,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
14291442
!> - Call mp_thompson()
14301443
call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
14311444
nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, &
1432-
pptrain, pptsnow, pptgraul, pptice, &
1445+
lsml, pptrain, pptsnow, pptgraul, pptice, &
14331446
#if ( WRF_CHEM == 1 )
14341447
rainprod1d, evapprod1d, &
14351448
#endif
@@ -1698,7 +1711,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
16981711
enddo
16991712
!> - Call calc_effectrad()
17001713
call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1701-
re_qc1d, re_qi1d, re_qs1d, kts, kte)
1714+
re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
17021715
do k = kts, kte
17031716
re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max))
17041717
re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max))
@@ -1841,7 +1854,7 @@ END SUBROUTINE thompson_finalize
18411854
!> @{
18421855
subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
18431856
nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, &
1844-
pptrain, pptsnow, pptgraul, pptice, &
1857+
lsml, pptrain, pptsnow, pptgraul, pptice, &
18451858
#if ( WRF_CHEM == 1 )
18461859
rainprod, evapprod, &
18471860
#endif
@@ -1879,6 +1892,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
18791892
REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq
18801893
REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
18811894
REAL, INTENT(IN):: dt
1895+
INTEGER, INTENT(IN):: lsml
18821896
REAL, INTENT(IN):: rand1, rand2, rand3
18831897
! Extended diagnostics, most arrays only allocated if ext_diag is true
18841898
LOGICAL, INTENT(IN) :: ext_diag
@@ -1982,6 +1996,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
19821996
REAL:: Ef_ra, Ef_sa, Ef_ga
19831997
REAL:: dtsave, odts, odt, odzq, hgt_agl, SR
19841998
REAL:: xslw1, ygra1, zans1, eva_factor
1999+
REAL:: av_i
19852000
INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
19862001
INTEGER, DIMENSION(5):: ksed1
19872002
INTEGER:: nir, nis, nig, nii, nic, niin
@@ -2006,6 +2021,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
20062021
odt = 1./dt
20072022
odts = 1./dtsave
20082023
iexfrq = 1
2024+
av_i = av_s * D0s ** (bv_s - bv_i)
20092025

20102026
!+---+-----------------------------------------------------------------+
20112027
!> - Initialize Source/sink terms. First 2 chars: "pr" represents source/sink of
@@ -2210,7 +2226,14 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
22102226
endif
22112227
nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
22122228
/ am_r*lamc**bm_r)
2213-
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
2229+
!if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
2230+
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
2231+
if (lsml == 0) then
2232+
nc(k) = Nt_c_o
2233+
else
2234+
nc(k) = Nt_c_l
2235+
endif
2236+
endif
22142237
else
22152238
qc1d(k) = 0.0
22162239
nc1d(k) = 0.0
@@ -2234,7 +2257,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
22342257
if (xDi.lt. 5.E-6) then
22352258
lami = cie(2)/5.E-6
22362259
ni(k) = MIN(4999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2237-
elseif (xDi.gt. 300.E-6) then
2260+
elseif (xDi.gt. D0s + 100.E-6) then
22382261
lami = cie(2)/300.E-6
22392262
ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
22402263
endif
@@ -2919,13 +2942,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
29192942

29202943
!> - Deposition nucleation of dust/mineral from DeMott et al (2010)
29212944
!! we may need to relax the temperature and ssati constraints.
2922-
if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps &
2945+
if ( (ssati(k).ge. 0.15) .or. (ssatw(k).gt. eps &
29232946
.and. temp(k).lt.253.15) ) then
29242947
if (dustyIce .AND. (is_aerosol_aware .or. merra2_aerosol_aware)) then
29252948
xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
29262949
xnc = xnc*(1.0 + 50.*rand3)
29272950
else
2928-
xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
2951+
xnc = MIN(1000.E3, TNO*EXP(ATO*(T_0-temp(k))))
29292952
endif
29302953
xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
29312954
pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
@@ -3273,7 +3296,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
32733296
lami = cie(2)/5.E-6
32743297
xni = MIN(4999.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
32753298
niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3276-
elseif (xDi.gt. 300.E-6) then
3299+
elseif (xDi.gt. D0s + 100.E-6) then
32773300
lami = cie(2)/300.E-6
32783301
xni = cig(1)*oig2*xri/am_i*lami**bm_i
32793302
niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
@@ -3389,7 +3412,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
33893412
if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
33903413
rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
33913414
nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
3392-
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
3415+
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
3416+
if(lsml == 0) then
3417+
nc(k) = Nt_c_o
3418+
else
3419+
nc(k) = Nt_c_l
3420+
endif
3421+
endif
33933422
L_qc(k) = .true.
33943423
else
33953424
rc(k) = R1
@@ -3560,7 +3589,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
35603589
if (is_aerosol_aware .or. merra2_aerosol_aware) then
35613590
xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k)))
35623591
else
3563-
xnc = Nt_c
3592+
if(lsml == 0) then
3593+
xnc = Nt_c_o
3594+
else
3595+
xnc = Nt_c_l
3596+
endif
35643597
endif
35653598
pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
35663599

@@ -3630,7 +3663,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
36303663
rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
36313664
if (rc(k).eq.R1) L_qc(k) = .false.
36323665
nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
3633-
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
3666+
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
3667+
if(lsml == 0) then
3668+
nc(k) = Nt_c_o
3669+
else
3670+
nc(k) = Nt_c_l
3671+
endif
3672+
endif
36343673
qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
36353674
temp(k) = t1d(k) + DT*tten(k)
36363675
rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
@@ -4235,7 +4274,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
42354274
xDi = (bm_i + mu_i + 1.) * ilami
42364275
if (xDi.lt. 5.E-6) then
42374276
lami = cie(2)/5.E-6
4238-
elseif (xDi.gt. 300.E-6) then
4277+
elseif (xDi.gt. D0s + 100.E-6) then
42394278
lami = cie(2)/300.E-6
42404279
endif
42414280
ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
@@ -5749,7 +5788,7 @@ END FUNCTION delta_p
57495788
!! distribution, not the second part, which is the larger sizes.
57505789

57515790
subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
5752-
& re_qc1d, re_qi1d, re_qs1d, kts, kte)
5791+
& re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
57535792

57545793
IMPLICIT NONE
57555794

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

@@ -5781,7 +5821,13 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
57815821
rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
57825822
rc(k) = MAX(R1, qc1d(k)*rho(k))
57835823
nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
5784-
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c
5824+
if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) then
5825+
if( lsml == 0) then
5826+
nc(k) = Nt_c_o
5827+
else
5828+
nc(k) = Nt_c_l
5829+
endif
5830+
endif
57855831
if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true.
57865832
ri(k) = MAX(R1, qi1d(k)*rho(k))
57875833
ni(k) = MAX(R2, ni1d(k)*rho(k))

0 commit comments

Comments
 (0)