ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions/Subroutines | Variables
mp_wsm6 Module Reference

Functions/Subroutines

subroutine, public mp_wsm6_init (den0, denr, dens, cl, cpv, hail_opt, errmsg, errflg)
 
subroutine, public mp_wsm6_finalize (errmsg, errflg)
 
subroutine, public mp_wsm6_run (t, q, qc, qi, qr, qs, qg, den, p, delz, delt, g, cpd, cpv, rd, rv, t0c, ep1, ep2, qmin, xls, xlv0, xlf0, den0, denr, cliq, cice, psat, rain, rainncv, sr, snow, snowncv, graupel, graupelncv, rainprod2d, evapprod2d, its, ite, kts, kte, errmsg, errflg)
 
real(kind=kind_phys) function rgmma (x)
 
real(kind=kind_phys) function fpvs (t, ice, rd, rv, cvap, cliq, cice, hvap, hsub, psat, t0c)
 
subroutine slope_wsm6 (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine slope_rain (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine slope_snow (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine slope_graup (qrs, den, denfac, t, rslope, rslopeb, rslope2, rslope3, vt, its, ite, kts, kte)
 
subroutine nislfv_rain_plm (im, km, denl, denfacl, tkl, dzl, wwl, rql, precip, dt, id, iter)
 
subroutine nislfv_rain_plm6 (im, km, denl, denfacl, tkl, dzl, wwl, rql, rql2, precip1, precip2, dt, id, iter)
 
subroutine, public refl10cm_wsm6 (qv1d, qr1d, qs1d, qg1d, t1d, p1d, dBZ, kts, kte)
 

Variables

integer, parameter kind_phys = c_double
 
real(kind=kind_phys), parameter, private dtcldcr = 120.
 
real(kind=kind_phys), parameter, private n0r = 8.e6
 
real(kind=kind_phys), parameter, private avtr = 841.9
 
real(kind=kind_phys), parameter, private bvtr = 0.8
 
real(kind=kind_phys), parameter, private r0 = .8e-5
 
real(kind=kind_phys), parameter, private peaut = .55
 
real(kind=kind_phys), parameter, private xncr = 3.e8
 
real(kind=kind_phys), parameter, private xmyu = 1.718e-5
 
real(kind=kind_phys), parameter, private avts = 11.72
 
real(kind=kind_phys), parameter, private bvts = .41
 
real(kind=kind_phys), parameter, private lamdarmax = 8.e4
 
real(kind=kind_phys), parameter, private lamdasmax = 1.e5
 
real(kind=kind_phys), parameter, private dicon = 11.9
 
real(kind=kind_phys), parameter, private dimax = 500.e-6
 
real(kind=kind_phys), parameter, private pfrz1 = 100.
 
real(kind=kind_phys), parameter, private pfrz2 = 0.66
 
real(kind=kind_phys), parameter, private qcrmin = 1.e-9
 
real(kind=kind_phys), parameter, private eacrc = 1.0
 
real(kind=kind_phys), parameter, private dens = 100.0
 
real(kind=kind_phys), parameter, private qs0 = 6.e-4
 
real(kind=kind_phys), parameter, public n0smax = 1.e11
 
real(kind=kind_phys), parameter, public n0s = 2.e6
 
real(kind=kind_phys), parameter, public alpha = .12
 
real(kind=kind_phys), save qc0
 
real(kind=kind_phys), save qck1
 
real(kind=kind_phys), save bvtr1
 
real(kind=kind_phys), save bvtr2
 
real(kind=kind_phys), save bvtr3
 
real(kind=kind_phys), save bvtr4
 
real(kind=kind_phys), save g1pbr
 
real(kind=kind_phys), save g3pbr
 
real(kind=kind_phys), save g4pbr
 
real(kind=kind_phys), save g5pbro2
 
real(kind=kind_phys), save pvtr
 
real(kind=kind_phys), save eacrr
 
real(kind=kind_phys), save pacrr
 
real(kind=kind_phys), save bvtr6
 
real(kind=kind_phys), save g6pbr
 
real(kind=kind_phys), save precr1
 
real(kind=kind_phys), save precr2
 
real(kind=kind_phys), save roqimax
 
real(kind=kind_phys), save bvts1
 
real(kind=kind_phys), save bvts2
 
real(kind=kind_phys), save bvts3
 
real(kind=kind_phys), save bvts4
 
real(kind=kind_phys), save g1pbs
 
real(kind=kind_phys), save g3pbs
 
real(kind=kind_phys), save g4pbs
 
real(kind=kind_phys), save n0g
 
real(kind=kind_phys), save avtg
 
real(kind=kind_phys), save bvtg
 
real(kind=kind_phys), save deng
 
real(kind=kind_phys), save lamdagmax
 
real(kind=kind_phys), save g5pbso2
 
real(kind=kind_phys), save pvts
 
real(kind=kind_phys), save pacrs
 
real(kind=kind_phys), save precs1
 
real(kind=kind_phys), save precs2
 
real(kind=kind_phys), save pidn0r
 
real(kind=kind_phys), save xlv1
 
real(kind=kind_phys), save pacrc
 
real(kind=kind_phys), save pi
 
real(kind=kind_phys), save bvtg1
 
real(kind=kind_phys), save bvtg2
 
real(kind=kind_phys), save bvtg3
 
real(kind=kind_phys), save bvtg4
 
real(kind=kind_phys), save g1pbg
 
real(kind=kind_phys), save g3pbg
 
real(kind=kind_phys), save g4pbg
 
real(kind=kind_phys), save g5pbgo2
 
real(kind=kind_phys), save pvtg
 
real(kind=kind_phys), save pacrg
 
real(kind=kind_phys), save precg1
 
real(kind=kind_phys), save precg2
 
real(kind=kind_phys), save pidn0g
 
real(kind=kind_phys), save rslopermax
 
real(kind=kind_phys), save rslopesmax
 
real(kind=kind_phys), save rslopegmax
 
real(kind=kind_phys), save rsloperbmax
 
real(kind=kind_phys), save rslopesbmax
 
real(kind=kind_phys), save rslopegbmax
 
real(kind=kind_phys), save rsloper2max
 
real(kind=kind_phys), save rslopes2max
 
real(kind=kind_phys), save rslopeg2max
 
real(kind=kind_phys), save rsloper3max
 
real(kind=kind_phys), save rslopes3max
 
real(kind=kind_phys), save rslopeg3max
 
real(kind=kind_phys), save, public pidn0s
 
real(kind=kind_phys), save, public pidnc
 

Function/Subroutine Documentation

◆ fpvs()

real(kind=kind_phys) function mp_wsm6::fpvs ( real(kind=kind_phys), intent(in)  t,
integer, intent(in)  ice,
real(kind=kind_phys), intent(in)  rd,
real(kind=kind_phys), intent(in)  rv,
real(kind=kind_phys), intent(in)  cvap,
real(kind=kind_phys), intent(in)  cliq,
real(kind=kind_phys), intent(in)  cice,
real(kind=kind_phys), intent(in)  hvap,
real(kind=kind_phys), intent(in)  hsub,
real(kind=kind_phys), intent(in)  psat,
real(kind=kind_phys), intent(in)  t0c 
)
private
1499 !=================================================================================================================
1500 
1501  integer,intent(in):: ice
1502  real(kind=kind_phys),intent(in):: cice,cliq,cvap,hsub,hvap,psat,rd,rv,t0c
1503  real(kind=kind_phys),intent(in):: t
1504 
1505  real(kind=kind_phys):: tr,ttp,dldt,dldti,xa,xb,xai,xbi
1506 
1507 !-----------------------------------------------------------------------------------------------------------------
1508 
1509  ttp=t0c+0.01
1510  dldt=cvap-cliq
1511  xa=-dldt/rv
1512  xb=xa+hvap/(rv*ttp)
1513  dldti=cvap-cice
1514  xai=-dldti/rv
1515  xbi=xai+hsub/(rv*ttp)
1516  tr=ttp/t
1517  if(t.lt.ttp.and.ice.eq.1) then
1518  fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1519  else
1520  fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1521  endif
1522 

◆ mp_wsm6_finalize()

subroutine, public mp_wsm6::mp_wsm6_finalize ( character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

arg_table_mp_wsm6_finalize

\html

197 !=================================================================================================================
198 
199 !--- output arguments:
200  character(len=*),intent(out):: errmsg
201  integer,intent(out):: errflg
202 
203 !-----------------------------------------------------------------------------------------------------------------
204 
205  errmsg = 'mp_wsm6_finalize OK'
206  errflg = 0
207 

◆ mp_wsm6_init()

subroutine, public mp_wsm6::mp_wsm6_init ( real(kind=kind_phys), intent(in)  den0,
real(kind=kind_phys), intent(in)  denr,
real(kind=kind_phys), intent(in)  dens,
real(kind=kind_phys), intent(in)  cl,
real(kind=kind_phys), intent(in)  cpv,
integer, intent(in)  hail_opt,
character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

arg_table_mp_wsm6_init

\html

75 !=================================================================================================================
76 
77 !input arguments:
78  integer,intent(in):: hail_opt ! RAS
79  real(kind=kind_phys),intent(in):: den0,denr,dens,cl,cpv
80 
81 !output arguments:
82  character(len=*),intent(out):: errmsg
83  integer,intent(out):: errflg
84 
85 !-----------------------------------------------------------------------------------------------------------------
86 
87 ! RAS13.1 define graupel parameters as graupel-like or hail-like,
88 ! depending on namelist option
89  if(hail_opt .eq. 1) then !Hail!
90  n0g = 4.e4
91  deng = 700.
92  avtg = 285.0
93  bvtg = 0.8
94  lamdagmax = 2.e4
95  else !Graupel!
96  n0g = 4.e6
97  deng = 500
98  avtg = 330.0
99  bvtg = 0.8
100  lamdagmax = 6.e4
101  endif
102 !
103  pi = 4.*atan(1.)
104  xlv1 = cl-cpv
105 !
106  qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
107  qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
108  pidnc = pi*denr/6. ! syb
109 !
110  bvtr1 = 1.+bvtr
111  bvtr2 = 2.5+.5*bvtr
112  bvtr3 = 3.+bvtr
113  bvtr4 = 4.+bvtr
114  bvtr6 = 6.+bvtr
115  g1pbr = rgmma(bvtr1)
116  g3pbr = rgmma(bvtr3)
117  g4pbr = rgmma(bvtr4) ! 17.837825
118  g6pbr = rgmma(bvtr6)
119  g5pbro2 = rgmma(bvtr2) ! 1.8273
120  pvtr = avtr*g4pbr/6.
121  eacrr = 1.0
122  pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
123  precr1 = 2.*pi*n0r*.78
124  precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
125  roqimax = 2.08e22*dimax**8
126 !
127  bvts1 = 1.+bvts
128  bvts2 = 2.5+.5*bvts
129  bvts3 = 3.+bvts
130  bvts4 = 4.+bvts
131  g1pbs = rgmma(bvts1) !.8875
132  g3pbs = rgmma(bvts3)
133  g4pbs = rgmma(bvts4) ! 12.0786
134  g5pbso2 = rgmma(bvts2)
135  pvts = avts*g4pbs/6.
136  pacrs = pi*n0s*avts*g3pbs*.25
137  precs1 = 4.*n0s*.65
138  precs2 = 4.*n0s*.44*avts**.5*g5pbso2
139  pidn0r = pi*denr*n0r
140  pidn0s = pi*dens*n0s
141 !
142  pacrc = pi*n0s*avts*g3pbs*.25*eacrc
143 !
144  bvtg1 = 1.+bvtg
145  bvtg2 = 2.5+.5*bvtg
146  bvtg3 = 3.+bvtg
147  bvtg4 = 4.+bvtg
148  g1pbg = rgmma(bvtg1)
149  g3pbg = rgmma(bvtg3)
150  g4pbg = rgmma(bvtg4)
151  pacrg = pi*n0g*avtg*g3pbg*.25
152  g5pbgo2 = rgmma(bvtg2)
153  pvtg = avtg*g4pbg/6.
154  precg1 = 2.*pi*n0g*.78
155  precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
156  pidn0g = pi*deng*n0g
157 !
158  rslopermax = 1./lamdarmax
159  rslopesmax = 1./lamdasmax
160  rslopegmax = 1./lamdagmax
161  rsloperbmax = rslopermax ** bvtr
162  rslopesbmax = rslopesmax ** bvts
163  rslopegbmax = rslopegmax ** bvtg
164  rsloper2max = rslopermax * rslopermax
165  rslopes2max = rslopesmax * rslopesmax
166  rslopeg2max = rslopegmax * rslopegmax
167  rsloper3max = rsloper2max * rslopermax
168  rslopes3max = rslopes2max * rslopesmax
169  rslopeg3max = rslopeg2max * rslopegmax
170 
171 !+---+-----------------------------------------------------------------+
172 !.. Set these variables needed for computing radar reflectivity. These
173 !.. get used within radar_init to create other variables used in the
174 !.. radar module.
175  xam_r = pi*denr/6.
176  xbm_r = 3.
177  xmu_r = 0.
178  xam_s = pi*dens/6.
179  xbm_s = 3.
180  xmu_s = 0.
181  xam_g = pi*deng/6.
182  xbm_g = 3.
183  xmu_g = 0.
184 
185  call radar_init
186 
187  errmsg = 'mp_wsm6_init OK'
188  errflg = 0
189 

Referenced by mp_wsm6_isohelper::mp_wsm6_init_c().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ mp_wsm6_run()

subroutine, public mp_wsm6::mp_wsm6_run ( real(kind=kind_phys), dimension(its:,:), intent(inout)  t,
real(kind=kind_phys), dimension(its:,:), intent(inout)  q,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qc,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qi,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qr,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qs,
real(kind=kind_phys), dimension(its:,:), intent(inout)  qg,
real(kind=kind_phys), dimension(its:,:), intent(in)  den,
real(kind=kind_phys), dimension(its:,:), intent(in)  p,
real(kind=kind_phys), dimension(its:,:), intent(in)  delz,
real(kind=kind_phys), intent(in)  delt,
real(kind=kind_phys), intent(in)  g,
real(kind=kind_phys), intent(in)  cpd,
real(kind=kind_phys), intent(in)  cpv,
real(kind=kind_phys), intent(in)  rd,
real(kind=kind_phys), intent(in)  rv,
real(kind=kind_phys), intent(in)  t0c,
real(kind=kind_phys), intent(in)  ep1,
real(kind=kind_phys), intent(in)  ep2,
real(kind=kind_phys), intent(in)  qmin,
real(kind=kind_phys), intent(in)  xls,
real(kind=kind_phys), intent(in)  xlv0,
real(kind=kind_phys), intent(in)  xlf0,
real(kind=kind_phys), intent(in)  den0,
real(kind=kind_phys), intent(in)  denr,
real(kind=kind_phys), intent(in)  cliq,
real(kind=kind_phys), intent(in)  cice,
real(kind=kind_phys), intent(in)  psat,
real(kind=kind_phys), dimension(its:), intent(inout)  rain,
real(kind=kind_phys), dimension(its:), intent(inout)  rainncv,
real(kind=kind_phys), dimension(its:), intent(inout)  sr,
real(kind=kind_phys), dimension(its:), intent(inout), optional  snow,
real(kind=kind_phys), dimension(its:), intent(inout), optional  snowncv,
real(kind=kind_phys), dimension(its:), intent(inout), optional  graupel,
real(kind=kind_phys), dimension(its:), intent(inout), optional  graupelncv,
real(kind=kind_phys), dimension(its:,:), intent(inout), optional  rainprod2d,
real(kind=kind_phys), dimension(its:,:), intent(inout), optional  evapprod2d,
integer, intent(in)  its,
integer, intent(in)  ite,
integer, intent(in)  kts,
integer, intent(in)  kte,
character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

arg_table_mp_wsm6_run

\html

221 !=================================================================================================================!
222 ! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the
223 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
224 ! number concentration is a function of temperature, and seperate assumption
225 ! is developed, in which ice crystal number concentration is a function
226 ! of ice amount. A theoretical background of the ice-microphysics and related
227 ! processes in the WSMMPs are described in Hong et al. (2004).
228 ! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
229 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
230 !
231 ! WSM6 cloud scheme
232 !
233 ! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
234 ! Summer 2003
235 !
236 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
237 ! Summer 2004
238 !
239 ! further modifications :
240 ! semi-lagrangian sedimentation (JH,2010),hong, aug 2009
241 ! ==> higher accuracy and efficient at lower resolutions
242 ! reflectivity computation from greg thompson, lim, jun 2011
243 ! ==> only diagnostic, but with removal of too large drops
244 ! add hail option from afwa, aug 2014
245 ! ==> switch graupel or hail by changing no, den, fall vel.
246 ! effective radius of hydrometeors, bae from kiaps, jan 2015
247 ! ==> consistency in solar insolation of rrtmg radiation
248 ! bug fix in melting terms, bae from kiaps, nov 2015
249 ! ==> density of air is divided, which has not been
250 !
251 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
252 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
253 ! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
254 ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
255 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
256 ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
257 ! Juang and Hong (JH, 2010) Mon. Wea. Rev.
258 !
259 
260 !input arguments:
261  integer,intent(in):: its,ite,kts,kte
262 
263  real(kind=kind_phys),intent(in),dimension(its:,:):: &
264  den, &
265  p, &
266  delz
267  real(kind=kind_phys),intent(in):: &
268  delt, &
269  g, &
270  cpd, &
271  cpv, &
272  t0c, &
273  den0, &
274  rd, &
275  rv, &
276  ep1, &
277  ep2, &
278  qmin, &
279  xls, &
280  xlv0, &
281  xlf0, &
282  cliq, &
283  cice, &
284  psat, &
285  denr
286 
287 !inout arguments:
288  real(kind=kind_phys),intent(inout),dimension(its:,:):: &
289  t
290  real(kind=kind_phys),intent(inout),dimension(its:,:):: &
291  q, &
292  qc, &
293  qi, &
294  qr, &
295  qs, &
296  qg
297  real(kind=kind_phys),intent(inout),dimension(its:):: &
298  rain, &
299  rainncv, &
300  sr
301 
302  real(kind=kind_phys),intent(inout),dimension(its:),optional:: &
303  snow, &
304  snowncv
305 
306  real(kind=kind_phys),intent(inout),dimension(its:),optional:: &
307  graupel, &
308  graupelncv
309 
310  real(kind=kind_phys),intent(inout),dimension(its:,:),optional:: &
311  rainprod2d, &
312  evapprod2d
313 
314 !output arguments:
315  character(len=*),intent(out):: errmsg
316  integer,intent(out):: errflg
317 
318 !local variables and arrays:
319  real(kind=kind_phys),dimension(its:ite,kts:kte,3):: &
320  rh, &
321  qsat, &
322  rslope, &
323  rslope2, &
324  rslope3, &
325  rslopeb, &
326  qrs_tmp, &
327  falk, &
328  fall, &
329  work1
330  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
331  fallc, &
332  falkc, &
333  work1c, &
334  work2c, &
335  workr, &
336  worka
337  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
338  den_tmp, &
339  delz_tmp
340  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
341  pigen, &
342  pidep, &
343  pcond, &
344  prevp, &
345  psevp, &
346  pgevp, &
347  psdep, &
348  pgdep, &
349  praut, &
350  psaut, &
351  pgaut, &
352  piacr, &
353  pracw, &
354  praci, &
355  pracs, &
356  psacw, &
357  psaci, &
358  psacr, &
359  pgacw, &
360  pgaci, &
361  pgacr, &
362  pgacs, &
363  paacw, &
364  psmlt, &
365  pgmlt, &
366  pseml, &
367  pgeml
368  real(kind=kind_phys),dimension(its:ite,kts:kte):: &
369  qsum, &
370  xl, &
371  cpm, &
372  work2, &
373  denfac, &
374  xni, &
375  denqrs1, &
376  denqrs2, &
377  denqrs3, &
378  denqci, &
379  n0sfac
380  real(kind=kind_phys),dimension(its:ite):: &
381  delqrs1, &
382  delqrs2, &
383  delqrs3, &
384  delqi
385  real(kind=kind_phys),dimension(its:ite):: &
386  tstepsnow, &
387  tstepgraup
388  integer,dimension(its:ite):: &
389  mstep, &
390  numdt
391  logical,dimension(its:ite):: flgcld
392  real(kind=kind_phys):: &
393  cpmcal, xlcal, diffus, &
394  viscos, xka, venfac, conden, diffac, &
395  x, y, z, a, b, c, d, e, &
396  qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
397  coeres, supsat, dtcld, xmi, eacrs, satdt, &
398  qimax, diameter, xni0, roqi0, &
399  fallsum, fallsum_qsi, fallsum_qg, &
400  vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
401  xlwork2, factor, source, value, &
402  xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
403  real(kind=kind_phys):: vt2ave
404  real(kind=kind_phys):: holdc, holdci
405  integer:: i, j, k, mstepmax, &
406  iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
407 
408 !Temporaries used for inlining fpvs function
409  real(kind=kind_phys):: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
410 
411 ! variables for optimization
412  real(kind=kind_phys),dimension(its:ite):: dvec1,tvec1
413  real(kind=kind_phys):: temp
414 
415 !-----------------------------------------------------------------------------------------------------------------
416 
417 ! compute internal functions
418 !
419  cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
420  xlcal(x) = xlv0-xlv1*(x-t0c)
421 !----------------------------------------------------------------
422 ! diffus: diffusion coefficient of the water vapor
423 ! viscos: kinematic viscosity(m2s-1)
424 ! Optimizatin : A**B => exp(log(A)*(B))
425 !
426  diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
427  viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
428  xka(x,y) = 1.414e3*viscos(x,y)*y
429  diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
430  venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
431  /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
432  conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
433 !
434 !
435  idim = ite-its+1
436  kdim = kte-kts+1
437 !
438 !----------------------------------------------------------------
439 ! paddint 0 for negative values generated by dynamics
440 !
441  do k = kts, kte
442  do i = its, ite
443  qc(i,k) = max(qc(i,k),0.0)
444  qr(i,k) = max(qr(i,k),0.0)
445  qi(i,k) = max(qi(i,k),0.0)
446  qs(i,k) = max(qs(i,k),0.0)
447  qg(i,k) = max(qg(i,k),0.0)
448  enddo
449  enddo
450 !
451 !----------------------------------------------------------------
452 ! latent heat for phase changes and heat capacity. neglect the
453 ! changes during microphysical process calculation emanuel(1994)
454 !
455  do k = kts, kte
456  do i = its, ite
457  cpm(i,k) = cpmcal(q(i,k))
458  xl(i,k) = xlcal(t(i,k))
459  enddo
460  enddo
461  do k = kts, kte
462  do i = its, ite
463  delz_tmp(i,k) = delz(i,k)
464  den_tmp(i,k) = den(i,k)
465  enddo
466  enddo
467 !
468 !----------------------------------------------------------------
469 ! initialize the surface rain, snow, graupel
470 !
471  do i = its, ite
472  rainncv(i) = 0.
473  if(present(snowncv) .and. present(snow)) snowncv(i) = 0.
474  if(present(graupelncv) .and. present(graupel)) graupelncv(i) = 0.
475  sr(i) = 0.
476 ! new local array to catch step snow and graupel
477  tstepsnow(i) = 0.
478  tstepgraup(i) = 0.
479  enddo
480 !
481 !----------------------------------------------------------------
482 ! compute the minor time steps.
483 !
484  loops = max(nint(delt/dtcldcr),1)
485  dtcld = delt/loops
486  if(delt.le.dtcldcr) dtcld = delt
487 !
488  do loop = 1,loops
489 !
490 !----------------------------------------------------------------
491 ! initialize the large scale variables
492 !
493  do i = its, ite
494  mstep(i) = 1
495  flgcld(i) = .true.
496  enddo
497 !
498 ! do k = kts, kte
499 ! do i = its, ite
500 ! denfac(i,k) = sqrt(den0/den(i,k))
501 ! enddo
502 ! enddo
503  do k = kts, kte
504  do i = its,ite
505  dvec1(i) = den(i,k)
506  enddo
507  call vrec(tvec1,dvec1,ite-its+1)
508  do i = its, ite
509  tvec1(i) = tvec1(i)*den0
510  enddo
511  call vsqrt(dvec1,tvec1,ite-its+1)
512  do i = its,ite
513  denfac(i,k) = dvec1(i)
514  enddo
515  enddo
516 !
517 ! Inline expansion for fpvs
518 ! qsat(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
519 ! qsat(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
520  hsub = xls
521  hvap = xlv0
522  cvap = cpv
523  ttp=t0c+0.01
524  dldt=cvap-cliq
525  xa=-dldt/rv
526  xb=xa+hvap/(rv*ttp)
527  dldti=cvap-cice
528  xai=-dldti/rv
529  xbi=xai+hsub/(rv*ttp)
530  do k = kts, kte
531  do i = its, ite
532  tr=ttp/t(i,k)
533  qsat(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
534  qsat(i,k,1) = min(qsat(i,k,1),0.99*p(i,k))
535  qsat(i,k,1) = ep2 * qsat(i,k,1) / (p(i,k) - qsat(i,k,1))
536  qsat(i,k,1) = max(qsat(i,k,1),qmin)
537  rh(i,k,1) = max(q(i,k) / qsat(i,k,1),qmin)
538  tr=ttp/t(i,k)
539  if(t(i,k).lt.ttp) then
540  qsat(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
541  else
542  qsat(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
543  endif
544  qsat(i,k,2) = min(qsat(i,k,2),0.99*p(i,k))
545  qsat(i,k,2) = ep2 * qsat(i,k,2) / (p(i,k) - qsat(i,k,2))
546  qsat(i,k,2) = max(qsat(i,k,2),qmin)
547  rh(i,k,2) = max(q(i,k) / qsat(i,k,2),qmin)
548  enddo
549  enddo
550 !
551 !----------------------------------------------------------------
552 ! initialize the variables for microphysical physics
553 !
554 !
555  do k = kts, kte
556  do i = its, ite
557  prevp(i,k) = 0.
558  psdep(i,k) = 0.
559  pgdep(i,k) = 0.
560  praut(i,k) = 0.
561  psaut(i,k) = 0.
562  pgaut(i,k) = 0.
563  pracw(i,k) = 0.
564  praci(i,k) = 0.
565  piacr(i,k) = 0.
566  psaci(i,k) = 0.
567  psacw(i,k) = 0.
568  pracs(i,k) = 0.
569  psacr(i,k) = 0.
570  pgacw(i,k) = 0.
571  paacw(i,k) = 0.
572  pgaci(i,k) = 0.
573  pgacr(i,k) = 0.
574  pgacs(i,k) = 0.
575  pigen(i,k) = 0.
576  pidep(i,k) = 0.
577  pcond(i,k) = 0.
578  psmlt(i,k) = 0.
579  pgmlt(i,k) = 0.
580  pseml(i,k) = 0.
581  pgeml(i,k) = 0.
582  psevp(i,k) = 0.
583  pgevp(i,k) = 0.
584  falk(i,k,1) = 0.
585  falk(i,k,2) = 0.
586  falk(i,k,3) = 0.
587  fall(i,k,1) = 0.
588  fall(i,k,2) = 0.
589  fall(i,k,3) = 0.
590  fallc(i,k) = 0.
591  falkc(i,k) = 0.
592  xni(i,k) = 1.e3
593  enddo
594  enddo
595 !-------------------------------------------------------------
596 ! Ni: ice crystal number concentraiton [HDC 5c]
597 !-------------------------------------------------------------
598  do k = kts, kte
599  do i = its, ite
600  temp = (den(i,k)*max(qi(i,k),qmin))
601  temp = sqrt(sqrt(temp*temp*temp))
602  xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
603  enddo
604  enddo
605 !
606 !----------------------------------------------------------------
607 ! compute the fallout term:
608 ! first, vertical terminal velosity for minor loops
609 !----------------------------------------------------------------
610  do k = kts, kte
611  do i = its, ite
612  qrs_tmp(i,k,1) = qr(i,k)
613  qrs_tmp(i,k,2) = qs(i,k)
614  qrs_tmp(i,k,3) = qg(i,k)
615  enddo
616  enddo
617  call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
618  work1,its,ite,kts,kte)
619 !
620  do k = kte, kts, -1
621  do i = its, ite
622  workr(i,k) = work1(i,k,1)
623  qsum(i,k) = max( (qs(i,k)+qg(i,k)), 1.e-15)
624  if( qsum(i,k) .gt. 1.e-15 ) then
625  worka(i,k) = (work1(i,k,2)*qs(i,k) + work1(i,k,3)*qg(i,k)) &
626  / qsum(i,k)
627  else
628  worka(i,k) = 0.
629  endif
630  denqrs1(i,k) = den(i,k)*qr(i,k)
631  denqrs2(i,k) = den(i,k)*qs(i,k)
632  denqrs3(i,k) = den(i,k)*qg(i,k)
633  if(qr(i,k).le.0.0) workr(i,k) = 0.0
634  enddo
635  enddo
636  call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
637  delqrs1,dtcld,1,1)
638  call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
639  denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
640  do k = kts, kte
641  do i = its, ite
642  qr(i,k) = max(denqrs1(i,k)/den(i,k),0.)
643  qs(i,k) = max(denqrs2(i,k)/den(i,k),0.)
644  qg(i,k) = max(denqrs3(i,k)/den(i,k),0.)
645  fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
646  fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
647  fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
648  enddo
649  enddo
650  do i = its, ite
651  fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
652  fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
653  fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
654  enddo
655  do k = kts, kte
656  do i = its, ite
657  qrs_tmp(i,k,1) = qr(i,k)
658  qrs_tmp(i,k,2) = qs(i,k)
659  qrs_tmp(i,k,3) = qg(i,k)
660  enddo
661  enddo
662  call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
663  work1,its,ite,kts,kte)
664 !
665  do k = kte, kts, -1
666  do i = its, ite
667  supcol = t0c-t(i,k)
668  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
669  if(t(i,k).gt.t0c) then
670 !---------------------------------------------------------------
671 ! psmlt: melting of snow [HL A33] [RH83 A25]
672 ! (T>T0: S->R)
673 !---------------------------------------------------------------
674  xlf = xlf0
675  work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
676  if(qs(i,k).gt.0.) then
677  coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
678  psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
679  *n0sfac(i,k)*(precs1*rslope2(i,k,2) &
680  +precs2*work2(i,k)*coeres)/den(i,k)
681  psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
682  -qs(i,k)/mstep(i)),0.)
683  qs(i,k) = qs(i,k) + psmlt(i,k)
684  qr(i,k) = qr(i,k) - psmlt(i,k)
685  t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
686  endif
687 !---------------------------------------------------------------
688 ! pgmlt: melting of graupel [HL A23] [LFO 47]
689 ! (T>T0: G->R)
690 !---------------------------------------------------------------
691  if(qg(i,k).gt.0.) then
692  coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
693  pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
694  *(t0c-t(i,k))*(precg1*rslope2(i,k,3) &
695  +precg2*work2(i,k)*coeres)/den(i,k)
696  pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
697  -qg(i,k)/mstep(i)),0.)
698  qg(i,k) = qg(i,k) + pgmlt(i,k)
699  qr(i,k) = qr(i,k) - pgmlt(i,k)
700  t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
701  endif
702  endif
703  enddo
704  enddo
705 !---------------------------------------------------------------
706 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
707 !---------------------------------------------------------------
708  do k = kte, kts, -1
709  do i = its, ite
710  if(qi(i,k).le.0.) then
711  work1c(i,k) = 0.
712  else
713  xmi = den(i,k)*qi(i,k)/xni(i,k)
714  diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
715  work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
716  endif
717  enddo
718  enddo
719 !
720 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
721 !
722  do k = kte, kts, -1
723  do i = its, ite
724  denqci(i,k) = den(i,k)*qi(i,k)
725  enddo
726  enddo
727  call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
728  delqi,dtcld,1,0)
729  do k = kts, kte
730  do i = its, ite
731  qi(i,k) = max(denqci(i,k)/den(i,k),0.)
732  enddo
733  enddo
734  do i = its, ite
735  fallc(i,1) = delqi(i)/delz(i,1)/dtcld
736  enddo
737 !
738 !----------------------------------------------------------------
739 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
740 !
741  do i = its, ite
742  fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
743  fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
744  fallsum_qg = fall(i,kts,3)
745  if(fallsum.gt.0.) then
746  rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
747  rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
748  endif
749  if(fallsum_qsi.gt.0.) then
750  tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
751  + tstepsnow(i)
752  if(present(snowncv) .and. present(snow)) then
753  snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
754  + snowncv(i)
755  snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
756  endif
757  endif
758  if(fallsum_qg.gt.0.) then
759  tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
760  + tstepgraup(i)
761  if(present (graupelncv) .and. present (graupel)) then
762  graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
763  + graupelncv(i)
764  graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
765  endif
766  endif
767  if(present (snowncv)) then
768  if(fallsum.gt.0.)sr(i)=(snowncv(i) + graupelncv(i))/(rainncv(i)+1.e-12)
769  else
770  if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
771  endif
772  enddo
773 !
774 !---------------------------------------------------------------
775 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
776 ! (T>T0: I->C)
777 !---------------------------------------------------------------
778  do k = kts, kte
779  do i = its, ite
780  supcol = t0c-t(i,k)
781  xlf = xls-xl(i,k)
782  if(supcol.lt.0.) xlf = xlf0
783  if(supcol.lt.0.and.qi(i,k).gt.0.) then
784  qc(i,k) = qc(i,k) + qi(i,k)
785  t(i,k) = t(i,k) - xlf/cpm(i,k)*qi(i,k)
786  qi(i,k) = 0.
787  endif
788 !---------------------------------------------------------------
789 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
790 ! (T<-40C: C->I)
791 !---------------------------------------------------------------
792  if(supcol.gt.40..and.qc(i,k).gt.0.) then
793  qi(i,k) = qi(i,k) + qc(i,k)
794  t(i,k) = t(i,k) + xlf/cpm(i,k)*qc(i,k)
795  qc(i,k) = 0.
796  endif
797 !---------------------------------------------------------------
798 ! pihtf: heterogeneous freezing of cloud water [HL A44]
799 ! (T0>T>-40C: C->I)
800 !---------------------------------------------------------------
801  if(supcol.gt.0..and.qc(i,k).gt.qmin) then
802 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
803 ! * den(i,k)/denr/xncr*qc(i,k)**2*dtcld,qc(i,k))
804  supcolt=min(supcol,50.)
805  pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
806  * den(i,k)/denr/xncr*qc(i,k)*qc(i,k)*dtcld,qc(i,k))
807  qi(i,k) = qi(i,k) + pfrzdtc
808  t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
809  qc(i,k) = qc(i,k)-pfrzdtc
810  endif
811 !---------------------------------------------------------------
812 ! pgfrz: freezing of rain water [HL A20] [LFO 45]
813 ! (T<T0, R->G)
814 !---------------------------------------------------------------
815  if(supcol.gt.0..and.qr(i,k).gt.0.) then
816 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
817 ! * (exp(pfrz2*supcol)-1.)*rslope3(i,k,1)**2 &
818 ! * rslope(i,k,1)*dtcld,qr(i,k))
819  temp = rslope3(i,k,1)
820  temp = temp*temp*rslope(i,k,1)
821  supcolt=min(supcol,50.)
822  pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
823  *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
824  qr(i,k))
825  qg(i,k) = qg(i,k) + pfrzdtr
826  t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
827  qr(i,k) = qr(i,k)-pfrzdtr
828  endif
829  enddo
830  enddo
831 !
832 !
833 !----------------------------------------------------------------
834 ! update the slope parameters for microphysics computation
835 !
836  do k = kts, kte
837  do i = its, ite
838  qrs_tmp(i,k,1) = qr(i,k)
839  qrs_tmp(i,k,2) = qs(i,k)
840  qrs_tmp(i,k,3) = qg(i,k)
841  enddo
842  enddo
843  call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
844  work1,its,ite,kts,kte)
845 !------------------------------------------------------------------
846 ! work1: the thermodynamic term in the denominator associated with
847 ! heat conduction and vapor diffusion
848 ! (ry88, y93, h85)
849 ! work2: parameter associated with the ventilation effects(y93)
850 !
851  do k = kts, kte
852  do i = its, ite
853  work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qsat(i,k,1))
854  work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qsat(i,k,2))
855  work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
856  enddo
857  enddo
858 !
859 !===============================================================
860 !
861 ! warm rain processes
862 !
863 ! - follows the processes in RH83 and LFO except for autoconcersion
864 !
865 !===============================================================
866 !
867  do k = kts, kte
868  do i = its, ite
869  supsat = max(q(i,k),qmin)-qsat(i,k,1)
870  satdt = supsat/dtcld
871 !---------------------------------------------------------------
872 ! praut: auto conversion rate from cloud to rain [HDC 16]
873 ! (C->R)
874 !---------------------------------------------------------------
875  if(qc(i,k).gt.qc0) then
876  praut(i,k) = qck1*qc(i,k)**(7./3.)
877  praut(i,k) = min(praut(i,k),qc(i,k)/dtcld)
878  endif
879 !---------------------------------------------------------------
880 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
881 ! (C->R)
882 !---------------------------------------------------------------
883  if(qr(i,k).gt.qcrmin.and.qc(i,k).gt.qmin) then
884  pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
885  * qc(i,k)*denfac(i,k),qc(i,k)/dtcld)
886  endif
887 !---------------------------------------------------------------
888 ! prevp: evaporation/condensation rate of rain [HDC 14]
889 ! (V->R or R->V)
890 !---------------------------------------------------------------
891  if(qr(i,k).gt.0.) then
892  coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
893  prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
894  + precr2*work2(i,k)*coeres)/work1(i,k,1)
895  if(prevp(i,k).lt.0.) then
896  prevp(i,k) = max(prevp(i,k),-qr(i,k)/dtcld)
897  prevp(i,k) = max(prevp(i,k),satdt/2)
898  else
899  prevp(i,k) = min(prevp(i,k),satdt/2)
900  endif
901  endif
902  enddo
903  enddo
904 !
905 !===============================================================
906 !
907 ! cold rain processes
908 !
909 ! - follows the revised ice microphysics processes in HDC
910 ! - the processes same as in RH83 and RH84 and LFO behave
911 ! following ice crystal hapits defined in HDC, inclduing
912 ! intercept parameter for snow (n0s), ice crystal number
913 ! concentration (ni), ice nuclei number concentration
914 ! (n0i), ice diameter (d)
915 !
916 !===============================================================
917 !
918  do k = kts, kte
919  do i = its, ite
920  supcol = t0c-t(i,k)
921  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
922  supsat = max(q(i,k),qmin)-qsat(i,k,2)
923  satdt = supsat/dtcld
924  ifsat = 0
925 !-------------------------------------------------------------
926 ! Ni: ice crystal number concentraiton [HDC 5c]
927 !-------------------------------------------------------------
928 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
929 ! * max(qi(i,k),qmin))**0.75,1.e3),1.e6)
930  temp = (den(i,k)*max(qi(i,k),qmin))
931  temp = sqrt(sqrt(temp*temp*temp))
932  xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
933  eacrs = exp(0.07*(-supcol))
934 !
935  xmi = den(i,k)*qi(i,k)/xni(i,k)
936  diameter = min(dicon * sqrt(xmi),dimax)
937  vt2i = 1.49e4*diameter**1.31
938  vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
939  vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
940  vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
941  qsum(i,k) = max( (qs(i,k)+qg(i,k)), 1.e-15)
942  if(qsum(i,k) .gt. 1.e-15) then
943  vt2ave=(vt2s*qs(i,k)+vt2g*qg(i,k))/(qsum(i,k))
944  else
945  vt2ave=0.
946  endif
947  if(supcol.gt.0.and.qi(i,k).gt.qmin) then
948  if(qr(i,k).gt.qcrmin) then
949 !-------------------------------------------------------------
950 ! praci: accretion of cloud ice by rain [HL A15] [LFO 25]
951 ! (T<T0: I->R)
952 !-------------------------------------------------------------
953  acrfac = 2.*rslope3(i,k,1)+2.*diameter*rslope2(i,k,1) &
954  + diameter**2*rslope(i,k,1)
955  praci(i,k) = pi*qi(i,k)*n0r*abs(vt2r-vt2i)*acrfac/4.
956 ! reduce collection efficiency (suggested by B. Wilt)
957  praci(i,k) = praci(i,k)*min(max(0.0,qr(i,k)/qi(i,k)),1.)**2
958  praci(i,k) = min(praci(i,k),qi(i,k)/dtcld)
959 !-------------------------------------------------------------
960 ! piacr: accretion of rain by cloud ice [HL A19] [LFO 26]
961 ! (T<T0: R->S or R->G)
962 !-------------------------------------------------------------
963  piacr(i,k) = pi**2*avtr*n0r*denr*xni(i,k)*denfac(i,k) &
964  * g6pbr*rslope3(i,k,1)*rslope3(i,k,1) &
965  * rslopeb(i,k,1)/24./den(i,k)
966 ! reduce collection efficiency (suggested by B. Wilt)
967  piacr(i,k) = piacr(i,k)*min(max(0.0,qi(i,k)/qr(i,k)),1.)**2
968  piacr(i,k) = min(piacr(i,k),qr(i,k)/dtcld)
969  endif
970 !-------------------------------------------------------------
971 ! psaci: accretion of cloud ice by snow [HDC 10]
972 ! (T<T0: I->S)
973 !-------------------------------------------------------------
974  if(qs(i,k).gt.qcrmin) then
975  acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
976  + diameter**2*rslope(i,k,2)
977  psaci(i,k) = pi*qi(i,k)*eacrs*n0s*n0sfac(i,k) &
978  * abs(vt2ave-vt2i)*acrfac/4.
979  psaci(i,k) = min(psaci(i,k),qi(i,k)/dtcld)
980  endif
981 !-------------------------------------------------------------
982 ! pgaci: accretion of cloud ice by graupel [HL A17] [LFO 41]
983 ! (T<T0: I->G)
984 !-------------------------------------------------------------
985  if(qg(i,k).gt.qcrmin) then
986  egi = exp(0.07*(-supcol))
987  acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
988  + diameter**2*rslope(i,k,3)
989  pgaci(i,k) = pi*egi*qi(i,k)*n0g*abs(vt2ave-vt2i)*acrfac/4.
990  pgaci(i,k) = min(pgaci(i,k),qi(i,k)/dtcld)
991  endif
992  endif
993 !-------------------------------------------------------------
994 ! psacw: accretion of cloud water by snow [HL A7] [LFO 24]
995 ! (T<T0: C->S, and T>=T0: C->R)
996 !-------------------------------------------------------------
997  if(qs(i,k).gt.qcrmin.and.qc(i,k).gt.qmin) then
998  psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
999 ! reduce collection efficiency (suggested by B. Wilt)
1000  * min(max(0.0,qs(i,k)/qc(i,k)),1.)**2 &
1001  * qc(i,k)*denfac(i,k),qc(i,k)/dtcld)
1002  endif
1003 !-------------------------------------------------------------
1004 ! pgacw: accretion of cloud water by graupel [HL A6] [LFO 40]
1005 ! (T<T0: C->G, and T>=T0: C->R)
1006 !-------------------------------------------------------------
1007  if(qg(i,k).gt.qcrmin.and.qc(i,k).gt.qmin) then
1008  pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3) &
1009 ! reduce collection efficiency (suggested by B. Wilt)
1010  * min(max(0.0,qg(i,k)/qc(i,k)),1.)**2 &
1011  * qc(i,k)*denfac(i,k),qc(i,k)/dtcld)
1012  endif
1013 !-------------------------------------------------------------
1014 ! paacw: accretion of cloud water by averaged snow/graupel
1015 ! (T<T0: C->G or S, and T>=T0: C->R)
1016 !-------------------------------------------------------------
1017  if(qsum(i,k) .gt. 1.e-15) then
1018  paacw(i,k) = (qs(i,k)*psacw(i,k)+qg(i,k)*pgacw(i,k)) &
1019  /(qsum(i,k))
1020  endif
1021 !-------------------------------------------------------------
1022 ! pracs: accretion of snow by rain [HL A11] [LFO 27]
1023 ! (T<T0: S->G)
1024 !-------------------------------------------------------------
1025  if(qs(i,k).gt.qcrmin.and.qr(i,k).gt.qcrmin) then
1026  if(supcol.gt.0) then
1027  acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,1) &
1028  + 2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1) &
1029  + .5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,1)
1030  pracs(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
1031  * (dens/den(i,k))*acrfac
1032 ! reduce collection efficiency (suggested by B. Wilt)
1033  pracs(i,k) = pracs(i,k)*min(max(0.0,qr(i,k)/qs(i,k)),1.)**2
1034  pracs(i,k) = min(pracs(i,k),qs(i,k)/dtcld)
1035  endif
1036 !-------------------------------------------------------------
1037 ! psacr: accretion of rain by snow [HL A10] [LFO 28]
1038 ! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
1039 !-------------------------------------------------------------
1040  acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,2) &
1041  + 2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
1042  +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,2)
1043  psacr(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
1044  * (denr/den(i,k))*acrfac
1045 ! reduce collection efficiency (suggested by B. Wilt)
1046  psacr(i,k) = psacr(i,k)*min(max(0.0,qs(i,k)/qr(i,k)),1.)**2
1047  psacr(i,k) = min(psacr(i,k),qr(i,k)/dtcld)
1048  endif
1049 !-------------------------------------------------------------
1050 ! pgacr: accretion of rain by graupel [HL A12] [LFO 42]
1051 ! (T<T0: R->G) (T>=T0: enhance melting of graupel)
1052 !-------------------------------------------------------------
1053  if(qg(i,k).gt.qcrmin.and.qr(i,k).gt.qcrmin) then
1054  acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,3) &
1055  + 2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
1056  + .5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,3)
1057  pgacr(i,k) = pi**2*n0r*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
1058  * acrfac
1059 ! reduce collection efficiency (suggested by B. Wilt)
1060  pgacr(i,k) = pgacr(i,k)*min(max(0.0,qg(i,k)/qr(i,k)),1.)**2
1061  pgacr(i,k) = min(pgacr(i,k),qr(i,k)/dtcld)
1062  endif
1063 !
1064 !-------------------------------------------------------------
1065 ! pgacs: accretion of snow by graupel [HL A13] [LFO 29]
1066 ! (S->G): This process is eliminated in V3.0 with the
1067 ! new combined snow/graupel fall speeds
1068 !-------------------------------------------------------------
1069  if(qg(i,k).gt.qcrmin.and.qs(i,k).gt.qcrmin) then
1070  pgacs(i,k) = 0.
1071  endif
1072  if(supcol.le.0) then
1073  xlf = xlf0
1074 !-------------------------------------------------------------
1075 ! pseml: enhanced melting of snow by accretion of water [HL A34]
1076 ! (T>=T0: S->R)
1077 !-------------------------------------------------------------
1078  if(qs(i,k).gt.0.) &
1079  pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
1080  / xlf,-qs(i,k)/dtcld),0.)
1081 !-------------------------------------------------------------
1082 ! pgeml: enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1083 ! (T>=T0: G->R)
1084 !-------------------------------------------------------------
1085  if(qg(i,k).gt.0.) &
1086  pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k)) &
1087  / xlf,-qg(i,k)/dtcld),0.)
1088  endif
1089  if(supcol.gt.0) then
1090 !-------------------------------------------------------------
1091 ! pidep: deposition/Sublimation rate of ice [HDC 9]
1092 ! (T<T0: V->I or I->V)
1093 !-------------------------------------------------------------
1094  if(qi(i,k).gt.0.and.ifsat.ne.1) then
1095  pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1096  supice = satdt-prevp(i,k)
1097  if(pidep(i,k).lt.0.) then
1098  pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1099  pidep(i,k) = max(pidep(i,k),-qi(i,k)/dtcld)
1100  else
1101  pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1102  endif
1103  if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1104  endif
1105 !-------------------------------------------------------------
1106 ! psdep: deposition/sublimation rate of snow [HDC 14]
1107 ! (T<T0: V->S or S->V)
1108 !-------------------------------------------------------------
1109  if(qs(i,k).gt.0..and.ifsat.ne.1) then
1110  coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1111  psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1112  + precs2*work2(i,k)*coeres)/work1(i,k,2)
1113  supice = satdt-prevp(i,k)-pidep(i,k)
1114  if(psdep(i,k).lt.0.) then
1115  psdep(i,k) = max(psdep(i,k),-qs(i,k)/dtcld)
1116  psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1117  else
1118  psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1119  endif
1120  if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1121  ifsat = 1
1122  endif
1123 !-------------------------------------------------------------
1124 ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1125 ! (T<T0: V->G or G->V)
1126 !-------------------------------------------------------------
1127  if(qg(i,k).gt.0..and.ifsat.ne.1) then
1128  coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1129  pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
1130  + precg2*work2(i,k)*coeres)/work1(i,k,2)
1131  supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1132  if(pgdep(i,k).lt.0.) then
1133  pgdep(i,k) = max(pgdep(i,k),-qg(i,k)/dtcld)
1134  pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1135  else
1136  pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1137  endif
1138  if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1139  abs(satdt)) ifsat = 1
1140  endif
1141 !-------------------------------------------------------------
1142 ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1143 ! (T<T0: V->I)
1144 !-------------------------------------------------------------
1145  if(supsat.gt.0.and.ifsat.ne.1) then
1146  supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1147  xni0 = 1.e3*exp(0.1*supcol)
1148  roqi0 = 4.92e-11*xni0**1.33
1149  pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qi(i,k),0.))/dtcld)
1150  pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1151  endif
1152 !
1153 !-------------------------------------------------------------
1154 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1155 ! (T<T0: I->S)
1156 !-------------------------------------------------------------
1157  if(qi(i,k).gt.0.) then
1158  qimax = roqimax/den(i,k)
1159  psaut(i,k) = max(0.,(qi(i,k)-qimax)/dtcld)
1160  endif
1161 !
1162 !-------------------------------------------------------------
1163 ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1164 ! (T<T0: S->G)
1165 !-------------------------------------------------------------
1166  if(qs(i,k).gt.0.) then
1167  alpha2 = 1.e-3*exp(0.09*(-supcol))
1168  pgaut(i,k) = min(max(0.,alpha2*(qs(i,k)-qs0)),qs(i,k)/dtcld)
1169  endif
1170  endif
1171 !
1172 !-------------------------------------------------------------
1173 ! psevp: evaporation of melting snow [HL A35] [RH83 A27]
1174 ! (T>=T0: S->V)
1175 !-------------------------------------------------------------
1176  if(supcol.lt.0.) then
1177  if(qs(i,k).gt.0..and.rh(i,k,1).lt.1.) then
1178  coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1179  psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1 &
1180  * rslope2(i,k,2)+precs2*work2(i,k) &
1181  * coeres)/work1(i,k,1)
1182  psevp(i,k) = min(max(psevp(i,k),-qs(i,k)/dtcld),0.)
1183  endif
1184 !-------------------------------------------------------------
1185 ! pgevp: evaporation of melting graupel [HL A25] [RH84 A19]
1186 ! (T>=T0: G->V)
1187 !-------------------------------------------------------------
1188  if(qg(i,k).gt.0..and.rh(i,k,1).lt.1.) then
1189  coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1190  pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
1191  + precg2*work2(i,k)*coeres)/work1(i,k,1)
1192  pgevp(i,k) = min(max(pgevp(i,k),-qg(i,k)/dtcld),0.)
1193  endif
1194  endif
1195  enddo
1196  enddo
1197 !
1198 !
1199 !----------------------------------------------------------------
1200 ! check mass conservation of generation terms and feedback to the
1201 ! large scale
1202 !
1203  do k = kts, kte
1204  do i = its, ite
1205 !
1206  delta2=0.
1207  delta3=0.
1208  if(qr(i,k).lt.1.e-4.and.qs(i,k).lt.1.e-4) delta2=1.
1209  if(qr(i,k).lt.1.e-4) delta3=1.
1210  if(t(i,k).le.t0c) then
1211 !
1212 ! cloud water
1213 !
1214  value = max(qmin,qc(i,k))
1215  source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1216  if (source.gt.value) then
1217  factor = value/source
1218  praut(i,k) = praut(i,k)*factor
1219  pracw(i,k) = pracw(i,k)*factor
1220  paacw(i,k) = paacw(i,k)*factor
1221  endif
1222 !
1223 ! cloud ice
1224 !
1225  value = max(qmin,qi(i,k))
1226  source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
1227  + pgaci(i,k))*dtcld
1228  if (source.gt.value) then
1229  factor = value/source
1230  psaut(i,k) = psaut(i,k)*factor
1231  pigen(i,k) = pigen(i,k)*factor
1232  pidep(i,k) = pidep(i,k)*factor
1233  praci(i,k) = praci(i,k)*factor
1234  psaci(i,k) = psaci(i,k)*factor
1235  pgaci(i,k) = pgaci(i,k)*factor
1236  endif
1237 !
1238 ! rain
1239 !
1240  value = max(qmin,qr(i,k))
1241  source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)+psacr(i,k) &
1242  + pgacr(i,k))*dtcld
1243  if (source.gt.value) then
1244  factor = value/source
1245  praut(i,k) = praut(i,k)*factor
1246  prevp(i,k) = prevp(i,k)*factor
1247  pracw(i,k) = pracw(i,k)*factor
1248  piacr(i,k) = piacr(i,k)*factor
1249  psacr(i,k) = psacr(i,k)*factor
1250  pgacr(i,k) = pgacr(i,k)*factor
1251  endif
1252 !
1253 ! snow
1254 !
1255  value = max(qmin,qs(i,k))
1256  source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)+piacr(i,k) &
1257  * delta3+praci(i,k)*delta3-pracs(i,k)*(1.-delta2) &
1258  + psacr(i,k)*delta2+psaci(i,k)-pgacs(i,k) )*dtcld
1259  if (source.gt.value) then
1260  factor = value/source
1261  psdep(i,k) = psdep(i,k)*factor
1262  psaut(i,k) = psaut(i,k)*factor
1263  pgaut(i,k) = pgaut(i,k)*factor
1264  paacw(i,k) = paacw(i,k)*factor
1265  piacr(i,k) = piacr(i,k)*factor
1266  praci(i,k) = praci(i,k)*factor
1267  psaci(i,k) = psaci(i,k)*factor
1268  pracs(i,k) = pracs(i,k)*factor
1269  psacr(i,k) = psacr(i,k)*factor
1270  pgacs(i,k) = pgacs(i,k)*factor
1271  endif
1272 !
1273 ! graupel
1274 !
1275  value = max(qmin,qg(i,k))
1276  source = -(pgdep(i,k)+pgaut(i,k) &
1277  + piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
1278  + psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
1279  + pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1280  if (source.gt.value) then
1281  factor = value/source
1282  pgdep(i,k) = pgdep(i,k)*factor
1283  pgaut(i,k) = pgaut(i,k)*factor
1284  piacr(i,k) = piacr(i,k)*factor
1285  praci(i,k) = praci(i,k)*factor
1286  psacr(i,k) = psacr(i,k)*factor
1287  pracs(i,k) = pracs(i,k)*factor
1288  paacw(i,k) = paacw(i,k)*factor
1289  pgaci(i,k) = pgaci(i,k)*factor
1290  pgacr(i,k) = pgacr(i,k)*factor
1291  pgacs(i,k) = pgacs(i,k)*factor
1292  endif
1293 !
1294  work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
1295 ! update
1296  q(i,k) = q(i,k)+work2(i,k)*dtcld
1297  qc(i,k) = max(qc(i,k)-(praut(i,k)+pracw(i,k) &
1298  + paacw(i,k)+paacw(i,k))*dtcld,0.)
1299  qr(i,k) = max(qr(i,k)+(praut(i,k)+pracw(i,k) &
1300  + prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1301  - psacr(i,k))*dtcld,0.)
1302  qi(i,k) = max(qi(i,k)-(psaut(i,k)+praci(i,k) &
1303  + psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1304  * dtcld,0.)
1305  qs(i,k) = max(qs(i,k)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1306  - pgaut(i,k)+piacr(i,k)*delta3 &
1307  + praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1308  - pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
1309  * dtcld,0.)
1310  qg(i,k) = max(qg(i,k)+(pgdep(i,k)+pgaut(i,k) &
1311  + piacr(i,k)*(1.-delta3) &
1312  + praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) &
1313  + pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1314  + pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1315  xlf = xls-xl(i,k)
1316  xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1317  -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
1318  +paacw(i,k)+pgacr(i,k)+psacr(i,k))
1319  t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1320  else
1321 !
1322 ! cloud water
1323 !
1324  value = max(qmin,qc(i,k))
1325  source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1326  if (source.gt.value) then
1327  factor = value/source
1328  praut(i,k) = praut(i,k)*factor
1329  pracw(i,k) = pracw(i,k)*factor
1330  paacw(i,k) = paacw(i,k)*factor
1331  endif
1332 !
1333 ! rain
1334 !
1335  value = max(qmin,qr(i,k))
1336  source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)-pracw(i,k) &
1337  -paacw(i,k)-prevp(i,k))*dtcld
1338  if (source.gt.value) then
1339  factor = value/source
1340  praut(i,k) = praut(i,k)*factor
1341  prevp(i,k) = prevp(i,k)*factor
1342  pracw(i,k) = pracw(i,k)*factor
1343  paacw(i,k) = paacw(i,k)*factor
1344  pseml(i,k) = pseml(i,k)*factor
1345  pgeml(i,k) = pgeml(i,k)*factor
1346  endif
1347 !
1348 ! snow
1349 !
1350  value = max(qcrmin,qs(i,k))
1351  source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1352  if (source.gt.value) then
1353  factor = value/source
1354  pgacs(i,k) = pgacs(i,k)*factor
1355  psevp(i,k) = psevp(i,k)*factor
1356  pseml(i,k) = pseml(i,k)*factor
1357  endif
1358 !
1359 ! graupel
1360 !
1361  value = max(qcrmin,qg(i,k))
1362  source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1363  if (source.gt.value) then
1364  factor = value/source
1365  pgacs(i,k) = pgacs(i,k)*factor
1366  pgevp(i,k) = pgevp(i,k)*factor
1367  pgeml(i,k) = pgeml(i,k)*factor
1368  endif
1369 !
1370  work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1371 ! update
1372  q(i,k) = q(i,k)+work2(i,k)*dtcld
1373  qc(i,k) = max(qc(i,k)-(praut(i,k)+pracw(i,k) &
1374  + paacw(i,k)+paacw(i,k))*dtcld,0.)
1375  qr(i,k) = max(qr(i,k)+(praut(i,k)+pracw(i,k) &
1376  + prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1377  - pgeml(i,k))*dtcld,0.)
1378  qs(i,k) = max(qs(i,k)+(psevp(i,k)-pgacs(i,k) &
1379  + pseml(i,k))*dtcld,0.)
1380  qg(i,k) = max(qg(i,k)+(pgacs(i,k)+pgevp(i,k) &
1381  + pgeml(i,k))*dtcld,0.)
1382  xlf = xls-xl(i,k)
1383  xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1384  -xlf*(pseml(i,k)+pgeml(i,k))
1385  t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1386  endif
1387  enddo
1388  enddo
1389 !
1390 ! Inline expansion for fpvs
1391 ! qsat(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1392 ! qsat(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1393  hsub = xls
1394  hvap = xlv0
1395  cvap = cpv
1396  ttp=t0c+0.01
1397  dldt=cvap-cliq
1398  xa=-dldt/rv
1399  xb=xa+hvap/(rv*ttp)
1400  dldti=cvap-cice
1401  xai=-dldti/rv
1402  xbi=xai+hsub/(rv*ttp)
1403  do k = kts, kte
1404  do i = its, ite
1405  tr=ttp/t(i,k)
1406  qsat(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1407  qsat(i,k,1) = min(qsat(i,k,1),0.99*p(i,k))
1408  qsat(i,k,1) = ep2 * qsat(i,k,1) / (p(i,k) - qsat(i,k,1))
1409  qsat(i,k,1) = max(qsat(i,k,1),qmin)
1410  tr=ttp/t(i,k)
1411  if(t(i,k).lt.ttp) then
1412  qsat(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1413  else
1414  qsat(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1415  endif
1416  qsat(i,k,2) = min(qsat(i,k,2),0.99*p(i,k))
1417  qsat(i,k,2) = ep2 * qsat(i,k,2) / (p(i,k) - qsat(i,k,2))
1418  qsat(i,k,2) = max(qsat(i,k,2),qmin)
1419  enddo
1420  enddo
1421 !
1422 !----------------------------------------------------------------
1423 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1424 ! if there exists additional water vapor condensated/if
1425 ! evaporation of cloud water is not enough to remove subsaturation
1426 !
1427  do k = kts, kte
1428  do i = its, ite
1429  work1(i,k,1) = conden(t(i,k),q(i,k),qsat(i,k,1),xl(i,k),cpm(i,k))
1430  work2(i,k) = qc(i,k)+work1(i,k,1)
1431  pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1432  if(qc(i,k).gt.0..and.work1(i,k,1).lt.0.) &
1433  pcond(i,k) = max(work1(i,k,1),-qc(i,k))/dtcld
1434  q(i,k) = q(i,k)-pcond(i,k)*dtcld
1435  qc(i,k) = max(qc(i,k)+pcond(i,k)*dtcld,0.)
1436  t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1437  enddo
1438  enddo
1439 !
1440 !
1441 !----------------------------------------------------------------
1442 ! padding for small values
1443 !
1444  do k = kts, kte
1445  do i = its, ite
1446  if(qc(i,k).le.qmin) qc(i,k) = 0.0
1447  if(qi(i,k).le.qmin) qi(i,k) = 0.0
1448  enddo
1449  enddo
1450  enddo ! big loops
1451 
1452  if(present(rainprod2d) .and. present(evapprod2d)) then
1453  do k = kts, kte
1454  do i = its,ite
1455  rainprod2d(i,k) = praut(i,k)+pracw(i,k)+praci(i,k)+psaci(i,k)+pgaci(i,k) &
1456  + psacw(i,k)+pgacw(i,k)+paacw(i,k)+psaut(i,k)
1457  evapprod2d(i,k) = -(prevp(i,k)+psevp(i,k)+pgevp(i,k)+psdep(i,k)+pgdep(i,k))
1458  enddo
1459  enddo
1460  endif
1461 !
1462 !----------------------------------------------------------------
1463 ! CCPP checks:
1464 !
1465 
1466  errmsg = 'mp_wsm6_run OK'
1467  errflg = 0
1468 

Referenced by mp_wsm6_isohelper::mp_wsm6_run_c().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ nislfv_rain_plm()

subroutine mp_wsm6::nislfv_rain_plm ( integer, intent(in)  im,
integer, intent(in)  km,
real(kind=kind_phys), dimension(im,km), intent(in)  denl,
real(kind=kind_phys), dimension(im,km), intent(in)  denfacl,
real(kind=kind_phys), dimension(im,km), intent(in)  tkl,
real(kind=kind_phys), dimension(im,km), intent(in)  dzl,
real(kind=kind_phys), dimension(im,km), intent(inout)  wwl,
real(kind=kind_phys), dimension(im,km), intent(inout)  rql,
real(kind=kind_phys), dimension(im), intent(inout)  precip,
real(kind=kind_phys), intent(in)  dt,
integer, intent(in)  id,
integer, intent(in)  iter 
)
private
1753 !=================================================================================================================
1754 !
1755 ! for non-iteration semi-Lagrangain forward advection for cloud
1756 ! with mass conservation and positive definite advection
1757 ! 2nd order interpolation with monotonic piecewise linear method
1758 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1759 !
1760 ! dzl depth of model layer in meter
1761 ! wwl terminal velocity at model layer m/s
1762 ! rql cloud density*mixing ration
1763 ! precip precipitation
1764 ! dt time step
1765 ! id kind of precip: 0 test case; 1 raindrop
1766 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1767 ! 0 : use departure wind for advection
1768 ! 1 : use mean wind for advection
1769 ! > 1 : use mean wind after iter-1 iterations
1770 !
1771 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1772 ! implemented by song-you hong
1773 !
1774 
1775 !--- input arguments:
1776  integer,intent(in):: im,km,id,iter
1777 
1778  real(kind=kind_phys),intent(in):: dt
1779  real(kind=kind_phys),intent(in),dimension(im,km):: dzl,denl,denfacl,tkl
1780 
1781 !--- inout arguments:
1782  real(kind=kind_phys),intent(inout),dimension(im):: precip
1783  real(kind=kind_phys),intent(inout),dimension(im,km):: rql,wwl
1784 
1785 !---- local variables and arrays:
1786  integer:: i,k,n,m,kk,kb,kt
1787  real(kind=kind_phys):: tl,tl2,qql,dql,qqd
1788  real(kind=kind_phys):: th,th2,qqh,dqh
1789  real(kind=kind_phys):: zsum,qsum,dim,dip,c1,con1,fa1,fa2
1790  real(kind=kind_phys):: allold,allnew,zz,dzamin,cflmax,decfl
1791  real(kind=kind_phys),dimension(km):: dz,ww,qq,wd,wa,was
1792  real(kind=kind_phys),dimension(km):: den,denfac,tk
1793  real(kind=kind_phys),dimension(km):: qn,qr,tmp,tmp1,tmp2,tmp3
1794  real(kind=kind_phys),dimension(km+1):: wi,zi,za
1795  real(kind=kind_phys),dimension(km+1):: dza,qa,qmi,qpi
1796 
1797 !-----------------------------------------------------------------------------------------------------------------
1798 
1799  precip(:) = 0.0
1800 
1801  i_loop: do i=1,im
1802  dz(:) = dzl(i,:)
1803  qq(:) = rql(i,:)
1804  ww(:) = wwl(i,:)
1805  den(:) = denl(i,:)
1806  denfac(:) = denfacl(i,:)
1807  tk(:) = tkl(i,:)
1808 ! skip for no precipitation for all layers
1809  allold = 0.0
1810  do k=1,km
1811  allold = allold + qq(k)
1812  enddo
1813  if(allold.le.0.0) then
1814  cycle i_loop
1815  endif
1816 !
1817 ! compute interface values
1818  zi(1)=0.0
1819  do k=1,km
1820  zi(k+1) = zi(k)+dz(k)
1821  enddo
1822 !
1823 ! save departure wind
1824  wd(:) = ww(:)
1825  n=1
1826  100 continue
1827 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1828 ! 2nd order interpolation to get wi
1829  wi(1) = ww(1)
1830  wi(km+1) = ww(km)
1831  do k=2,km
1832  wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1833  enddo
1834 ! 3rd order interpolation to get wi
1835  fa1 = 9./16.
1836  fa2 = 1./16.
1837  wi(1) = ww(1)
1838  wi(2) = 0.5*(ww(2)+ww(1))
1839  do k=3,km-1
1840  wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1841  enddo
1842  wi(km) = 0.5*(ww(km)+ww(km-1))
1843  wi(km+1) = ww(km)
1844 !
1845 ! terminate of top of raingroup
1846  do k=2,km
1847  if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1848  enddo
1849 !
1850 ! diffusivity of wi
1851  con1 = 0.05
1852  do k=km,1,-1
1853  decfl = (wi(k+1)-wi(k))*dt/dz(k)
1854  if( decfl .gt. con1 ) then
1855  wi(k) = wi(k+1) - con1*dz(k)/dt
1856  endif
1857  enddo
1858 ! compute arrival point
1859  do k=1,km+1
1860  za(k) = zi(k) - wi(k)*dt
1861  enddo
1862 !
1863  do k=1,km
1864  dza(k) = za(k+1)-za(k)
1865  enddo
1866  dza(km+1) = zi(km+1) - za(km+1)
1867 !
1868 ! computer deformation at arrival point
1869  do k=1,km
1870  qa(k) = qq(k)*dz(k)/dza(k)
1871  qr(k) = qa(k)/den(k)
1872  enddo
1873  qa(km+1) = 0.0
1874 ! call maxmin(km,1,qa,' arrival points ')
1875 !
1876 ! compute arrival terminal velocity, and estimate mean terminal velocity
1877 ! then back to use mean terminal velocity
1878  if( n.le.iter ) then
1879  call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1880  if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1881  do k=1,km
1882 !#ifdef DEBUG
1883 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1884 !#endif
1885 ! mean wind is average of departure and new arrival winds
1886  ww(k) = 0.5* ( wd(k)+wa(k) )
1887  enddo
1888  was(:) = wa(:)
1889  n=n+1
1890  go to 100
1891  endif
1892 !
1893 ! estimate values at arrival cell interface with monotone
1894  do k=2,km
1895  dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1896  dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1897  if( dip*dim.le.0.0 ) then
1898  qmi(k)=qa(k)
1899  qpi(k)=qa(k)
1900  else
1901  qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1902  qmi(k)=2.0*qa(k)-qpi(k)
1903  if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1904  qpi(k) = qa(k)
1905  qmi(k) = qa(k)
1906  endif
1907  endif
1908  enddo
1909  qpi(1)=qa(1)
1910  qmi(1)=qa(1)
1911  qmi(km+1)=qa(km+1)
1912  qpi(km+1)=qa(km+1)
1913 !
1914 ! interpolation to regular point
1915  qn = 0.0
1916  kb=1
1917  kt=1
1918  intp : do k=1,km
1919  kb=max(kb-1,1)
1920  kt=max(kt-1,1)
1921 ! find kb and kt
1922  if( zi(k).ge.za(km+1) ) then
1923  exit intp
1924  else
1925  find_kb : do kk=kb,km
1926  if( zi(k).le.za(kk+1) ) then
1927  kb = kk
1928  exit find_kb
1929  else
1930  cycle find_kb
1931  endif
1932  enddo find_kb
1933  find_kt : do kk=kt,km
1934  if( zi(k+1).le.za(kk) ) then
1935  kt = kk
1936  exit find_kt
1937  else
1938  cycle find_kt
1939  endif
1940  enddo find_kt
1941  kt = kt - 1
1942 ! compute q with piecewise constant method
1943  if( kt.eq.kb ) then
1944  tl=(zi(k)-za(kb))/dza(kb)
1945  th=(zi(k+1)-za(kb))/dza(kb)
1946  tl2=tl*tl
1947  th2=th*th
1948  qqd=0.5*(qpi(kb)-qmi(kb))
1949  qqh=qqd*th2+qmi(kb)*th
1950  qql=qqd*tl2+qmi(kb)*tl
1951  qn(k) = (qqh-qql)/(th-tl)
1952  else if( kt.gt.kb ) then
1953  tl=(zi(k)-za(kb))/dza(kb)
1954  tl2=tl*tl
1955  qqd=0.5*(qpi(kb)-qmi(kb))
1956  qql=qqd*tl2+qmi(kb)*tl
1957  dql = qa(kb)-qql
1958  zsum = (1.-tl)*dza(kb)
1959  qsum = dql*dza(kb)
1960  if( kt-kb.gt.1 ) then
1961  do m=kb+1,kt-1
1962  zsum = zsum + dza(m)
1963  qsum = qsum + qa(m) * dza(m)
1964  enddo
1965  endif
1966  th=(zi(k+1)-za(kt))/dza(kt)
1967  th2=th*th
1968  qqd=0.5*(qpi(kt)-qmi(kt))
1969  dqh=qqd*th2+qmi(kt)*th
1970  zsum = zsum + th*dza(kt)
1971  qsum = qsum + dqh*dza(kt)
1972  qn(k) = qsum/zsum
1973  endif
1974  cycle intp
1975  endif
1976 !
1977  enddo intp
1978 !
1979 ! rain out
1980  sum_precip: do k=1,km
1981  if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1982  precip(i) = precip(i) + qa(k)*dza(k)
1983  cycle sum_precip
1984  else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1985  precip(i) = precip(i) + qa(k)*(0.0-za(k))
1986  exit sum_precip
1987  endif
1988  exit sum_precip
1989  enddo sum_precip
1990 !
1991 ! replace the new values
1992  rql(i,:) = qn(:)
1993  enddo i_loop
1994 

Referenced by mp_wsm6_run().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ nislfv_rain_plm6()

subroutine mp_wsm6::nislfv_rain_plm6 ( integer, intent(in)  im,
integer, intent(in)  km,
real(kind=kind_phys), dimension(im,km), intent(in)  denl,
real(kind=kind_phys), dimension(im,km), intent(in)  denfacl,
real(kind=kind_phys), dimension(im,km), intent(in)  tkl,
real(kind=kind_phys), dimension(im,km), intent(in)  dzl,
real(kind=kind_phys), dimension(im,km), intent(inout)  wwl,
real(kind=kind_phys), dimension(im,km), intent(inout)  rql,
real(kind=kind_phys), dimension(im,km), intent(inout)  rql2,
real(kind=kind_phys), dimension(im), intent(inout)  precip1,
real(kind=kind_phys), dimension(im), intent(inout)  precip2,
real(kind=kind_phys), intent(in)  dt,
integer, intent(in)  id,
integer, intent(in)  iter 
)
private
1999 !=================================================================================================================
2000 !
2001 ! for non-iteration semi-Lagrangain forward advection for cloud
2002 ! with mass conservation and positive definite advection
2003 ! 2nd order interpolation with monotonic piecewise linear method
2004 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
2005 !
2006 ! dzl depth of model layer in meter
2007 ! wwl terminal velocity at model layer m/s
2008 ! rql cloud density*mixing ration
2009 ! precip precipitation
2010 ! dt time step
2011 ! id kind of precip: 0 test case; 1 raindrop
2012 ! iter how many time to guess mean terminal velocity: 0 pure forward.
2013 ! 0 : use departure wind for advection
2014 ! 1 : use mean wind for advection
2015 ! > 1 : use mean wind after iter-1 iterations
2016 !
2017 ! author: hann-ming henry juang <henry.juang@noaa.gov>
2018 ! implemented by song-you hong
2019 !
2020 
2021 !--- input arguments:
2022  integer,intent(in):: im,km,id,iter
2023 
2024  real(kind=kind_phys),intent(in):: dt
2025  real(kind=kind_phys),intent(in),dimension(im,km):: dzl,denl,denfacl,tkl
2026 
2027 !--- inout arguments:
2028  real(kind=kind_phys),intent(inout),dimension(im):: precip1,precip2
2029  real(kind=kind_phys),intent(inout),dimension(im,km):: rql,rql2,wwl
2030 
2031 !---- local variables and arrays:
2032  integer:: i,ist,k,n,m,kk,kb,kt
2033  real(kind=kind_phys):: tl,tl2,qql,dql,qqd
2034  real(kind=kind_phys):: th,th2,qqh,dqh
2035  real(kind=kind_phys):: zsum,qsum,dim,dip,c1,con1,fa1,fa2
2036  real(kind=kind_phys):: allold,allnew,zz,dzamin,cflmax,decfl
2037  real(kind=kind_phys),dimension(km):: dz,ww,qq,qq2,wd,wa,wa2,was
2038  real(kind=kind_phys),dimension(km):: den,denfac,tk
2039  real(kind=kind_phys),dimension(km):: qn,qr,qr2,tmp,tmp1,tmp2,tmp3
2040  real(kind=kind_phys),dimension(km+1):: wi,zi,za
2041  real(kind=kind_phys),dimension(km+1):: dza,qa,qa2,qmi,qpi
2042  real(kind=kind_phys),dimension(im):: precip
2043 
2044 !-----------------------------------------------------------------------------------------------------------------
2045 
2046  precip(:) = 0.0
2047  precip1(:) = 0.0
2048  precip2(:) = 0.0
2049 
2050  i_loop: do i=1,im
2051  dz(:) = dzl(i,:)
2052  qq(:) = rql(i,:)
2053  qq2(:) = rql2(i,:)
2054  ww(:) = wwl(i,:)
2055  den(:) = denl(i,:)
2056  denfac(:) = denfacl(i,:)
2057  tk(:) = tkl(i,:)
2058 ! skip for no precipitation for all layers
2059  allold = 0.0
2060  do k=1,km
2061  allold = allold + qq(k) + qq2(k)
2062  enddo
2063  if(allold.le.0.0) then
2064  cycle i_loop
2065  endif
2066 !
2067 ! compute interface values
2068  zi(1)=0.0
2069  do k=1,km
2070  zi(k+1) = zi(k)+dz(k)
2071  enddo
2072 !
2073 ! save departure wind
2074  wd(:) = ww(:)
2075  n=1
2076  100 continue
2077 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2078 ! 2nd order interpolation to get wi
2079  wi(1) = ww(1)
2080  wi(km+1) = ww(km)
2081  do k=2,km
2082  wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
2083  enddo
2084 ! 3rd order interpolation to get wi
2085  fa1 = 9./16.
2086  fa2 = 1./16.
2087  wi(1) = ww(1)
2088  wi(2) = 0.5*(ww(2)+ww(1))
2089  do k=3,km-1
2090  wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
2091  enddo
2092  wi(km) = 0.5*(ww(km)+ww(km-1))
2093  wi(km+1) = ww(km)
2094 !
2095 ! terminate of top of raingroup
2096  do k=2,km
2097  if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2098  enddo
2099 !
2100 ! diffusivity of wi
2101  con1 = 0.05
2102  do k=km,1,-1
2103  decfl = (wi(k+1)-wi(k))*dt/dz(k)
2104  if( decfl .gt. con1 ) then
2105  wi(k) = wi(k+1) - con1*dz(k)/dt
2106  endif
2107  enddo
2108 ! compute arrival point
2109  do k=1,km+1
2110  za(k) = zi(k) - wi(k)*dt
2111  enddo
2112 !
2113  do k=1,km
2114  dza(k) = za(k+1)-za(k)
2115  enddo
2116  dza(km+1) = zi(km+1) - za(km+1)
2117 !
2118 ! computer deformation at arrival point
2119  do k=1,km
2120  qa(k) = qq(k)*dz(k)/dza(k)
2121  qa2(k) = qq2(k)*dz(k)/dza(k)
2122  qr(k) = qa(k)/den(k)
2123  qr2(k) = qa2(k)/den(k)
2124  enddo
2125  qa(km+1) = 0.0
2126  qa2(km+1) = 0.0
2127 ! call maxmin(km,1,qa,' arrival points ')
2128 !
2129 ! compute arrival terminal velocity, and estimate mean terminal velocity
2130 ! then back to use mean terminal velocity
2131  if( n.le.iter ) then
2132  call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2133  call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
2134  do k = 1, km
2135  tmp(k) = max((qr(k)+qr2(k)), 1.e-15)
2136  if( tmp(k) .gt. 1.e-15 ) then
2137  wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
2138  else
2139  wa(k) = 0.
2140  endif
2141  enddo
2142  if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2143  do k=1,km
2144 !#ifdef DEBUG
2145 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
2146 ! ww(k),wa(k)
2147 !#endif
2148 ! mean wind is average of departure and new arrival winds
2149  ww(k) = 0.5* ( wd(k)+wa(k) )
2150  enddo
2151  was(:) = wa(:)
2152  n=n+1
2153  go to 100
2154  endif
2155 
2156  ist_loop : do ist = 1, 2
2157  if (ist.eq.2) then
2158  qa(:) = qa2(:)
2159  endif
2160 !
2161  precip(i) = 0.
2162 !
2163 ! estimate values at arrival cell interface with monotone
2164  do k=2,km
2165  dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2166  dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2167  if( dip*dim.le.0.0 ) then
2168  qmi(k)=qa(k)
2169  qpi(k)=qa(k)
2170  else
2171  qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2172  qmi(k)=2.0*qa(k)-qpi(k)
2173  if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2174  qpi(k) = qa(k)
2175  qmi(k) = qa(k)
2176  endif
2177  endif
2178  enddo
2179  qpi(1)=qa(1)
2180  qmi(1)=qa(1)
2181  qmi(km+1)=qa(km+1)
2182  qpi(km+1)=qa(km+1)
2183 !
2184 ! interpolation to regular point
2185  qn = 0.0
2186  kb=1
2187  kt=1
2188  intp : do k=1,km
2189  kb=max(kb-1,1)
2190  kt=max(kt-1,1)
2191 ! find kb and kt
2192  if( zi(k).ge.za(km+1) ) then
2193  exit intp
2194  else
2195  find_kb : do kk=kb,km
2196  if( zi(k).le.za(kk+1) ) then
2197  kb = kk
2198  exit find_kb
2199  else
2200  cycle find_kb
2201  endif
2202  enddo find_kb
2203  find_kt : do kk=kt,km
2204  if( zi(k+1).le.za(kk) ) then
2205  kt = kk
2206  exit find_kt
2207  else
2208  cycle find_kt
2209  endif
2210  enddo find_kt
2211  kt = kt - 1
2212 ! compute q with piecewise constant method
2213  if( kt.eq.kb ) then
2214  tl=(zi(k)-za(kb))/dza(kb)
2215  th=(zi(k+1)-za(kb))/dza(kb)
2216  tl2=tl*tl
2217  th2=th*th
2218  qqd=0.5*(qpi(kb)-qmi(kb))
2219  qqh=qqd*th2+qmi(kb)*th
2220  qql=qqd*tl2+qmi(kb)*tl
2221  qn(k) = (qqh-qql)/(th-tl)
2222  else if( kt.gt.kb ) then
2223  tl=(zi(k)-za(kb))/dza(kb)
2224  tl2=tl*tl
2225  qqd=0.5*(qpi(kb)-qmi(kb))
2226  qql=qqd*tl2+qmi(kb)*tl
2227  dql = qa(kb)-qql
2228  zsum = (1.-tl)*dza(kb)
2229  qsum = dql*dza(kb)
2230  if( kt-kb.gt.1 ) then
2231  do m=kb+1,kt-1
2232  zsum = zsum + dza(m)
2233  qsum = qsum + qa(m) * dza(m)
2234  enddo
2235  endif
2236  th=(zi(k+1)-za(kt))/dza(kt)
2237  th2=th*th
2238  qqd=0.5*(qpi(kt)-qmi(kt))
2239  dqh=qqd*th2+qmi(kt)*th
2240  zsum = zsum + th*dza(kt)
2241  qsum = qsum + dqh*dza(kt)
2242  qn(k) = qsum/zsum
2243  endif
2244  cycle intp
2245  endif
2246 !
2247  enddo intp
2248 !
2249 ! rain out
2250  sum_precip: do k=1,km
2251  if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2252  precip(i) = precip(i) + qa(k)*dza(k)
2253  cycle sum_precip
2254  else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2255  precip(i) = precip(i) + qa(k)*(0.0-za(k))
2256  exit sum_precip
2257  endif
2258  exit sum_precip
2259  enddo sum_precip
2260 !
2261 ! replace the new values
2262  if(ist.eq.1) then
2263  rql(i,:) = qn(:)
2264  precip1(i) = precip(i)
2265  else
2266  rql2(i,:) = qn(:)
2267  precip2(i) = precip(i)
2268  endif
2269  enddo ist_loop
2270 
2271  enddo i_loop
2272 

Referenced by mp_wsm6_run().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ refl10cm_wsm6()

subroutine, public mp_wsm6::refl10cm_wsm6 ( real(kind=kind_phys), dimension(kts:kte), intent(in)  qv1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  qr1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  qs1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  qg1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  t1d,
real(kind=kind_phys), dimension(kts:kte), intent(in)  p1d,
real(kind=kind_phys), dimension(kts:kte), intent(inout)  dBZ,
integer, intent(in)  kts,
integer, intent(in)  kte 
)
2277  implicit none
2278 !=================================================================================================================
2279 
2280 !..Sub arguments
2281  integer,intent(in):: kts,kte
2282  real(kind=kind_phys),intent(in),dimension(kts:kte):: qv1d,qr1d,qs1d,qg1d,t1d,p1d
2283  real(kind=kind_phys),intent(inout),dimension(kts:kte):: dbz
2284 
2285 !..Local variables
2286  logical:: melti
2287  logical,dimension(kts:kte):: l_qr,l_qs,l_qg
2288 
2289  INTEGER:: i,k,k_0,kbot,n
2290 
2291  real(kind=kind_phys),parameter:: r=287.
2292  real(kind=kind_phys):: temp_c
2293  real(kind=kind_phys),dimension(kts:kte):: temp,pres,qv,rho
2294  real(kind=kind_phys),dimension(kts:kte):: rr,rs,rg
2295  real(kind=kind_phys),dimension(kts:kte):: ze_rain,ze_snow,ze_graupel
2296 
2297  double precision:: fmelt_s,fmelt_g
2298  double precision:: cback,x,eta,f_d
2299  double precision,dimension(kts:kte):: ilamr,ilams,ilamg
2300  double precision,dimension(kts:kte):: n0_r, n0_s, n0_g
2301  double precision:: lamr,lams,lamg
2302 
2303 !-----------------------------------------------------------------------------------------------------------------
2304 
2305  do k = kts, kte
2306  dbz(k) = -35.0
2307  enddo
2308 
2309 !+---+-----------------------------------------------------------------+
2310 !..Put column of data into local arrays.
2311 !+---+-----------------------------------------------------------------+
2312  do k = kts, kte
2313  temp(k) = t1d(k)
2314  temp_c = min(-0.001, temp(k)-273.15)
2315  qv(k) = max(1.e-10, qv1d(k))
2316  pres(k) = p1d(k)
2317  rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
2318 
2319  if (qr1d(k) .gt. 1.e-9) then
2320  rr(k) = qr1d(k)*rho(k)
2321  n0_r(k) = n0r
2322  lamr = (xam_r*xcrg(3)*n0_r(k)/rr(k))**(1./xcre(1))
2323  ilamr(k) = 1./lamr
2324  l_qr(k) = .true.
2325  else
2326  rr(k) = 1.e-12
2327  l_qr(k) = .false.
2328  endif
2329 
2330  if (qs1d(k) .gt. 1.e-9) then
2331  rs(k) = qs1d(k)*rho(k)
2332  n0_s(k) = min(n0smax, n0s*exp(-alpha*temp_c))
2333  lams = (xam_s*xcsg(3)*n0_s(k)/rs(k))**(1./xcse(1))
2334  ilams(k) = 1./lams
2335  l_qs(k) = .true.
2336  else
2337  rs(k) = 1.e-12
2338  l_qs(k) = .false.
2339  endif
2340 
2341  if (qg1d(k) .gt. 1.e-9) then
2342  rg(k) = qg1d(k)*rho(k)
2343  n0_g(k) = n0g
2344  lamg = (xam_g*xcgg(3)*n0_g(k)/rg(k))**(1./xcge(1))
2345  ilamg(k) = 1./lamg
2346  l_qg(k) = .true.
2347  else
2348  rg(k) = 1.e-12
2349  l_qg(k) = .false.
2350  endif
2351  enddo
2352 
2353 !+---+-----------------------------------------------------------------+
2354 !..Locate K-level of start of melting (k_0 is level above).
2355 !+---+-----------------------------------------------------------------+
2356  melti = .false.
2357  k_0 = kts
2358  do k = kte-1, kts, -1
2359  if ( (temp(k).gt.273.15) .and. l_qr(k) &
2360  .and. (l_qs(k+1).or.l_qg(k+1)) ) then
2361  k_0 = max(k+1, k_0)
2362  melti=.true.
2363  goto 195
2364  endif
2365  enddo
2366  195 continue
2367 
2368 !+---+-----------------------------------------------------------------+
2369 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
2370 !.. and non-water-coated snow and graupel when below freezing are
2371 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
2372 !+---+-----------------------------------------------------------------+
2373 
2374  do k = kts, kte
2375  ze_rain(k) = 1.e-22
2376  ze_snow(k) = 1.e-22
2377  ze_graupel(k) = 1.e-22
2378  if (l_qr(k)) ze_rain(k) = n0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
2379  if (l_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
2380  * (xam_s/900.0)*(xam_s/900.0) &
2381  * n0_s(k)*xcsg(4)*ilams(k)**xcse(4)
2382  if (l_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
2383  * (xam_g/900.0)*(xam_g/900.0) &
2384  * n0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
2385  enddo
2386 
2387 
2388 !+---+-----------------------------------------------------------------+
2389 !..Special case of melting ice (snow/graupel) particles. Assume the
2390 !.. ice is surrounded by the liquid water. Fraction of meltwater is
2391 !.. extremely simple based on amount found above the melting level.
2392 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
2393 !.. routines).
2394 !+---+-----------------------------------------------------------------+
2395 
2396  if (melti .and. k_0.ge.kts+1) then
2397  do k = k_0-1, kts, -1
2398 
2399 !..Reflectivity contributed by melting snow
2400  if (l_qs(k) .and. l_qs(k_0) ) then
2401  fmelt_s = max(0.005d0, min(1.0d0-rs(k)/rs(k_0), 0.99d0))
2402  eta = 0.d0
2403  lams = 1./ilams(k)
2404  do n = 1, nrbins
2405  x = xam_s * xxds(n)**xbm_s
2406  call rayleigh_soak_wetgraupel (x,dble(xocms),dble(xobms), &
2407  fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
2408  cback, mixingrulestring_s, matrixstring_s, &
2409  inclusionstring_s, hoststring_s, &
2410  hostmatrixstring_s, hostinclusionstring_s)
2411  f_d = n0_s(k)*xxds(n)**xmu_s * dexp(-lams*xxds(n))
2412  eta = eta + f_d * cback * simpson(n) * xdts(n)
2413  enddo
2414  ze_snow(k) = sngl(lamda4 / (pi5 * k_w) * eta)
2415  endif
2416 
2417 
2418 !..Reflectivity contributed by melting graupel
2419 
2420  if (l_qg(k) .and. l_qg(k_0) ) then
2421  fmelt_g = max(0.005d0, min(1.0d0-rg(k)/rg(k_0), 0.99d0))
2422  eta = 0.d0
2423  lamg = 1./ilamg(k)
2424  do n = 1, nrbins
2425  x = xam_g * xxdg(n)**xbm_g
2426  call rayleigh_soak_wetgraupel (x,dble(xocmg),dble(xobmg), &
2427  fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
2428  cback, mixingrulestring_g, matrixstring_g, &
2429  inclusionstring_g, hoststring_g, &
2430  hostmatrixstring_g, hostinclusionstring_g)
2431  f_d = n0_g(k)*xxdg(n)**xmu_g * dexp(-lamg*xxdg(n))
2432  eta = eta + f_d * cback * simpson(n) * xdtg(n)
2433  enddo
2434  ze_graupel(k) = sngl(lamda4 / (pi5 * k_w) * eta)
2435  endif
2436 
2437  enddo
2438  endif
2439 
2440  do k = kte, kts, -1
2441  dbz(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
2442  enddo
2443 
2444 

◆ rgmma()

real(kind=kind_phys) function mp_wsm6::rgmma ( real(kind=kind_phys), intent(in)  x)
private
1473 !=================================================================================================================
1474 !rgmma function: use infinite product form
1475 
1476  real(kind=kind_phys),intent(in):: x
1477 
1478  integer:: i
1479  real(kind=kind_phys),parameter:: euler=0.577215664901532
1480  real(kind=kind_phys):: y
1481 
1482 !-----------------------------------------------------------------------------------------------------------------
1483 
1484  if(x.eq.1.)then
1485  rgmma=0.
1486  else
1487  rgmma=x*exp(euler*x)
1488  do i = 1,10000
1489  y = float(i)
1490  rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1491  enddo
1492  rgmma=1./rgmma
1493  endif
1494 

Referenced by mp_wsm6_init().

Here is the caller graph for this function:

◆ slope_graup()

subroutine mp_wsm6::slope_graup ( real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1704 !=================================================================================================================
1705 
1706 !--- input arguments:
1707  integer:: its,ite,kts,kte
1708 
1709  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1710 
1711 !--- inout arguments:
1712  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1713  rslope,rslopeb,rslope2,rslope3,vt
1714 
1715 !--- local variables and arrays:
1716  integer:: i,k
1717 
1718  real(kind=kind_phys),parameter:: t0c = 273.15
1719  real(kind=kind_phys):: lamdag,x,y
1720  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1721 
1722 !-----------------------------------------------------------------------------------------------------------------
1723 
1724 !size distributions: (x=mixing ratio, y=air density):
1725 !valid for mixing ratio > 1.e-9 kg/kg.
1726  lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
1727 
1728  do k = kts, kte
1729  do i = its, ite
1730 !---------------------------------------------------------------
1731 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1732 !---------------------------------------------------------------
1733  if(qrs(i,k).le.qcrmin)then
1734  rslope(i,k) = rslopegmax
1735  rslopeb(i,k) = rslopegbmax
1736  rslope2(i,k) = rslopeg2max
1737  rslope3(i,k) = rslopeg3max
1738  else
1739  rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
1740  rslopeb(i,k) = rslope(i,k)**bvtg
1741  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1742  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1743  endif
1744  vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
1745  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1746  enddo
1747  enddo
1748 

Referenced by nislfv_rain_plm6().

Here is the caller graph for this function:

◆ slope_rain()

subroutine mp_wsm6::slope_rain ( real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1607 !=================================================================================================================
1608 
1609 !--- input arguments:
1610  integer:: its,ite,kts,kte
1611 
1612  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1613 
1614 !--- inout arguments:
1615  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1616  rslope,rslopeb,rslope2,rslope3,vt
1617 
1618 !--- local variables and arrays:
1619  integer:: i,k
1620 
1621  real(kind=kind_phys),parameter:: t0c = 273.15
1622  real(kind=kind_phys):: lamdar,x,y
1623  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1624 
1625 !-----------------------------------------------------------------------------------------------------------------
1626 
1627 !size distributions: (x=mixing ratio, y=air density):
1628 !valid for mixing ratio > 1.e-9 kg/kg.
1629  lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1630 
1631  do k = kts, kte
1632  do i = its, ite
1633  if(qrs(i,k).le.qcrmin)then
1634  rslope(i,k) = rslopermax
1635  rslopeb(i,k) = rsloperbmax
1636  rslope2(i,k) = rsloper2max
1637  rslope3(i,k) = rsloper3max
1638  else
1639  rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
1640  rslopeb(i,k) = rslope(i,k)**bvtr
1641  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1642  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1643  endif
1644  vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1645  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1646  enddo
1647  enddo
1648 

Referenced by nislfv_rain_plm().

Here is the caller graph for this function:

◆ slope_snow()

subroutine mp_wsm6::slope_snow ( real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1653 !=================================================================================================================
1654 
1655 !--- input arguments:
1656  integer:: its,ite,kts,kte
1657 
1658  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: qrs,den,denfac,t
1659 
1660 !--- inout arguments:
1661  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte):: &
1662  rslope,rslopeb,rslope2,rslope3,vt
1663 
1664 !--- local variables and arrays:
1665  integer:: i,k
1666 
1667  real(kind=kind_phys),parameter:: t0c = 273.15
1668  real(kind=kind_phys):: lamdas,x,y,z,supcol
1669  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1670 
1671 !-----------------------------------------------------------------------------------------------------------------
1672 
1673 !size distributions: (x=mixing ratio, y=air density):
1674 !valid for mixing ratio > 1.e-9 kg/kg.
1675  lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1676 !
1677  do k = kts, kte
1678  do i = its, ite
1679  supcol = t0c-t(i,k)
1680 !---------------------------------------------------------------
1681 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1682 !---------------------------------------------------------------
1683  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1684  if(qrs(i,k).le.qcrmin)then
1685  rslope(i,k) = rslopesmax
1686  rslopeb(i,k) = rslopesbmax
1687  rslope2(i,k) = rslopes2max
1688  rslope3(i,k) = rslopes3max
1689  else
1690  rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1691  rslopeb(i,k) = rslope(i,k)**bvts
1692  rslope2(i,k) = rslope(i,k)*rslope(i,k)
1693  rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1694  endif
1695  vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1696  if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1697  enddo
1698  enddo
1699 

Referenced by nislfv_rain_plm6().

Here is the caller graph for this function:

◆ slope_wsm6()

subroutine mp_wsm6::slope_wsm6 ( real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(in)  qrs,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  den,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  denfac,
real(kind=kind_phys), dimension(its:ite,kts:kte), intent(in)  t,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslope,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslopeb,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslope2,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  rslope3,
real(kind=kind_phys), dimension(its:ite,kts:kte,3), intent(inout)  vt,
integer  its,
integer  ite,
integer  kts,
integer  kte 
)
private
1527 !=================================================================================================================
1528 
1529 !--- input arguments:
1530  integer:: its,ite,kts,kte
1531 
1532  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte):: den,denfac,t
1533  real(kind=kind_phys),intent(in),dimension(its:ite,kts:kte,3):: qrs
1534 
1535 !--- inout arguments:
1536  real(kind=kind_phys),intent(inout),dimension(its:ite,kts:kte,3):: &
1537  rslope,rslopeb,rslope2,rslope3,vt
1538 
1539 !--- local variables and arrays:
1540  integer:: i,k
1541 
1542  real(kind=kind_phys),parameter:: t0c = 273.15
1543  real(kind=kind_phys):: lamdar,lamdas,lamdag,x,y,z,supcol
1544  real(kind=kind_phys),dimension(its:ite,kts:kte):: n0sfac
1545 
1546 !-----------------------------------------------------------------------------------------------------------------
1547 
1548 !size distributions: (x=mixing ratio, y=air density):
1549 !valid for mixing ratio > 1.e-9 kg/kg.
1550  lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1551  lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1552  lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
1553 
1554  do k = kts, kte
1555  do i = its, ite
1556  supcol = t0c-t(i,k)
1557 !---------------------------------------------------------------
1558 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1559 !---------------------------------------------------------------
1560  n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1561  if(qrs(i,k,1).le.qcrmin)then
1562  rslope(i,k,1) = rslopermax
1563  rslopeb(i,k,1) = rsloperbmax
1564  rslope2(i,k,1) = rsloper2max
1565  rslope3(i,k,1) = rsloper3max
1566  else
1567  rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
1568  rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1569  rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1570  rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1571  endif
1572  if(qrs(i,k,2).le.qcrmin)then
1573  rslope(i,k,2) = rslopesmax
1574  rslopeb(i,k,2) = rslopesbmax
1575  rslope2(i,k,2) = rslopes2max
1576  rslope3(i,k,2) = rslopes3max
1577  else
1578  rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1579  rslopeb(i,k,2) = rslope(i,k,2)**bvts
1580  rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1581  rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1582  endif
1583  if(qrs(i,k,3).le.qcrmin)then
1584  rslope(i,k,3) = rslopegmax
1585  rslopeb(i,k,3) = rslopegbmax
1586  rslope2(i,k,3) = rslopeg2max
1587  rslope3(i,k,3) = rslopeg3max
1588  else
1589  rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
1590  rslopeb(i,k,3) = rslope(i,k,3)**bvtg
1591  rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
1592  rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
1593  endif
1594  vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1595  vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1596  vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
1597  if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1598  if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1599  if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
1600  enddo
1601  enddo
1602 

Referenced by mp_wsm6_run().

Here is the caller graph for this function:

Variable Documentation

◆ alpha

real(kind=kind_phys), parameter, public mp_wsm6::alpha = .12
44  real(kind=kind_phys),parameter,public :: alpha = .12 ! .122 exponen factor for n0s

Referenced by compute_max_reflectivity_dbz(), SurfaceLayer::fill_tsurf_with_sst_and_tsk(), interpolate_1d(), mp_wsm6_run(), realbdy_compute_interior_ghost_rhs(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ avtg

real(kind=kind_phys), save mp_wsm6::avtg
private

Referenced by mp_wsm6_init().

◆ avtr

real(kind=kind_phys), parameter, private mp_wsm6::avtr = 841.9
private
19  real(kind=kind_phys),parameter,private:: avtr = 841.9 ! a constant for terminal velocity of rain

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ avts

real(kind=kind_phys), parameter, private mp_wsm6::avts = 11.72
private
25  real(kind=kind_phys),parameter,private:: avts = 11.72 ! a constant for terminal velocity of snow

Referenced by mp_wsm6_init().

◆ bvtg

real(kind=kind_phys), save mp_wsm6::bvtg
private

Referenced by mp_wsm6_init(), slope_graup(), and slope_wsm6().

◆ bvtg1

real(kind=kind_phys), save mp_wsm6::bvtg1
private

Referenced by mp_wsm6_init().

◆ bvtg2

real(kind=kind_phys), save mp_wsm6::bvtg2
private

Referenced by mp_wsm6_init().

◆ bvtg3

real(kind=kind_phys), save mp_wsm6::bvtg3
private

Referenced by mp_wsm6_init().

◆ bvtg4

real(kind=kind_phys), save mp_wsm6::bvtg4
private

Referenced by mp_wsm6_init().

◆ bvtr

real(kind=kind_phys), parameter, private mp_wsm6::bvtr = 0.8
private
20  real(kind=kind_phys),parameter,private:: bvtr = 0.8 ! a constant for terminal velocity of rain

Referenced by mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ bvtr1

real(kind=kind_phys), save mp_wsm6::bvtr1
private

Referenced by mp_wsm6_init().

◆ bvtr2

real(kind=kind_phys), save mp_wsm6::bvtr2
private

Referenced by mp_wsm6_init().

◆ bvtr3

real(kind=kind_phys), save mp_wsm6::bvtr3
private

Referenced by mp_wsm6_init().

◆ bvtr4

real(kind=kind_phys), save mp_wsm6::bvtr4
private

Referenced by mp_wsm6_init().

◆ bvtr6

real(kind=kind_phys), save mp_wsm6::bvtr6
private

Referenced by mp_wsm6_init().

◆ bvts

real(kind=kind_phys), parameter, private mp_wsm6::bvts = .41
private
26  real(kind=kind_phys),parameter,private:: bvts = .41 ! a constant for terminal velocity of snow

Referenced by mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ bvts1

real(kind=kind_phys), save mp_wsm6::bvts1
private

Referenced by mp_wsm6_init().

◆ bvts2

real(kind=kind_phys), save mp_wsm6::bvts2
private

Referenced by mp_wsm6_init().

◆ bvts3

real(kind=kind_phys), save mp_wsm6::bvts3
private

Referenced by mp_wsm6_init().

◆ bvts4

real(kind=kind_phys), save mp_wsm6::bvts4
private

Referenced by mp_wsm6_init().

◆ deng

real(kind=kind_phys), save mp_wsm6::deng
private

Referenced by mp_wsm6_init().

◆ dens

real(kind=kind_phys), parameter, private mp_wsm6::dens = 100.0
private
39  real(kind=kind_phys),parameter,private:: dens = 100.0 ! Density of snow

Referenced by WSM6::Advance(), ERF::erf_enforce_hse(), ERF::MakeHorizontalAverages(), mp_wsm6_init(), and mp_wsm6_run().

◆ dicon

real(kind=kind_phys), parameter, private mp_wsm6::dicon = 11.9
private
33  real(kind=kind_phys),parameter,private:: dicon = 11.9 ! constant for the cloud-ice diamter

Referenced by mp_wsm6_run().

◆ dimax

real(kind=kind_phys), parameter, private mp_wsm6::dimax = 500.e-6
private
34  real(kind=kind_phys),parameter,private:: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ dtcldcr

real(kind=kind_phys), parameter, private mp_wsm6::dtcldcr = 120.
private
16  real(kind=kind_phys),parameter,private:: dtcldcr = 120. ! maximum time step for minor loops

Referenced by mp_wsm6_run().

◆ eacrc

real(kind=kind_phys), parameter, private mp_wsm6::eacrc = 1.0
private
38  real(kind=kind_phys),parameter,private:: eacrc = 1.0 ! Snow/cloud-water collection efficiency

Referenced by mp_wsm6_init().

◆ eacrr

real(kind=kind_phys), save mp_wsm6::eacrr
private

Referenced by mp_wsm6_init().

◆ g1pbg

real(kind=kind_phys), save mp_wsm6::g1pbg
private

Referenced by mp_wsm6_init().

◆ g1pbr

real(kind=kind_phys), save mp_wsm6::g1pbr
private

Referenced by mp_wsm6_init().

◆ g1pbs

real(kind=kind_phys), save mp_wsm6::g1pbs
private

Referenced by mp_wsm6_init().

◆ g3pbg

real(kind=kind_phys), save mp_wsm6::g3pbg
private

Referenced by mp_wsm6_init().

◆ g3pbr

real(kind=kind_phys), save mp_wsm6::g3pbr
private

Referenced by mp_wsm6_init().

◆ g3pbs

real(kind=kind_phys), save mp_wsm6::g3pbs
private

Referenced by mp_wsm6_init().

◆ g4pbg

real(kind=kind_phys), save mp_wsm6::g4pbg
private

Referenced by mp_wsm6_init().

◆ g4pbr

real(kind=kind_phys), save mp_wsm6::g4pbr
private

Referenced by mp_wsm6_init().

◆ g4pbs

real(kind=kind_phys), save mp_wsm6::g4pbs
private

Referenced by mp_wsm6_init().

◆ g5pbgo2

real(kind=kind_phys), save mp_wsm6::g5pbgo2
private

Referenced by mp_wsm6_init().

◆ g5pbro2

real(kind=kind_phys), save mp_wsm6::g5pbro2
private

Referenced by mp_wsm6_init().

◆ g5pbso2

real(kind=kind_phys), save mp_wsm6::g5pbso2
private

Referenced by mp_wsm6_init().

◆ g6pbr

real(kind=kind_phys), save mp_wsm6::g6pbr
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ kind_phys

integer, parameter mp_wsm6::kind_phys = c_double
9  integer, parameter :: kind_phys = c_double

◆ lamdagmax

real(kind=kind_phys), save mp_wsm6::lamdagmax
private

Referenced by mp_wsm6_init().

◆ lamdarmax

real(kind=kind_phys), parameter, private mp_wsm6::lamdarmax = 8.e4
private
30  real(kind=kind_phys),parameter,private:: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain

Referenced by mp_wsm6_init().

◆ lamdasmax

real(kind=kind_phys), parameter, private mp_wsm6::lamdasmax = 1.e5
private
31  real(kind=kind_phys),parameter,private:: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow

Referenced by mp_wsm6_init().

◆ n0g

real(kind=kind_phys), save mp_wsm6::n0g
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), and refl10cm_wsm6().

◆ n0r

real(kind=kind_phys), parameter, private mp_wsm6::n0r = 8.e6
private
17  real(kind=kind_phys),parameter,private:: n0r = 8.e6 ! intercept parameter rain

Referenced by mp_wsm6_init(), mp_wsm6_run(), and refl10cm_wsm6().

◆ n0s

real(kind=kind_phys), parameter, public mp_wsm6::n0s = 2.e6
43  real(kind=kind_phys),parameter,public :: n0s = 2.e6 ! temperature dependent intercept parameter snow

Referenced by mp_wsm6_init(), mp_wsm6_run(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ n0smax

real(kind=kind_phys), parameter, public mp_wsm6::n0smax = 1.e11
42  real(kind=kind_phys),parameter,public :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)

Referenced by mp_wsm6_run(), refl10cm_wsm6(), slope_snow(), and slope_wsm6().

◆ pacrc

real(kind=kind_phys), save mp_wsm6::pacrc
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pacrg

real(kind=kind_phys), save mp_wsm6::pacrg
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pacrr

real(kind=kind_phys), save mp_wsm6::pacrr
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pacrs

real(kind=kind_phys), save mp_wsm6::pacrs
private

Referenced by mp_wsm6_init().

◆ peaut

real(kind=kind_phys), parameter, private mp_wsm6::peaut = .55
private
22  real(kind=kind_phys),parameter,private:: peaut = .55 ! collection efficiency

Referenced by mp_wsm6_init().

◆ pfrz1

real(kind=kind_phys), parameter, private mp_wsm6::pfrz1 = 100.
private
35  real(kind=kind_phys),parameter,private:: pfrz1 = 100. ! constant in Biggs freezing

Referenced by mp_wsm6_run().

◆ pfrz2

real(kind=kind_phys), parameter, private mp_wsm6::pfrz2 = 0.66
private
36  real(kind=kind_phys),parameter,private:: pfrz2 = 0.66 ! constant in Biggs freezing

Referenced by mp_wsm6_run().

◆ pi

real(kind=kind_phys), save mp_wsm6::pi
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), and refl10cm_wsm6().

◆ pidn0g

real(kind=kind_phys), save mp_wsm6::pidn0g
private

Referenced by mp_wsm6_init(), slope_graup(), and slope_wsm6().

◆ pidn0r

real(kind=kind_phys), save mp_wsm6::pidn0r
private

Referenced by mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ pidn0s

real(kind=kind_phys), save, public mp_wsm6::pidn0s
64  real(kind=kind_phys),public,save:: pidn0s,pidnc

Referenced by mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ pidnc

real(kind=kind_phys), save, public mp_wsm6::pidnc

Referenced by mp_wsm6_init().

◆ precg1

real(kind=kind_phys), save mp_wsm6::precg1
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precg2

real(kind=kind_phys), save mp_wsm6::precg2
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precr1

real(kind=kind_phys), save mp_wsm6::precr1
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precr2

real(kind=kind_phys), save mp_wsm6::precr2
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precs1

real(kind=kind_phys), save mp_wsm6::precs1
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ precs2

real(kind=kind_phys), save mp_wsm6::precs2
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ pvtg

real(kind=kind_phys), save mp_wsm6::pvtg
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), slope_graup(), and slope_wsm6().

◆ pvtr

real(kind=kind_phys), save mp_wsm6::pvtr
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), slope_rain(), and slope_wsm6().

◆ pvts

real(kind=kind_phys), save mp_wsm6::pvts
private

Referenced by mp_wsm6_init(), mp_wsm6_run(), slope_snow(), and slope_wsm6().

◆ qc0

real(kind=kind_phys), save mp_wsm6::qc0
private
46  real(kind=kind_phys),save:: &
47  qc0,qck1, &
48  bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
49  g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
50  bvtr6,g6pbr, &
51  precr1,precr2,roqimax,bvts1, &
52  bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
53  n0g,avtg,bvtg,deng,lamdagmax, & !RAS13.3 - set these in wsm6init
54  g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
55  xlv1,pacrc,pi, &
56  bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
57  g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
58  precg1,precg2,pidn0g, &
59  rslopermax,rslopesmax,rslopegmax, &
60  rsloperbmax,rslopesbmax,rslopegbmax, &
61  rsloper2max,rslopes2max,rslopeg2max, &
62  rsloper3max,rslopes3max,rslopeg3max

Referenced by ERF::derive_diag_profiles_stag(), mp_wsm6_init(), and mp_wsm6_run().

◆ qck1

real(kind=kind_phys), save mp_wsm6::qck1
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ qcrmin

real(kind=kind_phys), parameter, private mp_wsm6::qcrmin = 1.e-9
private
37  real(kind=kind_phys),parameter,private:: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg

Referenced by mp_wsm6_run(), slope_graup(), slope_rain(), slope_snow(), and slope_wsm6().

◆ qs0

real(kind=kind_phys), parameter, private mp_wsm6::qs0 = 6.e-4
private
40  real(kind=kind_phys),parameter,private:: qs0 = 6.e-4 ! threshold amount for aggretion to occur

Referenced by mp_wsm6_run().

◆ r0

real(kind=kind_phys), parameter, private mp_wsm6::r0 = .8e-5
private
21  real(kind=kind_phys),parameter,private:: r0 = .8e-5 ! 8 microm in contrast to 10 micro m

Referenced by ERF::advance_dycore(), ApplySpongeZoneBCsForCC(), ApplySpongeZoneBCsForMom(), make_buoyancy(), make_fast_coeffs(), make_mom_sources(), make_sources(), mp_wsm6_init(), and ParallelFor().

◆ roqimax

real(kind=kind_phys), save mp_wsm6::roqimax
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ rslopeg2max

real(kind=kind_phys), save mp_wsm6::rslopeg2max
private

Referenced by mp_wsm6_init(), slope_graup(), and slope_wsm6().

◆ rslopeg3max

real(kind=kind_phys), save mp_wsm6::rslopeg3max
private

Referenced by mp_wsm6_init(), slope_graup(), and slope_wsm6().

◆ rslopegbmax

real(kind=kind_phys), save mp_wsm6::rslopegbmax
private

Referenced by mp_wsm6_init(), slope_graup(), and slope_wsm6().

◆ rslopegmax

real(kind=kind_phys), save mp_wsm6::rslopegmax
private

Referenced by mp_wsm6_init(), slope_graup(), and slope_wsm6().

◆ rsloper2max

real(kind=kind_phys), save mp_wsm6::rsloper2max
private

Referenced by mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ rsloper3max

real(kind=kind_phys), save mp_wsm6::rsloper3max
private

Referenced by mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ rsloperbmax

real(kind=kind_phys), save mp_wsm6::rsloperbmax
private

Referenced by mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ rslopermax

real(kind=kind_phys), save mp_wsm6::rslopermax
private

Referenced by mp_wsm6_init(), slope_rain(), and slope_wsm6().

◆ rslopes2max

real(kind=kind_phys), save mp_wsm6::rslopes2max
private

Referenced by mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ rslopes3max

real(kind=kind_phys), save mp_wsm6::rslopes3max
private

Referenced by mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ rslopesbmax

real(kind=kind_phys), save mp_wsm6::rslopesbmax
private

Referenced by mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ rslopesmax

real(kind=kind_phys), save mp_wsm6::rslopesmax
private

Referenced by mp_wsm6_init(), slope_snow(), and slope_wsm6().

◆ xlv1

real(kind=kind_phys), save mp_wsm6::xlv1
private

Referenced by mp_wsm6_init(), and mp_wsm6_run().

◆ xmyu

real(kind=kind_phys), parameter, private mp_wsm6::xmyu = 1.718e-5
private
24  real(kind=kind_phys),parameter,private:: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1

Referenced by mp_wsm6_init().

◆ xncr

real(kind=kind_phys), parameter, private mp_wsm6::xncr = 3.e8
private
23  real(kind=kind_phys),parameter,private:: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80

Referenced by mp_wsm6_init(), and mp_wsm6_run().