CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cldwat2m_micro.F
1
3
19
20!---------------------------------------------------------------------------------
21! Purpose:
22! CAM Interface for microphysics
23!---------------------------------------------------------------------------------
24
25#ifdef NEMS_GSM
26 use machine, only : r8 => kind_phys
27 use physcons, gravit => con_g, rair => con_rd, &
28 & rh2o => con_rv, epsilon => con_eps, &
29 & tmelt => con_tice, cpair => con_cp, &
30 & latvap => con_hvap, latice => con_hfus, &
31 & pi => con_pi
32 use wv_saturation, only : estblf, hlatv, tmin, hlatf, rgasv, pcf,&
33 & epsqs, ttrice, vqsatd2,cp, &
35 use funcphys, only : fpvs, fpvsl, fpvsi
36#endif
37
38#ifdef GEOS5
39 use mapl_constantsmod, r8 => mapl_r8
40 use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, &
41 & pcf, cp, epsqs, ttrice, vqsatd2, &
43#endif
44
45#ifdef CAM
46 use shr_kind_mod, only: r8=>shr_kind_r8
47 use spmd_utils, only: masterproc
48 use ppgrid, only: pcols, pver, pverp
49 use physconst, only: gravit, rair, tmelt, cpair, rh2o,
51 use physconst, only: latvap, latice
52 use abortutils, only: endrun
53 use error_function, only: erf,erfc
54 use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, &
55 & pcf, cp, epsqs, ttrice, vqsatd2, &
57 use cam_history, only: addfld, add_default, phys_decomp, &
58 & outfld
59 use cam_logfile, only: iulog
60 use rad_constituents, only: rad_cnst_get_clim_info,
61 & rad_cnst_get_clim_aer_props
62 use phys_control, only: phys_getopts
63 use cldwat2m_macro, only: rhmini, rhmaxi
64#endif
65
66
67 implicit none
68
69 real(r8), parameter :: zero=0.0_r8, one=1.0_r8, two=2.0_r8 &
70 &, three=3.0_r8, four=4.0_r8, five=5.0_r8 &
71 &, half=0.5_r8, oneb3=one/three &
73#ifdef NEMS_GSM
74!
75 integer, parameter :: iulog = 6
76
77 real(r8), parameter :: rhmini = 0.80_r8
78 real(r8), parameter :: rhmaxi = 1.1_r8
79! real(r8), parameter :: r_universal = 6.02214e26*1.38065e-23
80! real(r8), parameter :: r_universal = con_rgas * 1000.0
81 real(r8), parameter :: mwh2o = 18.016
82 real(r8), parameter :: rhoh2o = 1.000e3
83
84 logical :: ip = .true.
85 real(r8) :: tmn = 173.16_r8, tmx = 375.16_r8, trice = 35.00_r8
86#endif
87
88
89#ifdef GEOS5
90!++jtb
91!
92 real(r8), parameter :: rhmini = 0.80_r8
93 real(r8), parameter :: rhmaxi = 1.1_r8
94 integer, parameter :: iulog = 6
95
96
97
98 real(r8), parameter :: gravit = mapl_grav
99 real(r8), parameter :: rair = mapl_rgas
100 real(r8), parameter :: tmelt = mapl_tice
101 real(r8), parameter :: cpair = mapl_cp
102 real(r8), parameter :: rh2o = mapl_rvap
103 real(r8), parameter :: r_universal = mapl_runiv
104 real(r8), parameter :: mwh2o = mapl_h2omw
105 real(r8), parameter :: rhoh2o = mapl_rhowtr
106 real(r8), parameter :: latvap = mapl_alhl
107 real(r8), parameter :: latice = mapl_alhf
108 real(r8), parameter :: epsilon = mapl_vireps
109
110
111 logical :: ip = .true.
112 real(r8) :: tmn = 173.16_r8
113 real(r8) :: tmx = 375.16_r8
114 real(r8) :: trice = 35.00_r8
115#endif
116!--jtb
117
118
119 private
120! save
121
122 logical, public :: liu_in = .false.
123
124
126
127!constants remaped
128 real(r8), private:: g, ginv
129 real(r8), private:: r
130 real(r8), private:: rv
131! real(r8), private:: rr ! not used
132 real(r8), private:: cpp
133 real(r8), private:: rhow, pirhow
134 real(r8), private:: xxlv
135 real(r8), private:: xlf, xlfocp, cpoxlf
136 real(r8), private:: xxls
137
138 real(r8), private:: rhosn, pirhosn
139 real(r8), private:: rhoi, pirhoi
140
141 real(r8), private:: ac,bc,as,bs,ai,bi,ar,br
142 real(r8), private:: ci,di,oneodi
143 real(r8), private:: cs,ds
144 real(r8), private:: cr,dr
145 real(r8), private:: f1s,f2s
146 real(r8), private:: eii
147 real(r8), private:: ecc
148 real(r8), private:: ecr
149 real(r8), private:: f1r,f2r
150 real(r8), private:: dcs, ts_auto_ice
151 real(r8), private:: qsmall
152 real(r8), private:: qvsmall
153 real(r8), private:: bimm,aimm
154 real(r8), private:: rhosu
155 real(r8), private:: mi0
156 real(r8), private:: rin
157 real(r8), private:: qcvar
158! real(r8), private:: pi
159
160! Additional constants to help speed up code
161
162 real(r8), private:: cons1, cons2, cons3, cons4, cons5
163 &, cons6, cons7, cons8, cons9, cons10
164 &, cons11, cons12, cons13, cons14, cons15
165 &, cons16, cons17, cons18, cons19, cons20
166 &, cons21, cons22, cons23, cons24, cons25
167 &, cons27, cons28
168
169 real(r8), private:: lammini, lammaxi, lamminr, lammaxr
170 &, lammins, lammaxs
171
173 real(r8), private, parameter :: tmax_fsnow = tmelt
174 &, tmin_fsnow = tmelt-5._r8
175
176!needed for findsp
177 real(r8), private:: tt0
178
179!switch for specification of droplet and crystal number
180
181 real(r8), private:: csmin,csmax,minrefl,mindbz
182
183
184 contains
185
186!===============================================================================
190 subroutine ini_micro(Dcs_, QCVAR_, ts_auto_ice_)
191
192!-----------------------------------------------------------------------
193!
194! Purpose:
195! initialize constants for the morrison microphysics
196! called from stratiform.F90
197!
198! Author: Andrew Gettelman Dec 2005
199!
200!-----------------------------------------------------------------------
201
202#ifdef CAM
203 use cloud_fraction, only: cldfrc_getparams
204#endif
205 real(r8), intent(in) :: dcs_, qcvar_, ts_auto_ice_
206
207
208 integer k, l, m, iaer
209 real(r8) surften, arg, derf
210
211 character(len=16) :: eddy_scheme = ' '
212 logical :: history_microphysics
213
214
215
216#ifdef CAM
217
218 call phys_getopts(eddy_scheme_out = eddy_scheme,
219 & history_microphysics_out = history_microphysics )
220
221
222 call addfld ('QRAIN ','kg/kg ',pver, 'A','Diagnostic grid-
223 &mean rain mixing ratio' ,phys_decomp)
224 call addfld ('QSNOW ','kg/kg ',pver, 'A','Diagnostic grid-
225 &mean snow mixing ratio' ,phys_decomp)
226 call addfld ('NRAIN ','m-3 ',pver, 'A','Diagnostic grid-
227 &mean rain number conc' ,phys_decomp)
228 call addfld ('NSNOW ','m-3 ',pver, 'A','Diagnostic grid-
229 &mean snow number conc' ,phys_decomp)
230
231
232 call addfld ('RERCLD ','m ',pver, 'A',
233 &'Diagnostic effective radius of Liquid Cloud and Rain' ,
234 &phys_decomp)
235
236 call addfld ('DSNOW ','m ',pver, 'A','Diagnostic grid-
237 &mean snow diameter' ,phys_decomp)
238
239
240 call addfld ('MGFLXPRC ','kg/m2/s ',pver+1, 'A','Diagnostic grid-
241 &mean rain flux at layer interface', phys_decomp)
242 call addfld ('MGFLXSNW ','kg/m2/s ',pver+1, 'A','Diagnostic grid-
243 &mean snow flux at layer interface', phys_decomp)
244
245
246 call addfld ('REFL ','DBz ',pver, 'A',
247 &'94 GHz radar reflectivity' ,phys_decomp)
248 call addfld ('AREFL ','DBz ',pver, 'A',
249 &'Average 94 GHz radar reflectivity' ,phys_decomp)
250 call addfld ('FREFL ','fraction ',pver, 'A',
251 &'Fractional occurance of radar reflectivity' ,phys_decomp)
252
253 call addfld ('CSRFL ','DBz ',pver, 'A',
254 &'94 GHz radar reflectivity (CloudSat thresholds)' ,phys_decomp)
255 call addfld ('ACSRFL ','DBz ',pver, 'A',
256 &'Average 94 GHz radar reflectivity (CloudSat thresholds)' ,
257 &phys_decomp)
258 call addfld ('FCSRFL ','fraction ',pver, 'A',
259 &'Fractional occurance of radar reflectivity (CloudSat thresholds)'
260 & ,phys_decomp)
261
262 call addfld ('AREFLZ ','mm^6/m^3 ',pver, 'A',
263 &'Average 94 GHz radar reflectivity' ,phys_decomp)
264
265
266 call addfld ('NCAL ','#/m3 ',pver, 'A',
267 &'Number Concentation Activated for Liquid',phys_decomp)
268 call addfld ('NCAI ','#/m3 ',pver, 'A',
269 &'Number Concentation Activated for Ice',phys_decomp)
270
271
272
273 call addfld ('AQRAIN ','kg/kg ',pver, 'A',
274 &'Average rain mixing ratio' ,phys_decomp)
275 call addfld ('AQSNOW ','kg/kg ',pver, 'A',
276 &'Average snow mixing ratio' ,phys_decomp)
277 call addfld ('ANRAIN ','m-3 ',pver, 'A',
278 &'Average rain number conc' ,phys_decomp)
279 call addfld ('ANSNOW ','m-3 ',pver, 'A',
280 &'Average snow number conc' ,phys_decomp)
281 call addfld ('ADRAIN ','Micron ',pver, 'A',
282 &'Average rain effective Diameter' ,phys_decomp)
283 call addfld ('ADSNOW ','Micron ',pver, 'A',
284 &'Average snow effective Diameter' ,phys_decomp)
285 call addfld ('FREQR ','fraction ',pver, 'A',
286 &'Fractional occurance of rain' ,phys_decomp)
287 call addfld ('FREQS ','fraction ',pver, 'A',
288 &'Fractional occurance of snow' ,phys_decomp)
289
290 if ( history_microphysics) then
291 call add_default ('AQSNOW ', 1, ' ')
292 call add_default ('FREQR ', 1, ' ')
293 call add_default ('FREQS ', 1, ' ')
294 call add_default ('AQRAIN ', 1, ' ')
295 call add_default ('AQSNOW ', 1, ' ')
296 call add_default ('ANRAIN ', 1, ' ')
297 call add_default ('ANSNOW ', 1, ' ')
298 end if
299
300#endif
301!--jtb commented out when GEOS5
302
303!declarations for morrison codes (transforms variable names)
304
305 g = gravit
306 ginv = one / g
307 r = rair
308 rv = rh2o
309! rr = r_universal
310 cpp = cpair
311 cp = cpair
312 rhow = rhoh2o
313
314! latent heats
315
316 xxlv = latvap
317 xlf = latice
318 xxls = xxlv + xlf
319 xlfocp = xlf / cpair
320 cpoxlf = cpair / xlf
321
322! write(0,*)' xlfocp=',xlfocp,' cpoxlf=',cpoxlf
323! parameters below from Reisner et al. (1998)
324! density parameters (kg/m3)
325
326 rhosn = 100._r8
327 rhoi = 500._r8
328 rhow = 1000._r8
329
330
331! fall speed parameters, V = aD^b
332! V is in m/s
333
334! droplets
335 ac = 3.e7_r8
336 bc = two
337
338! snow
339 as = 11.72_r8
340 bs = 0.41_r8
341
342! cloud ice
343 ai = 700._r8
344 bi = one
345
346! rain
347 ar = 841.99667_r8
348 br = 0.8_r8
349
350! particle mass-diameter relationship
351! currently we assume spherical particles for cloud ice/snow
352! m = cD^d
353
354! pi= 3.1415927_r8
355! pi= 3.1415926535897931_r8
356
357! pi = four*atan(one)
358
359 pirhow = pi * rhow
360 pirhosn = pi * rhosn
361 pirhoi = pi * rhoi
362
363! cloud ice mass-diameter relationship
364
365 ci = pirhoi/6._r8
366 di = three
367 oneodi = one / di
368
369! snow mass-diameter relationship
370
371 cs = pirhosn/6._r8
372 ds = three
373
374! drop mass-diameter relationship
375
376 cr = pirhow/6._r8
377 dr = three
378
379! ventilation parameters for snow
380! hall and prupacher
381
382 f1s = 0.86_r8
383 f2s = 0.28_r8
384
385! collection efficiency, aggregation of cloud ice and snow
386
387 eii = 0.1_r8
388
389! collection efficiency, accretion of cloud water by rain
390
391 ecr = one
392
393! ventilation constants for rain
394
395 f1r = 0.78_r8
396 f2r = 0.32_r8
397
398! autoconversion size threshold for cloud ice to snow (m)
399
400 dcs = dcs_ * 1.0e-6_r8
401! ice autoconversion time scale (default 180s)
402
403 ts_auto_ice = ts_auto_ice_
404
405! smallest mixing ratio considered in microphysics
406
407 qsmall = 1.e-18_r8
408 qvsmall = 1.e-6_r8
409
410! immersion freezing parameters, bigg 1953
411
412 bimm = 100._r8
413 aimm = 0.66_r8
414
415! typical air density at 850 mb
416
417 rhosu = 85000._r8/(rair * tmelt)
418
419! mass of new crystal due to aerosol freezing and growth (kg)
420
421 mi0 = (four/three)*pirhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8)
422
423
424! radius of contact nuclei aerosol (m)
425
426 rin = 0.1e-6_r8
427
428! 1 / relative variance of sub-grid cloud water distribution
429! see morrison and gettelman, 2008, J. Climate for details
430
431
432
433 qcvar = qcvar_
434
435! freezing temperature
436 tt0 = 273.15_r8
437
438
439!++jtb
440 tmn = 173.16_r8
441 tmx = 375.16_r8
442 trice = 35.00_r8
443 ip = .true.
444
445 call gestbl(tmn ,tmx ,trice ,ip ,epsilon , latvap ,latice ,rh2o ,
446 & cpair ,tmelt )
447
448!--jtb (08/18/10)
449
450! pi = 4._r8*atan(1.0_r8)
451
452
453 csmin = -30._r8
454 csmax = 26._r8
455 mindbz = -99._r8
456
457 minrefl = 1.26e-10_r8
458
459! Define constants to help speed up code (limit calls to gamma function)
460
461 cons1 = gamma(one+di)
462! cons1 = gamma(one+miu_ice+ di)/gamma(one+miu_ice)
463
464 cons2 = gamma(qcvar+2.47_r8)
465 cons3 = gamma(qcvar)
466 cons4 = gamma(one+br)
467 cons5 = gamma(four+br)
468 cons6 = gamma(one+ds)
469 cons7 = gamma(one+bs)
470 cons8 = gamma(four+bs)
471 cons9 = gamma(qcvar+two)
472 cons10 = gamma(qcvar+one)
473 cons11 = gamma(three+bs)
474 cons12 = gamma(qcvar+1.15_r8)
475 cons13 = gamma((five+br)*half)
476 cons14 = gamma((five+bs)*half)
477 cons15 = gamma(qcvar+bc/three)
478
479
480
481
482
483 cons18 = qcvar**2.47_r8
484 cons19 = qcvar*qcvar
485 cons20 = qcvar**1.15_r8
486 cons21 = qcvar**(bc/three)
487 cons22 = (four/three)*pirhow*(25.e-6_r8)**3
488 cons23 = dcs*dcs*dcs
489 cons24 = dcs*dcs
490 cons25 = dcs**bs
491 cons27 = xxlv*xxlv
492 cons28 = xxls*xxls
493
494 lammaxi = one / 1.e-6_r8
495 lammini = one / (one*dcs)
496 lammaxr = one / 20.e-6_r8
497 lamminr = one / 500.e-6_r8
498 lammaxs = one / 10.e-6_r8
499 lammins = one / 2000.e-6_r8
500
501 return
502 end subroutine ini_micro
503
504!===============================================================================
508 subroutine mmicro_pcond ( lchnk, ncol, deltatin, tn, ttend, &
509 & pcols, pver, &
510 & qn, qtend, cwtend, qc, qi, nc, ni,fprcp,qrn,qsnw,nrn,nsnw, &
511 & p, pdel, cldn, liqcldf, &
512 & icecldf, cldo, pint, rpdel, zm, rate1ord_cw2pr_st, naai, &
513 & npccnin, rndst,nacon, rhdfda, rhu00, fice, tlat, qvlat, qctend, &
514 & qitend, nctend, nitend, effc, effc_fn, effi, prect, preci, &
515 & nevapr, evapsnow, prain, prodsnow, cmeout, deffi, pgamrad, &
516 & lamcrad,qsout2,dsout2,qrout2, drout2, qcsevap,qisevap,qvres, &
517 & cmeiout, vtrmc,vtrmi,qcsedten,qisedten, prao, &
518 & prco,mnuccco,mnuccto, &
519 & msacwio,psacwso, bergso,bergo,melto,homoo,qcreso,prcio,praio, &
520 & qireso, mnuccro,pracso,meltsdt,frzrdt, ncal, ncai, mnuccdo, &
521 & nnuccto, &
522 & nsout2, nrout2, ncnst, ninst, nimm, miu_disp, nsoot, rnsoot, &
523 & ui_scale, dcrit, nnuccdo, nnuccco, nsacwio, nsubio, nprcio, &
524 & npraio, npccno, npsacwso, nsubco, nprao, nprc1o, tlataux, &
525 & nbincontactdust, lprint, xlat, xlon, rhc)
526! & nbincontactdust, ts_auto_ice,xlat,xlon)
527
528
529!Author: Hugh Morrison, Andrew Gettelman, NCAR
530! e-mail: morrison@ucar.edu, andrew@ucar.edu
531
533
534#ifdef CAM
535 use constituents, only: pcnst
536 real(r8), parameter :: ncnst = 100.0e6
537 real(r8), parameter :: ninst = 0.10e6
538#endif
539
540!--jtb (08/18/10) Is this still needed? Pcnst is not used.
541
542 integer, intent (in) :: pcols, pver,fprcp
543 real(r8), intent (in) :: ncnst
544 real(r8), intent (in) :: ninst
545 real(r8), intent (in) :: nimm (pcols,pver)
546 real(r8), intent (in) :: miu_disp , ui_scale, dcrit
547! real(r8), intent (in) :: miu_disp , ui_scale, dcrit, ts_auto_ice
548 real(r8), intent (in) :: nsoot (pcols,pver) , rnsoot (pcols,pver)
549 & ,xlon,xlat
550! integer, intent(in) :: ktrop_min
551
552
553 logical lprint
554
555
556 integer, intent(in) :: lchnk
557 integer, intent(in) :: ncol
558 real(r8), intent(in) :: deltatin
559 real(r8), intent(in) :: tn(pcols,pver)
560 real(r8), intent(in) :: ttend(pcols,pver)
561 real(r8), intent(in) :: qn(pcols,pver)
562 real(r8), intent(in) :: qtend(pcols,pver)
563 real(r8), intent(in) :: cwtend(pcols,pver)
564
565 real(r8), intent(inout) :: qc(pcols,pver)
566 real(r8), intent(inout) :: qi(pcols,pver)
567 real(r8), intent(inout) :: nc(pcols,pver)
568 real(r8), intent(inout) :: ni(pcols,pver)
569 real(r8), intent(inout) :: qrn(pcols,pver)
570 real(r8), intent(inout) :: qsnw(pcols,pver)
571 real(r8), intent(inout) :: nrn(pcols,pver)
572 real(r8), intent(inout) :: nsnw(pcols,pver)
573 real(r8), intent(in) :: p(pcols,pver)
574 real(r8), intent(in) :: pdel(pcols,pver)
575 real(r8), intent(in) :: cldn(pcols,pver)
576 real(r8), intent(in) :: icecldf(pcols,pver)
577 real(r8), intent(in) :: liqcldf(pcols,pver)
578 real(r8), intent(inout) :: cldo(pcols,pver)
579 real(r8), intent(in) :: pint(pcols,pver+1)
580
581 real(r8), intent(in) :: rpdel(pcols,pver)
582 real(r8), intent(in) :: zm(pcols,pver)
583! real(r8), intent(in) :: omega(pcols,pver)
584 real(r8), intent(in) :: rhc(pcols,pver)
585
586 real(r8), intent(out) :: rate1ord_cw2pr_st(pcols,pver)
587
588! Inputs for aerosol activation
589 real(r8), intent(in) :: naai(pcols,pver)
590 real(r8), intent(inout) :: npccnin(pcols,pver)
591 integer :: nbincontactdust
592 real(r8), intent(in), dimension(pcols,pver, 10) :: rndst, nacon
593
594
595 real(r8), intent(in) :: rhdfda(pcols,pver)
596 real(r8), intent(in) :: rhu00(pcols,pver)
597 real(r8), intent(in) :: fice(pcols,pver)
598
599 real(r8), intent(out) :: tlat(pcols,pver)
600
601 real(r8), intent(out) :: tlataux(pcols,pver)
602
603 real(r8), intent(out) :: qvlat(pcols,pver)
604 real(r8), intent(out) :: qctend(pcols,pver)
605 real(r8), intent(out) :: qitend(pcols,pver)
606 real(r8), intent(out) :: nctend(pcols,pver)
607 real(r8), intent(out) :: nitend(pcols,pver)
608 real(r8), intent(out) :: effc(pcols,pver)
609 real(r8), intent(out) :: effc_fn(pcols,pver)
610 real(r8), intent(out) :: effi(pcols,pver)
611 real(r8), intent(out) :: prect(pcols)
612 real(r8), intent(out) :: preci(pcols)
613 real(r8), intent(out) :: nevapr(pcols,pver)
614 real(r8), intent(out) :: evapsnow(pcols,pver)
615 real(r8), intent(out) :: prain(pcols,pver)
616 real(r8), intent(out) :: prodsnow(pcols,pver)
617 real(r8), intent(out) :: cmeout(pcols,pver)
618 real(r8), intent(out) :: deffi(pcols,pver)
619 real(r8), intent(out) :: pgamrad(pcols,pver)
620 real(r8), intent(out) :: lamcrad(pcols,pver)
621
622 real(r8), intent(out) :: qcsevap(pcols,pver)
623 real(r8), intent(out) :: qisevap(pcols,pver)
624 real(r8), intent(out) :: qvres(pcols,pver)
625 real(r8), intent(out) :: cmeiout(pcols,pver)
626 real(r8), intent(out) :: vtrmc(pcols,pver)
627 real(r8), intent(out) :: vtrmi(pcols,pver)
628 real(r8), intent(out) :: qcsedten(pcols,pver)
629 real(r8), intent(out) :: qisedten(pcols,pver)
630! microphysical process rates for output (mixing ratio tendencies)
631 real(r8), intent(out) :: prao(pcols,pver)
632 real(r8), intent(out) :: prco(pcols,pver)
633 real(r8), intent(out) :: mnuccco(pcols,pver)
634 real(r8), intent(out) :: mnuccto(pcols,pver)
635 real(r8), intent(out) :: msacwio(pcols,pver)
636 real(r8), intent(out) :: psacwso(pcols,pver)
637 real(r8), intent(out) :: bergso(pcols,pver)
638 real(r8), intent(out) :: bergo(pcols,pver)
639 real(r8), intent(out) :: melto(pcols,pver)
640 real(r8), intent(out) :: homoo(pcols,pver)
641 real(r8), intent(out) :: qcreso(pcols,pver)
642 real(r8), intent(out) :: prcio(pcols,pver)
643 real(r8), intent(out) :: praio(pcols,pver)
644 real(r8), intent(out) :: qireso(pcols,pver)
645 real(r8), intent(out) :: mnuccro(pcols,pver)
646 real(r8), intent(out) :: pracso (pcols,pver)
647 real(r8), intent(out) :: meltsdt(pcols,pver)
648 real(r8), intent(out) :: frzrdt (pcols,pver)
649 real(r8), intent(out) :: mnuccdo(pcols,pver)
650 real(r8), intent(out) :: nnuccto(pcols,pver)
651
652 real(r8), intent(out) :: nnuccdo(pcols,pver)
653 real(r8), intent(out) :: nnuccco(pcols,pver)
654 real(r8), intent(out) :: nsacwio(pcols,pver)
655 real(r8), intent(out) :: nsubio(pcols,pver)
656 real(r8), intent(out) :: nprcio(pcols,pver)
657 real(r8), intent(out) :: npraio(pcols,pver)
658
659 real(r8), intent(out) :: npccno(pcols,pver)
660 real(r8), intent(out) :: npsacwso(pcols,pver)
661 real(r8), intent(out) :: nsubco(pcols,pver)
662 real(r8), intent(out) :: nprao(pcols,pver)
663 real(r8), intent(out) :: nprc1o(pcols,pver)
664
665! local workspace
666! all units mks unless otherwise stated
667
668! temporary variables for sub-stepping
669
670 real(r8) :: t1(pcols,pver)
671 real(r8) :: q1(pcols,pver)
672 real(r8) :: qc1(pcols,pver)
673 real(r8) :: qi1(pcols,pver)
674 real(r8) :: nc1(pcols,pver)
675 real(r8) :: ni1(pcols,pver)
676 real(r8) :: tlat1(pcols,pver)
677 real(r8) :: qvlat1(pcols,pver)
678 real(r8) :: qctend1(pcols,pver)
679 real(r8) :: qitend1(pcols,pver)
680 real(r8) :: nctend1(pcols,pver)
681 real(r8) :: nitend1(pcols,pver)
682 real(r8) :: prect1(pcols)
683 real(r8) :: preci1(pcols)
684
685 real(r8) :: tlat1_aux(pcols,pver)
686
687
688! hm 5/12/11
689! temporary variable for old nc before updating with activation tendency
690 real(r8) :: ncold(pcols,pver)
691!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
692
693 real(r8) :: deltat, dti, deltam
694 real(r8) :: omsm
695 real(r8) :: dto2
696 real(r8) :: mincld
697 real(r8), dimension(pcols,pver) :: q, t, rho, irho, rhof, dv
698 &, mu, sc, kap, cldmax, cldm
699 &, icldm, lcldm, cme, cmei
700 &, cwml, cwmi, lcldn, lcldo
701 &, nctend_mixnuc, npccn
702
703 real(r8), dimension(pcols) :: icwc, calpha, cbeta, cbetah
704 &, cgamma, cgamah, rcgama
705 &, cmec1, cmec2, cmec3, cmec4
706
707 real(r8), dimension(pver) :: nnuccd, mnuccd
708 &, qcsinksum_rate1ord
709 &, qcsum_rate1ord
710
711 real(r8) :: qtmp, dum, dum1, dum2, qcld
712 &, arg, alpha
713
714
715 real(r8) :: qcic(pcols,pver)
716 real(r8) :: qiic(pcols,pver)
717 real(r8) :: qniic(pcols,pver)
718 real(r8) :: qric(pcols,pver)
719 real(r8) :: ncic(pcols,pver)
720 real(r8) :: niic(pcols,pver)
721 real(r8) :: nsic(pcols,pver)
722 real(r8) :: nric(pcols,pver)
723 real(r8) :: lami(pver)
724 real(r8) :: n0i(pver)
725 real(r8) :: lamc(pver)
726 real(r8) :: n0c(pver)
727 real(r8) :: lams(pver)
728 real(r8) :: n0s(pver)
729 real(r8) :: lamr(pver)
730 real(r8) :: n0r(pver)
731 real(r8) :: cdist1(pver)
732! combined size of precip & cloud drops
733 real(r8) :: rercld(pcols,pver)
734 real(r8) :: arcld(pcols,pver)
735 real(r8) :: actmp
736 real(r8) :: artmp
737
738 real(r8) :: pgam(pver)
739 real(r8) :: lammax
740 real(r8) :: lammin
741 real(r8) :: nacnt
742 real(r8) :: mnuccc(pver)
743 real(r8) :: nnuccc(pver)
744
745 real(r8) :: mnucct(pver)
746 real(r8) :: nnucct(pver)
747 real(r8) :: msacwi(pver)
748 real(r8) :: nsacwi(pver)
749
750 real(r8) :: prc(pver)
751 real(r8) :: nprc(pver)
752 real(r8) :: nprc1(pver)
753 real(r8) :: nsagg(pver)
754 real(r8) :: dc0
755 real(r8) :: ds0
756 real(r8) :: eci
757 real(r8) :: psacws(pver)
758 real(r8) :: npsacws(pver)
759 real(r8) :: uni
760 real(r8) :: umi
761 real(r8) :: uns(pver)
762 real(r8) :: ums(pver)
763 real(r8) :: unr(pver)
764 real(r8) :: umr(pver)
765 real(r8) :: unc
766 real(r8) :: umc
767 real(r8) :: pracs(pver)
768 real(r8) :: npracs(pver)
769 real(r8) :: mnuccr(pver)
770 real(r8) :: nnuccr(pver)
771 real(r8) :: pra(pver)
772 real(r8) :: npra(pver)
773 real(r8) :: nragg(pver)
774 real(r8) :: prci(pver)
775 real(r8) :: nprci(pver)
776 real(r8) :: prai(pver)
777 real(r8) :: nprai(pver)
778 real(r8) :: qvs
779 real(r8) :: qvi
780 real(r8) :: dqsdt
781 real(r8) :: dqsidt
782 real(r8) :: ab
783 real(r8) :: qclr
784 real(r8) :: abi,oneoabi
785 real(r8) :: epss
786 real(r8) :: epsr
787 real(r8) :: pre(pver)
788 real(r8) :: prds(pver)
789 real(r8) :: qce
790 real(r8) :: qie
791 real(r8) :: nce
792 real(r8) :: nie
793 real(r8) :: ratio
794 real(r8) :: dumc(pcols,pver)
795 real(r8) :: dumnc(pcols,pver)
796 real(r8) :: dumi(pcols,pver)
797 real(r8) :: dumni(pcols,pver)
798 real(r8) :: dums(pcols,pver)
799 real(r8) :: dumns(pcols,pver)
800 real(r8) :: dumr(pcols,pver)
801 real(r8) :: dumnr(pcols,pver)
802! below are parameters for cloud water and cloud ice sedimentation calculations
803 real(r8) :: fr(pver)
804 real(r8) :: fnr(pver)
805 real(r8) :: fc(pver)
806 real(r8) :: fnc(pver)
807 real(r8) :: fi(pver)
808 real(r8) :: fni(pver)
809 real(r8) :: fs(pver)
810 real(r8) :: fns(pver)
811 real(r8) :: faloutr(pver)
812 real(r8) :: faloutnr(pver)
813 real(r8) :: faloutc(pver)
814 real(r8) :: faloutnc(pver)
815 real(r8) :: falouti(pver)
816 real(r8) :: faloutni(pver)
817 real(r8) :: falouts(pver)
818 real(r8) :: faloutns(pver)
819 real(r8) :: faltndr
820 real(r8) :: faltndnr
821 real(r8) :: faltndc
822 real(r8) :: faltndnc
823 real(r8) :: faltndi
824 real(r8) :: faltndni
825 real(r8) :: faltnds
826 real(r8) :: faltndns
827 real(r8) :: faltndqie
828 real(r8) :: faltndqce
829!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
830 real(r8) :: relhum(pcols,pver)
831 real(r8) :: csigma(pcols)
832 real(r8) :: rgvm
833 real(r8) :: arn(pcols,pver)
834 real(r8) :: asn(pcols,pver)
835 real(r8) :: acn(pcols,pver)
836 real(r8) :: ain(pcols,pver)
837 real(r8) :: nsubi(pver)
838 real(r8) :: nsubc(pver)
839 real(r8) :: nsubs(pver)
840 real(r8) :: nsubr(pver)
841 real(r8) :: mtime
842 real(r8) :: dz(pcols,pver)
843
844!fice variable
845 real(r8) :: nfice(pcols,pver)
846
847!add variables for rain and snow flux at layer interfaces
848 real(r8) :: rflx(pcols,pver+1)
849 real(r8) :: sflx(pcols,pver+1)
850
851 real(r8) :: rflx1(pcols,pver+1)
852 real(r8) :: sflx1(pcols,pver+1)
853
854! returns from function/subroutine calls
855 real(r8) :: tsp(pcols,pver)
856 real(r8) :: qsp(pcols,pver)
857 real(r8) :: qsphy(pcols,pver)
858 real(r8) :: qs(pcols)
859 real(r8) :: es(pcols)
860 real(r8) :: esl(pcols,pver)
861 real(r8) :: esi(pcols,pver)
862! real(r8) :: gammas(pcols)
863
864! sum of source/sink terms for diagnostic precip
865
866 real(r8) :: qnitend(pcols,pver)
867 real(r8) :: nstend(pcols,pver)
868 real(r8) :: qrtend(pcols,pver)
869 real(r8) :: nrtend(pcols,pver)
870 real(r8) :: qrtot
871 real(r8) :: nrtot
872 real(r8) :: qstot
873 real(r8) :: nstot
874
875! new terms for Bergeron process
876
877 real(r8) :: dumnnuc
878 real(r8) :: ninew
879 real(r8) :: qinew
880 real(r8) :: qvl
881 real(r8) :: epsi
882 real(r8) :: prd
883 real(r8) :: berg(pcols,pver)
884 real(r8) :: bergs(pver)
885
886!bergeron terms
887 real(r8) :: bergtsf
888 real(r8) :: rhin
889
890! diagnostic rain/snow for output to history
891! values are in-precip (local) !!!!
892
893 real(r8) :: qrout(pcols,pver)
894 real(r8) :: nrout(pcols,pver)
895 real(r8) :: nsout(pcols,pver)
896 real(r8) :: dsout(pcols,pver)
897 real(r8) :: qsout(pcols,pver)
898
899
900!averageed rain/snow for history
901 real(r8) , intent(out) :: qrout2(pcols,pver)
902 real(r8) , intent(out) :: qsout2(pcols,pver)
903 real(r8) , intent(out) :: nrout2(pcols,pver)
904 real(r8) , intent(out) :: nsout2(pcols,pver)
905 real(r8) :: freqs(pcols,pver)
906 real(r8) :: freqr(pcols,pver)
907 real(r8) :: dumfice
908 real(r8), intent(out), dimension(pcols,pver) :: drout2, dsout2
909
910!ice nucleation, droplet activation
911 real(r8) :: dum2i(pcols,pver)
912 real(r8) :: dum2l(pcols,pver)
913 real(r8) :: ncmax(pcols,pver)
914 real(r8) :: nimax
915
916!output fields for number conc
917 real(r8) :: ncai(pcols,pver)
918 real(r8) :: ncal(pcols,pver)
919
920! loop array variables
921 integer i,k,nstep,n, l
922 integer ii,kk, m, ind_aux, km, kp
923
924! loop variables for sub-step solution
925 integer iter,it,ltrue(pcols)
926
927! used in contact freezing via dust particles
928 real(r8) tcnt, viscosity, mfp, nslipsoot, ndfaersoot
929 real(r8), dimension(nbincontactdust) :: ndfaer, nslip, slip
930
931
932! used in ice effective radius
933 real(r8) bbi, cci, ak, iciwc, rvi, riter
934
935! used in Bergeron processe and water vapor deposition
936 real(r8) tk, deles, aprpr, bprpr, cice, qi0, crate, qidep
937
938! mean cloud fraction over the time step
939 real(r8) cldmw(pcols,pver)
940
941! used in secondary ice production
942 real(r8) ni_secp
943
944! variabels to check for RH after rain evap
945
946 real(r8) :: esn
947 real(r8) :: qsn
948 real(r8) :: ttmp
949
950
951 real(r8) :: refl(pcols,pver)
952 real(r8) :: rainrt(pcols,pver)
953 real(r8) :: rainrt1(pcols,pver)
954 real(r8) :: csrfl(pcols,pver)
955 real(r8) :: arefl(pcols,pver)
956 real(r8) :: acsrfl(pcols,pver)
957 real(r8) :: frefl(pcols,pver)
958 real(r8) :: fcsrfl(pcols,pver)
959 real(r8) :: areflz(pcols,pver)
960 real(r8) :: tmp, miu_ice(pver)
961 real(r8) :: cons16
962 real(r8) :: cons17
963
964 real(r8) dmc,ssmc,dstrn
965
966 real(r8), parameter :: cdnl = 0.e6_r8
967
968! integer, parameter :: auto_option=3 ! dcrit only used with auto_option=4
969! integer, parameter :: auto_option=2 ! dcrit only used with auto_option=4
970 integer, parameter :: auto_option=1 ! dcrit only used with auto_option=4
971! integer, parameter :: auto_option=4 ! dcrit only used with auto_option=4
972 real(r8) :: beta6, xs, nssoot, nsdust, taux, psc, bh, vaux, aux,
973 & lw, nw, tx1, tx2, tx3, tx4, tx5, omeps, esloesi
974 &, rdz, rdzi
975
976
977! note: number will be adjusted as needed to keep mean size within bounds,
978! even when cosntant droplet or ice number is used
979
980! ***note: Even if constant cloud ice number is set, ice number is allowed
981! to evolve based on process rates. This is needed in order to calculate
982! the change in mass due to ice nucleation. All other ice microphysical
983! processes are consistent with the specified constant ice number if
984! this switch is turned on.
985
986 logical :: nccons,nicons
987
988
989 omeps = one - epsqs
990
991 nccons = .false.
992 nicons = .false.
993!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
994
995! initialize output fields for number conc qand ice nucleation
996 do k=1,pver
997 miu_ice(k) = zero
998 enddo
999 do k=1,pver
1000 do i=1,ncol
1001 ncai(i,k) = zero
1002 ncal(i,k) = zero
1003
1004!Initialize rain size
1005 rercld(i,k) = zero
1006 arcld(i,k) = zero
1007
1008!initialize radiation output variables
1009 pgamrad(i,k) = zero
1010 lamcrad(i,k) = zero
1011 deffi(i,k) = zero
1012!initialize radiation output variables
1013!initialize water vapor tendency term output
1014 qcsevap(i,k) = zero
1015 qisevap(i,k) = zero
1016 qvres(i,k) = zero
1017 cmeiout(i,k) = zero
1018 vtrmc(i,k) = zero
1019 vtrmi(i,k) = zero
1020 qcsedten(i,k) = zero
1021 qisedten(i,k) = zero
1022
1023 prao(i,k) = zero
1024 prco(i,k) = zero
1025
1026 mnuccco(i,k) = zero
1027 mnuccto(i,k) = zero
1028 msacwio(i,k) = zero
1029 psacwso(i,k) = zero
1030 bergso(i,k) = zero
1031 bergo(i,k) = zero
1032 melto(i,k) = zero
1033 homoo(i,k) = zero
1034 qcreso(i,k) = zero
1035 prcio(i,k) = zero
1036 praio(i,k) = zero
1037 qireso(i,k) = zero
1038 mnuccro(i,k) = zero
1039 pracso(i,k) = zero
1040 meltsdt(i,k) = zero
1041 frzrdt(i,k) = zero
1042 mnuccdo(i,k) = zero
1043 nnuccto(i,k) = zero
1044
1045
1046 nnuccdo(i,k) = zero
1047 nnuccco(i,k) = zero
1048 nsacwio(i,k) = zero
1049 nsubio(i,k) = zero
1050 nprcio(i,k) = zero
1051 npraio(i,k) = zero
1052
1053 npccno(i,k) = zero
1054 npsacwso(i,k) = zero
1055 nsubco(i,k) = zero
1056 nprao(i,k) = zero
1057 nprc1o(i,k) = zero
1058 enddo
1059 enddo
1060
1061! assign variable deltat for sub-stepping...
1062 deltat = deltatin
1063 dti = one / deltat
1064 deltam = one / max(deltat, 150.0_r8)
1065
1066! parameters for scheme
1067
1068! omsm = 0.99999_r8
1069 omsm = 0.99999999_r8
1070 dto2 = 0.5_r8*deltat
1071 mincld = 0.00001_r8
1072
1073! initialize time-varying parameters
1074
1075 do k=1,pver
1076 do i=1,ncol
1077
1078! initialize multi-level fields
1079 t(i,k) = tn(i,k)
1080 q(i,k) = qn(i,k)
1081 rho(i,k) = p(i,k) / (r*t(i,k))
1082 irho(i,k) = one / rho(i,k)
1083 dv(i,k) = 8.794e-5_r8*t(i,k)**1.81_r8/p(i,k)
1084 tx1 = t(i,k) * sqrt(t(i,k)) / (t(i,k)+120._r8)
1085 mu(i,k) = 1.496e-6_r8 * tx1
1086 sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
1087 kap(i,k) = 1.414e3_r8 * 1.496e-6_r8 * tx1
1088
1089! air density adjustment for fallspeed parameters
1090! includes air density correction factor to the
1091! power of 0.54 following Heymsfield and Bansemer 2007
1092
1093 rhof(i,k) = (rhosu*irho(i,k))**0.54
1094
1095 arn(i,k) = ar * rhof(i,k)
1096 asn(i,k) = as * rhof(i,k)
1097 acn(i,k) = ac * rhof(i,k)
1098 ain(i,k) = ai * rhof(i,k)
1099
1100! get dz from dp and hydrostatic approx
1101! keep dz positive (define as layer k-1 - layer k)
1102
1103 dz(i,k) = pdel(i,k) * irho(i,k) * ginv
1104
1105!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1106! droplet activation
1107! hm, modify 5/12/11
1108! get provisional droplet number after activation. This is used for
1109! all microphysical process calculations, for consistency with update of
1110! droplet mass before microphysics
1111
1112! calculate potential for droplet activation if cloud water is present
1113! tendency from activation (npccnin) is read in from companion routine
1114
1115! hm note: npccn is no longer needed below this code - so this can
1116! be rewwritten and this parameters can be removed
1117
1118 lcldm(i,k) = max(liqcldf(i,k),mincld)
1119 if (qc(i,k) >= qsmall) then
1120 npccn(i,k) = max( (lcldm(i,k)*npccnin(i,k)-nc(i,k))
1121 & * deltam, zero)
1122 ncold(i,k) = nc(i,k)
1123 nc(i,k) = nc(i,k) + npccn(i,k)*deltat
1124 else
1125 npccn(i,k) = zero
1126 end if
1127
1128! initialization
1129 t1(i,k) = t(i,k)
1130 q1(i,k) = q(i,k)
1131 qc1(i,k) = qc(i,k)
1132 qi1(i,k) = qi(i,k)
1133 nc1(i,k) = nc(i,k)
1134 ni1(i,k) = ni(i,k)
1135
1136! initialize tendencies to zero
1137 tlat1(i,k) = zero
1138 tlat1_aux(i,k) = zero
1139
1140 qvlat1(i,k) = zero
1141 qctend1(i,k) = zero
1142 qitend1(i,k) = zero
1143 nctend1(i,k) = zero
1144 nitend1(i,k) = zero
1145
1146! initialize precip output
1147 qrout(i,k) = zero
1148 qsout(i,k) = zero
1149 nrout(i,k) = zero
1150 nsout(i,k) = zero
1151 dsout(i,k) = zero
1152
1153! initialize variables for trop_mozart
1154 nevapr(i,k) = zero
1155 evapsnow(i,k) = zero
1156 prain(i,k) = zero
1157 prodsnow(i,k) = zero
1158 cmeout(i,k) = zero
1159
1160! for refl calc
1161 rainrt1(i,k) = zero
1162
1163! initialize precip fraction and output tendencies
1164 cldmax(i,k) = mincld
1165
1166!initialize aerosol number
1167! naer2(i,k,:) = zero
1168 dum2l(i,k) = zero
1169 dum2i(i,k) = zero
1170 ncmax(i,k) = zero
1171
1172! for debug purpose
1173! prect(1:ncol)=0._r8
1174! preci(1:ncol)=0._r8
1175! tlat(1:ncol,1:pver)=0._r8
1176! qvlat(1:ncol,1:pver)=0._r8
1177! qctend(1:ncol,1:pver)=0._r8
1178! qitend(1:ncol,1:pver)=0._r8
1179! nctend(1:ncol,1:pver)=0._r8
1180
1181 enddo
1182 enddo
1183
1184! initialize avg precip rate
1185 do i=1,ncol
1186 prect1(i) = zero
1187 preci1(i) = zero
1188 end do
1189!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1190!Get humidity and saturation vapor pressures big loop1
1191
1192 do k=1,pver ! big loop1 - k loop
1193
1194! find wet bulk temperature and saturation value for provisional t and q without
1195! condensation
1196
1197! call vqsatd_water(t(1,k),p(1,k),es,qs,gammas,ncol)
1198
1199 do i=1,ncol ! big i loop1
1200
1201#ifdef NEMS_GSM
1202 esl(i,k) = min(fpvsl(t(i,k)), p(i,k))
1203 esi(i,k) = min(fpvsi(t(i,k)), p(i,k))
1204#else
1205 esl(i,k) = min(polysvp(t(i,k),0), p(i,k))
1206 esi(i,k) = min(polysvp(t(i,k),1), p(i,k))
1207#endif
1208
1209 esloesi = esl(i,k) / esi(i,k)
1210! hm fix, make sure when above freezing that esi=esl, not active yet
1211 if (t(i,k) > tmelt) esi(i,k) = esl(i,k)
1212
1213 qs(i) = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
1214
1215 relhum(i,k) = min(q(i,k)/qs(i), one) !Anning limiting relhum
1216
1217! if (lprint .and. k==29) write(0,*)' esl=',esl(i,k)
1218! &,' esi=',esi(i,k),' pres=',p(i,k),' t=',t(i,k)
1219! &,' relhum=',relhum(i,k),' q=',q(i,k),' qs=',qs(i)
1220
1221! get cloud fraction, check for minimum
1222
1223 cldm(i,k) = max(cldn(i,k), mincld)
1224 cldmw(i,k) = max(cldn(i,k), mincld)
1225
1226 icldm(i,k) = max(icecldf(i,k), mincld)
1227 lcldm(i,k) = max(liqcldf(i,k), mincld)
1228
1229
1230 if (qc(i,k) >= qsmall) then
1231 tx1 = one / lcldm(i,k)
1232 dum2l(i,k) = (ncold(i,k)+npccn(i,k)*deltat) * tx1
1233 dum2l(i,k) = max(dum2l(i,k),cdnl*irho(i,k))
1234 ncmax(i,k) = max(dum2l(i,k)*lcldm(i,k), zero)
1235 dum2l(i,k) = npccn(i,k)*deltat*tx1
1236 else
1237 dum2l(i,k) = zero
1238 end if
1239
1240
1241
1242! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet)
1243
1244 nfice(i,k) = zero
1245 dumfice = qc(i,k) + qi(i,k)
1246 if (dumfice > qsmall .and. qi(i,k) > qsmall) then
1247 nfice(i,k) = qi(i,k)/dumfice
1248 endif
1249
1250 if (t(i,k) < tmelt - five) then
1251 if (liu_in) then
1252
1253! if aerosols interact with ice set number of activated ice nuclei
1254 dum2 = naai(i,k)
1255
1256 else
1257! cooper curve (factor of 1000 is to convert from L-1 to m-3)
1258 dum2 = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8
1259! put limit on number of nucleated crystals, set to number at T=-30 C
1260! cooper (limit to value at -35 C)
1261 dum2 = min(dum2, 208.9e3_r8)*irho(i,k)
1262 endif
1263
1264 dumnnuc = max((dum2*icldm(i,k)-ni(i,k))*deltam, zero)
1265
1266
1267! get provisional ni and qi after nucleation in order to calculate
1268! Bergeron process below
1269 ninew = ni(i,k) + dumnnuc*deltat
1270 qinew = qi(i,k) + dumnnuc*deltat*mi0
1271!T>268
1272 else
1273 ninew = ni(i,k)
1274 qinew = qi(i,k)
1275 end if
1276
1277! Initialize CME components
1278
1279 cme(i,k) = zero
1280 cmei(i,k) = zero
1281
1282
1283!-------------------------------------------------------------------
1284!Bergeron process
1285! make sure to initialize bergeron process to zero
1286 berg(i,k) = zero
1287 prd = zero
1288
1289!condensation loop.
1290
1291! get in-cloud qi and ni after nucleation
1292 if (icldm(i,k) > zero) then
1293 tx1 = one / icldm(i,k)
1294 qiic(i,k) = qinew * tx1
1295 niic(i,k) = ninew * tx1
1296 else
1297 qiic(i,k) = zero
1298 niic(i,k) = zero
1299 endif
1300
1301
1302 if (nicons) then
1303 niic(i,k) = ninst*irho(i,k)
1304 end if
1305
1306
1307
1308!if T < 0 C then bergeron.
1309 if (t(i,k) < 273.15) then
1310!if ice exists
1311 if (qi(i,k) > qsmall) then
1312
1313 bergtsf = zero
1314
1315 qvi = epsqs*esi(i,k) / (p(i,k)-omeps*esi(i,k))
1316 qvl = epsqs*esl(i,k) / (p(i,k)-omeps*esl(i,k))
1317
1318 dqsidt = xxls*qvi / (rv*t(i,k)*t(i,k))
1319 abi = one + dqsidt*xxls/cpp
1320 oneoabi = one / abi
1321
1322! get ice size distribution parameters
1323
1324 if (qiic(i,k) >= qsmall) then
1325
1326 lami(k) = (cons1*ci*niic(i,k)/qiic(i,k))**oneodi
1327
1328! miu_ice = mui_hemp_l(lami(k))
1329 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
1330 & 10.0_r8), 0.1_r8)
1331 tx1 = one + miu_ice(k)
1332 tx2 = one / gamma(tx1)
1333 aux = (gamma(tx1+di)*tx2) ** oneodi
1334 lami(k) = aux*lami(k)
1335 n0i(k) = niic(i,k) * lami(k)**tx1 * tx2
1336
1337
1338! check for slope
1339! adjust vars
1340 if (lami(k) < lammini*aux) then
1341 miu_ice(k) = zero
1342 lami(k) = lammini
1343 n0i(k) = lami(k)**(di+one)*qiic(i,k)/(ci*cons1)
1344 niic(i,k) = n0i(k)/lami(k)
1345 else if (lami(k) > lammaxi*aux) then
1346 miu_ice(k) = zero
1347 lami(k) = lammaxi
1348 n0i(k) = lami(k)**(di+one)*qiic(i,k)/(ci*cons1)
1349 niic(i,k) = n0i(k)/lami(k)
1350 end if
1351
1352 epsi = (miu_ice(k)+two)*(miu_ice(k)+one)*pi*
1353 & niic(i,k)*rho(i,k)*dv(i,k)/lami(k)
1354
1355!if liquid exists
1356 if (qc(i,k) > qsmall) then
1357
1358!begin bergeron process
1359! do bergeron (vapor deposition with RHw=1)
1360! code to find berg (a rate) goes here
1361
1362! calculate Bergeron process
1363
1364 prd = epsi*(qvl-qvi)*oneoabi
1365
1366 else
1367 prd = zero
1368 end if
1369
1370! multiply by cloud fraction and transfer of existing cloud liquid to ice
1371
1372 berg(i,k) = max(zero, prd*min(icldm(i,k),lcldm(i,k)))
1373
1374 end if ! qiic(i,k) >= qsmall
1375
1376
1377 if (berg(i,k) > zero) then
1378 tx1 = qc(i,k) * dti
1379 bergtsf = max(zero, tx1/berg(i,k))
1380 if(bergtsf < one) berg(i,k) = max(zero, tx1)
1381 endif
1382
1383!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1384
1385 if (bergtsf < one .or. icldm(i,k) > lcldm(i,k)) then
1386
1387 if (qiic(i,k) >= qsmall) then
1388
1389! first case is for case when liquid water is present, but is completely depleted in time step, i.e., bergrsf > 0 but < 1
1390
1391 if (qc(i,k) >= qsmall) then
1392 rhin = (one + relhum(i,k)) * half
1393 if (rhin*esloesi > one) then
1394 prd = epsi*(rhin*qvl-qvi)*oneoabi
1395
1396! multiply by cloud fraction assuming liquid/ice maximum overlap and add to cmei
1397 cmei(i,k) = cmei(i,k)
1398 & + prd * min(icldm(i,k),lcldm(i,k))
1399 & * (one- bergtsf)
1400! if (lprint .and. k == 29) write(0,*)' cmei1=',cmei(i,k),
1401! &' prd=',prd,' bergtsf=',bergtsf
1402
1403 end if
1404 end if
1405
1406! moved here by Moorthi
1407! note: for case of no liquid, need to set liquid cloud fraction to zero
1408! store liquid cloud fraction in 'dum'
1409 if (qc(i,k) < qsmall) then
1410 dum = zero
1411 else
1412 dum = lcldm(i,k)
1413 end if
1414
1415! second case is for pure ice cloud, either no liquid, or icldm > lcldm
1416 if (qc(i,k) < qsmall) then
1417
1418! note: for case of no liquid, need to set liquid cloud fraction to zero
1419! store liquid cloud fraction in 'dum'
1420
1421!Moorthi if (qc(i,k) < qsmall) then
1422! dum = 0._r8
1423! else
1424! dum = lcldm(i,k)
1425! end if
1426
1427! set RH to grid-mean value for pure ice cloud
1428 rhin = relhum(i,k)
1429
1430 if (rhin*esloesi > one) then
1431
1432 prd = epsi*(rhin*qvl-qvi)*oneoabi
1433
1434! multiply by relevant cloud fraction for pure ice cloud
1435! assuming maximum overlap of liquid/ice
1436 cmei(i,k) = cmei(i,k)
1437 & + prd * max((icldm(i,k)-dum), zero)
1438! if (lprint .and. k == 29) write(0,*)' cmei2=',cmei(i,k),
1439! &' prd=',prd,' icldm=',icldm(i,k),' dum=',dum
1440
1441 end if
1442 end if
1443 end if
1444 end if
1445
1446! if deposition, it should not reduce grid mean rhi below 1.0
1447 if(cmei(i,k) > zero .and. relhum(i,k)*esloesi > one)
1448 & cmei(i,k) = min(cmei(i,k),(q(i,k)-qs(i)/esloesi)
1449 & * oneoabi * dti)
1450
1451! if (lprint .and. k == 29) write(0,*)' cmei3=',cmei(i,k)
1452
1453 end if ! if (qi(i,k) > qsmall)
1454
1455
1456!-------------------------------------------------------------------
1457 end if
1458!..............................................................
1459
1460! evaporation should not exceed available water
1461
1462 tx1 = qc(i,k)*dti
1463 if (-berg(i,k) < -tx1) berg(i,k) = max(tx1, zero) ! is this correct?????
1464
1465! berg(i,k) = min(berg(i,k), qc(i,k)/deltat) ! Moorthi - ask anning
1466
1467!sublimation process...
1468
1469 if (relhum(i,k)*esloesi < one .and.
1470 & qiic(i,k) >= qsmall ) then
1471
1472 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
1473 qvl = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
1474 dqsidt = xxls*qvi/(rv*t(i,k)*t(i,k))
1475 abi = one + dqsidt*xxls/cpp
1476 oneoabi = one / abi
1477
1478! get ice size distribution parameters
1479
1480 lami(k) = (cons1*ci*niic(i,k)/qiic(i,k))**oneodi
1481
1482! miu_ice(k)=mui_hemp_l(lami(k))
1483 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
1484 & 10.0_r8), 0.1_r8)
1485 tx1 = one + miu_ice(k)
1486 tx2 = one / gamma(tx1)
1487 aux = (gamma(tx1+di) * tx2) ** oneodi
1488 lami(k) = aux*lami(k)
1489
1490 n0i(k) = niic(i,k) * lami(k)**tx1 * tx2
1491! check for slope
1492! adjust vars
1493 if (lami(k) < lammini*aux) then
1494 miu_ice(k) = zero
1495 lami(k) = lammini
1496 n0i(k) = lami(k)**(di+one)*qiic(i,k)/(ci*cons1)
1497 else if (lami(k) > lammaxi*aux) then
1498 miu_ice(k) = zero
1499 lami(k) = lammaxi
1500 n0i(k) = lami(k)**(di+one)*qiic(i,k)/(ci*cons1)
1501 end if
1502
1503 epsi = (miu_ice(k)+two)*(miu_ice(k)+one)*pi*
1504 & niic(i,k)*rho(i,k)*dv(i,k) / lami(k)
1505
1506! modify for ice fraction below
1507 prd = epsi*(relhum(i,k)*qvl-qvi)*oneoabi * icldm(i,k)
1508 cmei(i,k) = min(prd, zero)
1509
1510! if (lprint .and. k == 29) write(0,*)' cmei3a=',cmei(i,k)
1511! &,' prd=',prd
1512 endif
1513
1514! sublimation should not exceed available ice
1515
1516 cmei(i,k) = max(cmei(i,k), -qi(i,k)*dti)
1517
1518! if (lprint .and. k == 29) write(0,*)' cmei3b=',cmei(i,k)
1519! &,' qi=',qi(i,k),' deltat=',deltat
1520
1521! sublimation should not increase grid mean rhi above 1.0
1522 if(cmei(i,k) < zero .and. relhum(i,k)*esloesi < one)
1523 & cmei(i,k) = min(zero, max(cmei(i,k),(q(i,k)-qs(i)
1524 & /esloesi) * oneoabi * dti ))
1525! if (lprint .and. k == 29) write(0,*)' cmei3c=',cmei(i,k)
1526! &,' q=',q(i,k),' qs=',qs(i),' esi=',esi(i,k)
1527! &,' esl=',esl(i,k),' abi=',abi
1528
1529! limit cmei due for roundoff error
1530
1531 cmei(i,k) = cmei(i,k)*omsm
1532! if (lprint .and. k == 29) write(0,*)' cmei4=',cmei(i,k),
1533! &' omsm=',omsm
1534
1535! conditional for ice nucleation
1536
1537 if (t(i,k) < (tmelt - five)) then
1538 if ( liu_in ) then
1539
1540! using Liu et al. (2007)/Barahona & Nenes (2009) ice nucleation with hooks into simulated aerosol
1541! ice nucleation rate (dum2) has already been calculated and read in (naai)
1542
1543 dum2i(i,k) = naai(i,k)
1544 else
1545
1546! cooper curve (factor of 1000 is to convert from L-1 to m-3)
1547 dum2i(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))
1548 & * 1000._r8
1549! put limit on number of nucleated crystals, set to number at T=-30 C
1550! cooper (limit to value at -35 C)
1551 dum2i(i,k) = min(dum2i(i,k),208.9e3_r8)*irho(i,k)
1552 endif
1553 else
1554 dum2i(i,k) = zero
1555 end if
1556
1557 end do ! end big i loop1
1558 end do !end big loop1 - k loop
1559
1560!!
1561 do i=1,ncol
1562 rflx1(i,1) = zero ! initialize sub-step precip flux variables
1563 sflx1(i,1) = zero
1564
1565 rflx(i,1) = zero ! initialize final precip flux variables.
1566 sflx(i,1) = zero
1567 ltrue(i) = 0
1568 end do
1569 do k=1,pver
1570 do i=1,ncol
1571 cldo(i,k) = cldn(i,k)
1572
1573 rflx1(i,k+1) = zero ! initialize sub-step precip flux variables
1574 sflx1(i,k+1) = zero
1575
1576 rflx(i,k+1) = zero ! initialize final precip flux variables.
1577 sflx(i,k+1) = zero
1578
1579! skip microphysical calculations if no cloud water
1580
1581 if ((qc(i,k) >= qsmall .or. qi(i,k) >= qsmall .or.
1582 & cmei(i,k) >= qsmall).and.q(i,k)>=qvsmall) ltrue(i) = 1
1583
1584 rate1ord_cw2pr_st(i,k) = zero
1585 end do
1586 end do
1587!!
1588
1589
1590! assign number of sub-steps to iter
1591! use 2 sub-steps, following tests described in MG2008
1592! Anning Cheng 9/17/2016
1593 if (fprcp == 1) then
1594 iter = 1
1595 else
1596 iter = 2
1597 end if
1598
1599 riter = one / float(iter)
1600! get sub-step time step
1601 deltat = deltat * riter
1602 dti = one / deltat
1603
1604! since activation/nucleation processes are fast, need to take into account
1605! factor mtime = mixing timescale in cloud / model time step
1606! mixing time can be interpreted as cloud depth divided by sub-grid vertical velocity
1607! for now mixing timescale is assumed to be 1 timestep for modal aerosols, 20 min bulk
1608
1609! note: mtime for bulk aerosols was set to: mtime=deltat/1200._r8
1610
1611 mtime = one
1612
1613!!!! skip calculations if no cloud water
1614
1615 do i=1,ncol !big i loop2
1616 if (ltrue(i) == 0) then
1617 prect(i) = zero
1618 preci(i) = zero
1619 do k=1,pver
1620 tlat(i,k) = zero
1621 qvlat(i,k) = zero
1622 qctend(i,k) = zero
1623 qitend(i,k) = zero
1624 qnitend(i,k) = zero
1625 qrtend(i,k) = zero
1626 nctend(i,k) = zero
1627 nitend(i,k) = zero
1628 nrtend(i,k) = zero
1629 nstend(i,k) = zero
1630 qniic(i,k) = zero
1631 qric(i,k) = zero
1632 nsic(i,k) = zero
1633 nric(i,k) = zero
1634 rainrt(i,k) = zero
1635 enddo
1636
1637 else
1638
1639 do k=1,pver
1640 qcsinksum_rate1ord(k) = zero
1641 qcsum_rate1ord(k) = zero
1642 enddo
1643
1644
1645!!!!!!!!! begin sub-step!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1646!.....................................................................................................
1647 do it=1,iter ! big iter loop
1648
1649! initialize sub-step microphysical tendencies
1650
1651 do k=1,pver
1652 tlat(i,k) = zero
1653 qvlat(i,k) = zero
1654 qctend(i,k) = zero
1655 qitend(i,k) = zero
1656 qnitend(i,k) = zero
1657 qrtend(i,k) = zero
1658 nctend(i,k) = zero
1659 nitend(i,k) = zero
1660 nrtend(i,k) = zero
1661 nstend(i,k) = zero
1662
1663! initialize diagnostic precipitation to zero
1664
1665 qniic(i,k) = zero
1666 qric(i,k) = zero
1667 nsic(i,k) = zero
1668 nric(i,k) = zero
1669
1670 rainrt(i,k) = zero
1671
1672 enddo
1673
1674
1675! begin new i,k loop, calculate new cldmax after adjustment to cldm above
1676
1677! initialize vertically-integrated rain and snow tendencies
1678
1679 qrtot = zero
1680 nrtot = zero
1681 qstot = zero
1682 nstot = zero
1683
1684! initialize precip at surface
1685
1686 prect(i) = zero
1687 preci(i) = zero
1688
1689 do k=1,pver
1690
1691 km = k - 1
1692! set cwml and cwmi to current qc and qi
1693
1694 cwml(i,k) = qc(i,k)
1695 cwmi(i,k) = qi(i,k)
1696
1697! initialize precip fallspeeds to zero
1698
1699 ums(k) = zero
1700 uns(k) = zero
1701 umr(k) = zero
1702 unr(k) = zero
1703
1704! calculate precip fraction based on maximum overlap assumption
1705
1706 if (k == 1) then
1707 cldmax(i,k) = cldm(i,k)
1708 else
1709! if rain or snow mix ratio is smaller than
1710! threshold, then set cldmax to cloud fraction at current level
1711 if (qric(i,km) >= qsmall .or.
1712 & qniic(i,km) >= qsmall) then
1713 cldmax(i,k) = max(cldmax(i,km),cldm(i,k))
1714 else
1715 cldmax(i,k) = cldm(i,k)
1716 end if
1717 end if
1718
1719 rdz = rho(i,k) * dz(i,k)
1720 rdzi = one / rdz
1721
1722!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1723!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1724
1725! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
1726
1727 if (dum2i(i,k) > zero .and. t(i,k) < (tmelt-five) .and.
1728 & relhum(i,k)*esl(i,k)/esi(i,k) > rhmini+0.05_r8) then
1729
1730!if NCAI > 0. then set numice = ncai (as before)
1731!note: this is gridbox averaged
1732
1733 nimax = dum2i(i,k)*icldm(i,k)
1734! nnuccd(k) = (dum2i(i,k)-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
1735 nnuccd(k) = max(zero, (nimax-ni(i,k))*dti)
1736
1737!Calc mass of new particles using new crystal mass...
1738!also this will be multiplied by mtime as nnuccd is...
1739
1740 mnuccd(k) = nnuccd(k) * mi0
1741
1742! add mnuccd to cmei....
1743 cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime
1744
1745! if (lprint .and. k == 29) write(0,*)' cmei5=',cmei(i,k)
1746
1747! limit cmei
1748
1749 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
1750 dqsidt = xxls*qvi/(rv*t(i,k)*t(i,k))
1751 abi = one + dqsidt*xxls/cpp
1752 cmei(i,k) = min(cmei(i,k),(q(i,k)-qvi)/(abi*deltat))
1753
1754! limit for roundoff error
1755 cmei(i,k) = cmei(i,k)*omsm
1756
1757! if (lprint .and. k == 29) write(0,*)' cmei6=',cmei(i,k)
1758 else
1759 nnuccd(k) = zero
1760 nimax = zero
1761 mnuccd(k) = zero
1762 end if
1763
1764
1765!c............................................................................
1766!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1767! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
1768! for microphysical process calculations
1769! units are kg/kg for mixing ratio, 1/kg for number conc
1770
1771! limit in-cloud values to 0.005 kg/kg
1772
1773 tx1 = one / lcldm(i,k)
1774 tx2 = one / icldm(i,k)
1775 qcic(i,k) = min(cwml(i,k)*tx1, 5.e-3_r8)
1776 qiic(i,k) = min(cwmi(i,k)*tx2, 5.e-3_r8)
1777 ncic(i,k) = max(nc(i,k)*tx1, zero)
1778 niic(i,k) = max(ni(i,k)*tx2, zero)
1779
1780
1781 if (nccons) then
1782 ncic(i,k) = ncnst*irho(i,k)
1783 end if
1784
1785 if (nicons) then
1786 niic(i,k) = ninst*irho(i,k)
1787 end if
1788
1789 tx1 = qc(i,k) - berg(i,k)*deltat
1790 if (tx1 < qsmall) then
1791 qcic(i,k) = zero
1792 ncic(i,k) = zero
1793 if (tx1 < zero) then
1794 berg(i,k) = qc(i,k)*dti*omsm
1795 end if
1796 end if
1797
1798 tx1 = qi(i,k) + (cmei(i,k)+berg(i,k))*deltat
1799 if (tx1 < qsmall) then
1800 qiic(i,k) = zero
1801 niic(i,k) = zero
1802 if (tx1 < zero) then
1803 cmei(i,k) = (-qi(i,k)*dti-berg(i,k))*omsm
1804! if (lprint .and. k == 29) write(0,*)' cmei7=',cmei(i,k)
1805 end if
1806 end if
1807
1808! add to cme output
1809
1810 cmeout(i,k) = cmeout(i,k) + cmei(i,k)
1811
1812! decrease in number concentration due to sublimation/evap
1813! divide by cloud fraction to get in-cloud decrease
1814! don't reduce Nc due to bergeron process
1815
1816!Moved it here since nsubi since cmei was not limited before(DONIF 03/13/2015)
1817 if (cmei(i,k) < zero .and. qi(i,k) > qsmall
1818 & .and. cldm(i,k) > mincld) then
1819 nsubi(k) = cmei(i,k)*ni(i,k) / (qi(i,k)*cldm(i,k))
1820 else
1821 nsubi(k) = zero
1822 end if
1823
1824 nsubc(k) = zero
1825
1826
1827
1828
1829!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1830! get size distribution parameters based on in-cloud cloud water/ice
1831! these calculations also ensure consistency between number and mixing ratio
1832!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1833
1834!......................................................................
1835! cloud ice
1836
1837 if (qiic(i,k) >= qsmall) then
1838
1839! add upper limit to in-cloud number concentration to prevent numerical error
1840 niic(i,k) = min(niic(i,k),qiic(i,k)*1.e20_r8)
1841
1842
1843 lami(k) = (cons1*ci*niic(i,k)/qiic(i,k))**oneodi
1844
1845! miu_ice(k) = mui_hemp_l(lami(k))
1846 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
1847 & 10.0_r8), 0.1_r8)
1848 tx1 = one + miu_ice(k)
1849 tx2 = one / gamma(tx1)
1850 aux = (gamma(tx1+di) * tx2) **oneodi
1851 lami(k) = aux*lami(k)
1852
1853 n0i(k) = niic(i,k)*lami(k)**tx1 * tx2
1854
1855! check for slope
1856! adjust vars
1857
1858 if (lami(k) < lammini) then
1859 lami(k) = lammini
1860 n0i(k) = lami(k)**(di+one)*qiic(i,k)/(ci*cons1)
1861 niic(i,k) = n0i(k)/lami(k)
1862 else if (lami(k) > lammaxi) then
1863 lami(k) = lammaxi
1864 n0i(k) = lami(k)**(di+one)*qiic(i,k)/(ci*cons1)
1865 niic(i,k) = n0i(k)/lami(k)
1866 end if
1867
1868 else
1869 lami(k) = zero
1870 n0i(k) = zero
1871 end if
1872
1873 if (qcic(i,k) >= qsmall) then
1874
1875! add upper limit to in-cloud number concentration to prevent numerical error
1876 ncic(i,k) = min(ncic(i,k),qcic(i,k)*1.e20_r8)
1877
1878 ncic(i,k) = max(ncic(i,k),cdnl*irho(i,k))
1879
1880! get pgam from fit to observations of martin et al. 1994
1881
1882 pgam(k) = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))
1883 & + 0.2714_r8
1884
1885
1886 if (.true.) then
1887 if ((ncic(i,k) > 1.0e-3) .and.
1888 & (qcic(i,k) > 1.0e-11)) then
1889 xs = 0.07_r8*(1000._r8*qcic(i,k)/ncic(i,k))
1890 & **(-0.14_r8)
1891 else
1892 xs = 1.2
1893 end if
1894
1895 xs = max(min(xs, 1.7_r8), 1.1_r8)
1896 xs = xs*xs*xs
1897 xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - four)/8.0_r8
1898 pgam(k) = sqrt(xs)
1899
1900 end if
1901
1902
1903
1904 pgam(k) = one / (pgam(k)*pgam(k)) - one
1905 pgam(k) = max(two, min(15._r8, pgam(k)))
1906
1907
1908! calculate lamc
1909 tx1 = pirhow * gamma(pgam(k)+four)
1910 tx3 = gamma(pgam(k)+one)
1911 tx2 = 6._r8 * qcic(i,k) * tx3
1912
1913 lamc(k) = (tx1*ncic(i,k)/tx2) ** oneb3
1914
1915! lammin, 50 micron diameter max mean size
1916
1917 lammin = (pgam(k)+one) / 50.e-6_r8
1918 lammax = (pgam(k)+one) / 2.e-6_r8
1919
1920 if (lamc(k) < lammin) then
1921 lamc(k) = lammin
1922 ncic(i,k) = tx2*lamc(k)*lamc(k)*lamc(k) / tx1
1923 else if (lamc(k) > lammax) then
1924 lamc(k) = lammax
1925 ncic(i,k) = tx2*lamc(k)*lamc(k)*lamc(k) / tx1
1926 end if
1927
1928! parameter to calculate droplet freezing
1929
1930 cdist1(k) = ncic(i,k) / tx3
1931
1932 else
1933 lamc(k) = zero
1934 cdist1(k) = zero
1935 end if
1936
1937!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1938!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1939! begin micropysical process calculations
1940!.................................................................
1941! autoconversion of cloud liquid water to rain
1942! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
1943! minimum qc of 1 x 10^-8 prevents floating point error
1944
1945 xs = 0.0
1946
1947 if (qcic(i,k) >= 1.e-8_r8) then
1948
1949! nprc is increase in rain number conc due to autoconversion
1950! nprc1 is decrease in cloud droplet conc due to autoconversion
1951
1952! assume exponential sub-grid distribution of qc, resulting in additional
1953! factor related to qcvar below
1954
1955
1956! prc(k) = cons2/(cons3*cons18)*1350._r8
1957! & * qcic(i,k)**2.47_r8
1958! & * (1.e-6_r8*ncic(i,k)*rho(i,k))**(-1.79_r8)
1959
1960 if (auto_option == 1) then
1961
1962 tx1 = qcic(i,k)*rho(i, k)/3.0e-4
1963 prc(k) = 1.0e-3 * qcic(i,k) * (one-exp(-tx1*tx1))
1964 & * gamma(one+qcvar)/(cons3*qcvar)
1965
1966 elseif (auto_option == 2) then
1967
1968 tx1 = qcic(i,k)*rho(i, k)/7.5e-4
1969 tx1 = tx1 * tx1
1970 prc(k) = 1.0e-4 * qcic(i,k) * (one-exp(-tx1*tx1))
1971 & * gamma(one+qcvar)/(cons3*qcvar)
1972
1973 elseif (auto_option == 3) then
1974
1975 tx1 = qcic(i,k) / 3.0e-4
1976
1977 prc(k) = 1.0e-3 * (one-exp(-tx1*tx1))
1978 & * cons2/(cons3*cons18)
1979 & / (1.e-6_r8/50.0*ncic(i,k)*rho(i,k))**1.79_r8
1980
1981 elseif (auto_option == 4) then
1982
1983 xs = one / (one+pgam(k))
1984
1985 beta6 = (one+3.0*xs)*(one+4.0*xs)*(one+5.0*xs)
1986 & / ((one+xs)*(one+xs+xs))
1987
1988 lw = 1.0e-3_r8 * qcic(i,k) * rho(i,k)
1989 nw = ncic(i,k) * rho(i,k) * 1.e-6_r8
1990
1991 xs = min(20.0, 1.03e16*(lw*lw)/(nw*sqrt(nw)))
1992 prc(k) = 1.1e10*beta6*lw*lw*lw
1993 & * (one-exp(-(xs**miu_disp))) / nw
1994 prc(k) = prc(k)*1.0e3*irho(i,k)
1995 prc(k) = prc(k) * gamma(two+qcvar)
1996 & / (gamma(qcvar)*(qcvar*qcvar))
1997
1998 prc(k) = prc(k)*dcrit
1999
2000 xs = 1/xs
2001 else
2002 prc(k) = cons2/(cons3*cons18)*1350._r8
2003 & * qcic(i,k)**2.47_r8
2004 & * (1.e-6_r8*ncic(i,k)*rho(i,k))**(-1.79_r8)
2005 endif
2006
2007 nprc(k) = prc(k)/cons22
2008
2009 nprc1(k) = ncic(i,k)*prc(k) / (qcic(i,k)*(one+xs))
2010
2011 else
2012 prc(k) = zero
2013 nprc(k) = zero
2014 nprc1(k) = zero
2015 endif
2016
2017
2018
2019! add autoconversion to precip from above to get provisional rain mixing ratio
2020! and number concentration (qric and nric)
2021
2022! 0.45 m/s is fallspeed of new rain drop (80 micron diameter)
2023
2024 if (fprcp == 1) then
2025 tx1 = one / lcldm(i,k)
2026 qric(i,k) = min(qrn(i,k)*tx1, 10.e-3_r8)
2027 nric(i,k) = max(nrn(i,k)*tx1, zero)
2028 else
2029 dum = 0.45_r8
2030 dum1 = 0.45_r8
2031
2032 if (k == 1) then
2033 tx1 = lcldm(i,k)*dz(i,k)/(cldmax(i,k)*dum)
2034 qric(i,k) = prc(k) * tx1
2035 nric(i,k) = nprc(k) * tx1
2036 else
2037 if (qric(i,km) >= qsmall) then
2038! dum=umr(k-1)
2039! dum1=unr(k-1)
2040! Anning Cheng find a possible untable case here
2041 dum = max(umr(km),dum)
2042 dum1 = max(unr(km),dum1)
2043 endif
2044
2045! no autoconversion of rain number if rain/snow falling from above
2046! this assumes that new drizzle drops formed by autoconversion are rapidly collected
2047! by the existing rain/snow particles from above
2048
2049 if (qric(i,km) >= 1.e-9_r8 .or.
2050 & qniic(i,km) >= 1.e-9_r8) then
2051 nprc(k) = zero
2052 endif
2053
2054 tx1 = rho(i,km) * cldmax(i,km)
2055 tx3 = rho(i,k) * cldmax(i,k)
2056 qric(i,k) = (tx1*umr(km)*qric(i,km)
2057 & + (rdz*((pra(km)+prc(k))*lcldm(i,k)
2058 & + (pre(km)-pracs(km)-mnuccr(km))*cldmax(i,k))))
2059 & / (dum*tx3)
2060
2061
2062! nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+
2063! & (rho(i,k)*dz(i,k)*(nprc(k)*lcldm(i,k)+(nsubr(k-1)-npracs(k-1)-
2064! &nnuccr(k-1)+nragg(k-1))*cldmax(i,k)))) /(dum1*rho(i,k)*cldmax(i,
2065! &k))
2066
2067! Anning nsubr never given a value before
2068 nric(i,k) = (tx1*unr(km)*nric(i,km)
2069 & + (rdz*(nprc(k)*lcldm(i,k)
2070 & +(-npracs(km)-nnuccr(km)+nragg(km))*cldmax(i,k))))
2071 & / (dum1*tx3)
2072
2073 endif
2074 endif !fprcp
2075
2076!.......................................................................
2077! Autoconversion of cloud ice to snow (Ice_aut)
2078! similar to Ferrier (1994)
2079
2080 if (t(i,k) <= 273.15_r8 .and. qiic(i,k) >= qsmall) then
2081
2082! note: assumes autoconversion timescale of 180 sec
2083!vaux = 180.0_r8*10.0_r8
2084
2085 if (.false.) then
2086 vaux = ts_auto_ice * 10.0_r8
2087
2088 nprci(k) = (niic(i, k)/vaux)*exp(-lami(k)*dcs)
2089 tx1 = one / lami(k)
2090 tx2 = tx1 * tx1
2091 prci(k) = pi*irho(i,k)*niic(i,k)*lami(k)
2092 & / (6._r8*vaux)
2093 & * (cons23*tx1+three*cons24*tx2
2094 & + 6._r8*dcs*tx1*tx2+6._r8*tx2*tx2)
2095 & * exp(-lami(k)*dcs)
2096
2097 else
2098
2099! miu_ice(k) = mui_hemp_l(lami(k))
2100 miu_ice(k) =max(min(0.008_r8*(lami(k)*0.01)**0.87_r8
2101 &, 10.0_r8), 0.1_r8)
2102 tx1 = lami(k)*dcs
2103 nprci(k) = (niic(i,k)/ts_auto_ice)
2104 & * (one - gamma_incomp(miu_ice(k), tx1))
2105
2106 prci(k) = (qiic(i,k)/ts_auto_ice)
2107 & * (one - gamma_incomp(miu_ice(k)+three, tx1))
2108
2109 end if
2110 else
2111 prci(k) = zero
2112 nprci(k) = zero
2113 end if
2114
2115
2116! add autoconversion to flux from level above to get provisional snow mixing ratio
2117! and number concentration (qniic and nsic)
2118! Anning Cheng 9/16/2016 forecasting rain and snow, corresponding to MG2
2119 if (fprcp == 1) then
2120 tx1 = one / icldm(i,k)
2121 qniic(i,k) = min(qsnw(i,k)*tx1, 10.e-3_r8)
2122 nsic(i,k) = max(nsnw(i,k)*tx1, 0._r8)
2123 else
2124
2125 dum = (asn(i,k)*cons25)
2126 dum1 = (asn(i,k)*cons25)
2127
2128 if (k == 1) then
2129 tx1 = icldm(i,k)*dz(i,k)/(cldmax(i,k)*dum)
2130 qniic(i,k) = prci(k) * tx1
2131 nsic(i,k) = nprci(k) * tx1
2132 else
2133 if (qniic(i,km) >= qsmall) then
2134 dum = ums(km)
2135 dum1 = uns(km)
2136 end if
2137
2138 tx1 = rho(i,km) * cldmax(i,km)
2139 tx3 = rho(i,k) * cldmax(i,k)
2140 qniic(i,k) = (tx1*ums(km)*qniic(i,km) + (rdz*
2141 & ((prci(k)+prai(km)+psacws(km)+bergs(km))*icldm(i,k)
2142 & +(prds(km)+pracs(km)+mnuccr(km))*cldmax(i,k))))
2143 & / (dum*tx3)
2144
2145! nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+
2146! & (rho(i,k)*dz(i,k)*(nprci(k)*icldm(i,k)+(nsubs(k-1)+nsagg(k-1)+
2147! &nnuccr(k-1))*cldmax(i,k)))) /(dum1*rho(i,k)*cldmax(i,k))
2148
2149! nsubs never given a value before
2150 nsic(i,k) = (tx1*uns(km)*nsic(i,km)
2151 & + (rdz*(nprci(k)*icldm(i,k)+(nsagg(km)
2152 & + nnuccr(km))*cldmax(i,k)))) /(dum1*tx3)
2153
2154 end if
2155 end if !fprcp
2156
2157! if precip mix ratio is zero so should number concentration
2158
2159 if (qniic(i,k) < qsmall) then
2160 qniic(i,k) = zero
2161 nsic(i,k) = zero
2162 end if
2163
2164 if (qric(i,k) < qsmall) then
2165 qric(i,k) = zero
2166 nric(i,k) = zero
2167 end if
2168
2169! make sure number concentration is a positive number to avoid
2170! taking root of negative later
2171
2172 nric(i,k) = max(nric(i,k),zero)
2173 nsic(i,k) = max(nsic(i,k),zero)
2174
2175!.......................................................................
2176! get size distribution parameters for precip
2177!......................................................................
2178! rain
2179
2180 if (qric(i,k) >= qsmall) then
2181 lamr(k) = (pirhow*nric(i,k)/qric(i,k))**oneb3
2182 n0r(k) = nric(i,k)*lamr(k)
2183
2184! check for slope
2185! adjust vars
2186
2187 if (lamr(k) < lamminr) then
2188 lamr(k) = lamminr
2189 tx1 = lamminr * lamminr
2190 n0r(k) = tx1 * tx1 * qric(i,k)/pirhow
2191 nric(i,k) = n0r(k)/lamr(k)
2192 else if (lamr(k) > lammaxr) then
2193 lamr(k) = lammaxr
2194 tx1 = lammaxr * lammaxr
2195 n0r(k) = tx1 * tx1 * qric(i,k)/pirhow
2196 nric(i,k) = n0r(k)/lamr(k)
2197 end if
2198
2199! provisional rain number and mass weighted mean fallspeed (m/s)
2200
2201 tx1 = arn(i,k) / lamr(k) ** br
2202 tx2 = 9.1_r8*rhof(i,k)
2203 unr(k) = min(tx1*cons4, tx2)
2204 umr(k) = min(tx1*(cons5/6._r8), tx2)
2205
2206 else
2207 lamr(k) = zero
2208 n0r(k) = zero
2209 umr(k) = zero
2210 unr(k) = zero
2211 end if
2212
2213!......................................................................
2214! snow
2215
2216 if (qniic(i,k) >= qsmall) then
2217 lams(k) = (cons6*cs*nsic(i,k) / qniic(i,k))**(one/ds)
2218 n0s(k) = nsic(i,k)*lams(k)
2219
2220! check for slope
2221! adjust vars
2222
2223 if (lams(k) < lammins) then
2224 lams(k) = lammins
2225 n0s(k) = lams(k)**(ds+one)*qniic(i,k)/(cs*cons6)
2226 nsic(i,k) = n0s(k)/lams(k)
2227 else if (lams(k) > lammaxs) then
2228 lams(k) = lammaxs
2229 n0s(k) = lams(k)**(ds+one)*qniic(i,k)/(cs*cons6)
2230 nsic(i,k) = n0s(k)/lams(k)
2231 end if
2232
2233! provisional snow number and mass weighted mean fallspeed (m/s)
2234
2235 tx1 = asn(i,k) / lams(k)**bs
2236 tx2 = 1.2_r8*rhof(i,k)
2237 ums(k) = min(tx1*(cons8/6._r8), tx2)
2238 uns(k) = min(tx1*cons7, tx2)
2239
2240 else
2241 lams(k) = zero
2242 n0s(k) = zero
2243 ums(k) = zero
2244 uns(k) = zero
2245 end if
2246
2247!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2248checked up to here moorthi
2249
2250! heterogeneous freezing of cloud water
2251
2252 if (qcic(i,k) >= qsmall .and. t(i,k) < 269.15_r8) then
2253
2254! immersion freezing (Bigg, 1953)
2255
2256 tx2 = one / (lamc(k) * lamc(k) * lamc(k))
2257 tx3 = min(aimm*(273.15_r8-t(i,k)), 25.0)
2258 tx1 = cdist1(k) * tx2 * bimm * exp(tx3)
2259 mnuccc(k) = cons9/(cons3*cons19)* pi*pirhow/36._r8
2260 & * gamma(7._r8+pgam(k))*tx1*tx2
2261
2262 nnuccc(k) = cons10/(cons3*qcvar)* pi/6._r8
2263 & * gamma(pgam(k)+four) * tx1
2264
2265 if (.true.) then
2266 nnuccc(k) = nimm(i,k)
2267 mnuccc(k) = nimm(i,k)*qcic(i,k)/max(cdist1(k), one)
2268 end if
2269
2270! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
2271! dust size and number in 4 bins are read in from companion routine
2272
2273!DONIFFFF this implementation would oveestimate strongly contact IN
2274! According to Young 1974 the base concentration must be the active contact IN at -4 C
2275! This implementation is using all the dust at all levels
2276!Freezing fraction of collected IN is assumed to be 1.
2277
2278 tcnt = (270.16_r8-t(i,k))**1.3_r8
2279 viscosity = 1.8e-5_r8*(t(i,k)/298.0_r8)**0.85_r8
2280 mfp = two*viscosity
2281 & / (p(i,k) *sqrt(8.0_r8*28.96e-3_r8/(pi*
2282 & 8.314409_r8*t(i,k))))
2283 taux = t(i,k) - three
2284 taux = max(taux-273.15, -40.0)
2285
2286 nsdust = max(exp(-0.517*taux + 8.934) -3.76e6, 0.0)
2287 nssoot = max(1.0e4*exp((-0.0101*taux-0.8525)*taux
2288 & +0.7667) -3.77e9, 0.0)
2289
2290 do ind_aux = 1, nbincontactdust
2291 tx1 = rndst(i,k,ind_aux)
2292 nslip(ind_aux) = one + (mfp/tx1) *
2293 & (1.257_r8+(0.4_r8*exp(-(1.1_r8*tx1/mfp))))
2294 ndfaer(ind_aux) = 1.381e-23_r8*t(i,k)* nslip(ind_aux)
2295 & / (6._r8*pi*viscosity*tx1)
2296
2297 ndfaer(ind_aux) = ndfaer(ind_aux)
2298 & * (1.0-exp(-nsdust*12.5664*tx1*tx1))
2299 end do
2300
2301 nslipsoot = one + (mfp/rnsoot(i,k)) *
2302 & (1.257_r8+(0.4_r8*exp(-(1.1_r8*rnsoot(i,k)/mfp))))
2303 ndfaersoot = 1.381e-23_r8*t(i,k)*nslipsoot
2304 & / (6._r8*pi*viscosity*rnsoot(i,k))
2305 ndfaersoot = ndfaersoot *
2306 & (1.0-exp(-nssoot*12.5664*rnsoot(i,k)*rnsoot(i,k)))
2307
2308
2309
2310 tx1 = sum(ndfaer(:)*nacon(i,k,:))+ndfaersoot*nsoot(i,k)
2311 tx1 = tx1 * pi * cdist1(k)
2312 tx2 = one / lamc(k)
2313 tx3 = tx2 * tx2
2314 mnucct(k) = (pi*cons10/(cons3*qcvar*three)) * tx1
2315 & * rhow*gamma(pgam(k)+five)*tx3*tx3
2316
2317
2318 nnucct(k) = (tx1+tx1) * gamma(pgam(k)+two) * tx2
2319
2320! make sure number of droplets frozen does not exceed available ice nuclei concentration
2321! this prevents 'runaway' droplet freezing
2322#ifdef CAM
2323 if (nnuccc(k)*lcldm(i,k) > nnuccd(k)) then
2324 dum = (nnuccd(k)/(nnuccc(k)*lcldm(i,k)))
2325! scale mixing ratio of droplet freezing with limit
2326 mnuccc(k) = mnuccc(k)*dum
2327 nnuccc(k) = nnuccd(k)/lcldm(i,k)
2328 end if
2329#endif
2330
2331
2332 if ((nnucct(k)+nnuccc(k))*deltat > ncic(i, k)) then
2333
2334 tx1 = tx2 * tx3
2335 nnuccc(k) = ncic(i, k)*dti
2336 mnuccc(k) = cons9/(cons3*cons19)* pi*pirhow/36._r8
2337 & * cdist1(k)*gamma(7._r8+pgam(k))*ncic(i,k)
2338 & * tx1 * tx1
2339
2340 nnucct(k) = zero
2341 mnucct(k) = zero
2342 end if
2343
2344
2345 else
2346 mnuccc(k) = zero
2347 nnuccc(k) = zero
2348 mnucct(k) = zero
2349 nnucct(k) = zero
2350 end if
2351
2352!.......................................................................
2353! snow self-aggregation from passarelli, 1978, used by reisner, 1998
2354! this is hard-wired for bs = 0.4 for now
2355! ignore self-collection of cloud ice
2356
2357 if (qniic(i,k) >= qsmall .and. t(i,k) <= 273.15_r8) then
2358 tx1 = (two+bs)*oneb3
2359 nsagg(k) = -1108._r8*asn(i,k)*eii* pi**((one-bs)*oneb3)
2360 & * (rho(i,k)*qniic(i,k)/rhosn) ** tx1
2361 & * (nsic(i,k)*rho(i,k))**((four-bs)*oneb3)
2362 & / (four*720._r8*rho(i,k))
2363 else
2364 nsagg(k) = zero
2365 end if
2366
2367!.......................................................................
2368! accretion of cloud droplets onto snow/graupel
2369! here use continuous collection equation with
2370! simple gravitational collection kernel
2371! ignore collisions between droplets/cloud ice
2372! since minimum size ice particle for accretion is 50 - 150 micron
2373
2374! ignore collision of snow with droplets above freezing
2375
2376 if (qniic(i,k) >= qsmall .and. t(i,k) <= tmelt
2377 & .and. qcic(i,k) >= qsmall) then
2378
2379! put in size dependent collection efficiency
2380! mean diameter of snow is area-weighted, since
2381! accretion is function of crystal geometric area
2382! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)
2383
2384 dc0 = (pgam(k)+one)/lamc(k)
2385 ds0 = one/lams(k)
2386 dum = dc0*dc0*uns(k)*rhow/(9._r8*mu(i,k)*ds0)
2387 eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
2388
2389 eci = min(one, max(eci, zero))
2390
2391
2392! no impact of sub-grid distribution of qc since psacws
2393! is linear in qc
2394
2395 tx1 = pi/four*asn(i,k)*rho(i,k)*n0s(k)*eci*cons11
2396 & / lams(k)**(bs+three)
2397 psacws(k) = tx1 * qcic(i,k)
2398 npsacws(k) = tx1 * ncic(i,k)
2399 else
2400 psacws(k) = zero
2401 npsacws(k) = zero
2402 end if
2403
2404! add secondary ice production due to accretion of droplets by snow
2405! (Hallet-Mossop process) (from Cotton et al., 1986)
2406
2407 if((t(i,k) < 270.16_r8) .and. (t(i,k) >= 268.16_r8)) then
2408 ni_secp = 0.5*3.5e8_r8*(270.16_r8-t(i,k))*psacws(k)
2409 nsacwi(k) = ni_secp
2410 msacwi(k) = min(ni_secp*mi0, psacws(k))
2411 else if((t(i,k) < 268.16_r8) .and.
2412 & (t(i,k) >= 265.16_r8)) then
2413 ni_secp = oneb3*3.5e8_r8*(t(i,k)-265.16_r8)*psacws(k)
2414 nsacwi(k) = ni_secp
2415 msacwi(k) = min(ni_secp*mi0, psacws(k))
2416 else
2417 ni_secp = zero
2418 nsacwi(k) = zero
2419 msacwi(k) = zero
2420 endif
2421 psacws(k) = max(zero, psacws(k)-ni_secp*mi0)
2422
2423!.......................................................................
2424! accretion of rain water by snow
2425! formula from ikawa and saito, 1991, used by reisner et al., 1998
2426
2427 if (qric(i,k) >= 1.e-8_r8 .and. qniic(i,k) >= 1.e-8_r8
2428 & .and. t(i,k) <= 273.15_r8) then
2429! Anning decrease pracs it can reach 2.3e5 so/1.e6
2430 tx1 = 1.2_r8*umr(k) - 0.95_r8*ums(k)
2431 tx2 = sqrt(tx1*tx1+0.08_r8*ums(k)*umr(k))
2432 tx1 = one / lamr(k)
2433 tx3 = one / lams(k)
2434 tx4 = tx1 * tx1
2435 tx5 = pi * ecr * rho(i,k) *n0r(k) * n0s(k)
2436
2437 pracs(k) = pirhow*tx2*tx5 *
2438 & (tx4*tx4*tx3*(five*tx4+tx3*(two*tx1+half*tx3)))
2439
2440
2441 tx2 = unr(k) - uns(k)
2442 tx2 = sqrt(1.7_r8*tx2*tx2 + 0.3_r8*unr(k)*uns(k))
2443
2444 npracs(k) = half*tx2*tx5*tx1*tx3*(tx4+tx3*(tx1+tx3))
2445
2446 else
2447 pracs(k) = zero
2448 npracs(k) = zero
2449 end if
2450
2451!.......................................................................
2452! heterogeneous freezing of rain drops
2453! follows from Bigg (1953)
2454
2455! if (t(i,k).lt.269.15_r8 .and. qric(i,k).ge.qsmall) then
2456! Anning change to prevent huge value of mnuccr
2457
2458 if (t(i,k) < 269.15_r8 .and. qric(i,k) >= qsmall
2459 & .and. t(i,k) > 260.0_r8) then
2460
2461! tx1 = exp(min(aimm*(273.15_r8-t(i,k))-one, 50.0))
2462 tx1 = exp(min(aimm*(273.15_r8-t(i,k))-one, 25.0))
2463 tx2 = 1.0 / (lamr(k)*lamr(k)*lamr(k))
2464
2465 nnuccr(k) = pi * nric(i,k) * bimm * tx1 * tx2
2466 mnuccr(k) = 20._r8 * pirhow * nnuccr(k) * tx2
2467
2468 else
2469 mnuccr(k) = zero
2470 nnuccr(k) = zero
2471
2472 end if
2473
2474!.......................................................................
2475! accretion of cloud liquid water by rain
2476! formula from Khrouditnov and Kogan (2000)
2477! gravitational collection kernel, droplet fall speed neglected
2478
2479 if (qric(i,k) >= qsmall .and. qcic(i,k) >= qsmall) then
2480
2481! include sub-grid distribution of cloud water
2482
2483 pra(k) = cons12/(cons3*cons20)
2484 & * 67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
2485 npra(k) = pra(k) * (ncic(i,k)/qcic(i,k))
2486
2487 else
2488 pra(k) = zero
2489 npra(k) = zero
2490 end if
2491
2492!Not done above
2493!.......................................................................
2494! Self-collection of rain drops
2495! from Beheng(1994)
2496
2497 if (qric(i,k) >= qsmall) then
2498 nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k)
2499 else
2500 nragg(k) = zero
2501 end if
2502
2503!.......................................................................
2504! Accretion of cloud ice by snow
2505! For this calculation, it is assumed that the Vs >> Vi
2506! and Ds >> Di for continuous collection
2507
2508 if (qniic(i,k) >= qsmall .and. qiic(i,k) >= qsmall
2509 & .and.t(i,k) <= 273.15_r8) then
2510
2511 tx1 = (pi/four)*asn(i,k)*rho(i,k)*n0s(k)*eii*cons11
2512 & / lams(k)**(bs+three)
2513 prai(k) = tx1 * qiic(i,k)
2514 nprai(k) = tx1 * niic(i,k)
2515
2516 nprai(k)= min(nprai(k), 1.0e10)
2517
2518 else
2519 prai(k) = zero
2520 nprai(k) = zero
2521 end if
2522
2523!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2524! calculate evaporation/sublimation of rain and snow
2525! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
2526! in-cloud condensation/deposition of rain and snow is neglected
2527! except for transfer of cloud water to snow through bergeron process
2528
2529! initialize evap/sub tendncies
2530 pre(k) = zero
2531 prds(k) = zero
2532
2533! evaporation of rain
2534! only calculate if there is some precip fraction > cloud fraction
2535
2536 if (qcic(i,k)+qiic(i,k) < 1.e-7_r8 .or.
2537 & cldmax(i,k) > lcldm(i,k)) then
2538
2539! set temporary cloud fraction to zero if cloud water + ice is very small
2540! this will ensure that evaporation/sublimation of precip occurs over
2541! entire grid cell, since min cloud fraction is specified otherwise
2542 if (qcic(i,k)+qiic(i,k) < 1.e-7_r8) then
2543 dum = zero
2544 else
2545 dum = lcldm(i,k)
2546 end if
2547
2548! saturation vapor pressure
2549! esn = polysvp(t(i,k),0)
2550 esn = min(fpvsl(t(i,k)), p(i,k))
2551 qsn = min(epsqs*esn/(p(i,k)-omeps*esn), one)
2552
2553! recalculate saturation vapor pressure for liquid and ice
2554 esl(i,k) = esn
2555! esi(i,k) = polysvp(t(i,k),1)
2556 esi(i,k) = min(fpvsi(t(i,k)), p(i,k))
2557! hm fix, make sure when above freezing that esi=esl, not active yet
2558 if (t(i,k) > tmelt) esi(i,k) = esl(i,k)
2559
2560! calculate q for out-of-cloud region
2561
2562 qclr = (q(i,k)-dum*qsn) / (one-dum)
2563
2564 if (qric(i,k) >= qsmall) then
2565
2566 qvs = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
2567 dqsdt = xxlv*qvs/(rv*t(i,k)*t(i,k))
2568 ab = one + dqsdt*xxlv/cpp
2569 epsr = (pi+pi)*n0r(k)*rho(i,k)*dv(i,k)
2570 & * (f1r/(lamr(k)*lamr(k))
2571 & + f2r*sqrt(arn(i,k)*rho(i,k)/mu(i,k))
2572 & * sc(i,k)**oneb3
2573 & * cons13 / lamr(k)**((five+br)*half))
2574
2575 pre(k) = epsr*(qclr-qvs)/ab
2576
2577! only evaporate in out-of-cloud region
2578! and distribute across cldmax
2579 pre(k) = min(pre(k)*(cldmax(i,k)-dum), zero)
2580 pre(k) = pre(k) / cldmax(i,k)
2581
2582 end if
2583
2584! sublimation of snow
2585 if (qniic(i,k) >= qsmall) then
2586 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
2587 dqsidt = xxls*qvi/(rv*t(i,k)*t(i,k))
2588 abi = one + dqsidt*xxls/cpp
2589 epss = (pi+pi)*n0s(k)*rho(i,k)*dv(i,k)
2590 & * (f1s/(lams(k)*lams(k))
2591 & + f2s*sqrt((asn(i,k)*rho(i,k)/mu(i,k)))
2592 & * sc(i,k)**oneb3
2593 & * cons14/ lams(k)**((five+bs)*half))
2594 prds(k) = epss*(qclr-qvi)/abi
2595
2596! only sublimate in out-of-cloud region and distribute over cldmax
2597 prds(k) = min(prds(k)*(cldmax(i,k)-dum), zero)
2598 prds(k) = prds(k)/cldmax(i,k)
2599 end if
2600
2601
2602! make sure RH not pushed above 100% due to rain evaporation/snow sublimation
2603! get updated RH at end of time step based on cloud water/ice condensation/evap
2604
2605 tx1 = pre(k) * cldmax(i,k)
2606 tx2 = prds(k) * cldmax(i,k)
2607 qtmp = q(i,k) - (cmei(i,k)+tx1+tx2) * deltat
2608 ttmp = t(i,k) + (tx1*xxlv + (cmei(i,k)+tx2)*xxls)
2609 & * (deltat/cpp)
2610
2611!limit range of temperatures!
2612 ttmp = max(180._r8,min(ttmp,323._r8))
2613
2614! esn = polysvp(ttmp,0)
2615 esn = min(fpvsl(ttmp), p(i,k))
2616 qsn = min(epsqs*esn/(p(i,k)-omeps*esn), one)
2617
2618! modify precip evaporation rate if q > qsat
2619 if (qtmp > qsn) then
2620 if (pre(k)+prds(k) < -1.e-20) then
2621 dum1 = pre(k) / (pre(k)+prds(k))
2622! recalculate q and t after cloud water cond but without precip evap
2623 qtmp = q(i,k) - cmei(i,k)*deltat
2624 ttmp = t(i,k) + cmei(i,k)*xxls*deltat/cpp
2625! esn = polysvp(ttmp,0)
2626 esn = min(fpvsl(ttmp), p(i,k))
2627 qsn = min(epsqs*esn/(p(i,k)-omeps*esn), one)
2628 tx1 = one / (cpp*rv*ttmp*ttmp)
2629 dum = min(zero, (qtmp-qsn)/(one+cons27*qsn*tx1))
2630
2631! modify rates if needed, divide by cldmax to get local (in-precip) value
2632 pre(k) = dum*dum1/(deltat*cldmax(i,k))
2633
2634! do separately using RHI for prds....
2635! esn = polysvp(ttmp,1)
2636 esn = min(fpvsi(ttmp), p(i,k))
2637 qsn = min(epsqs*esn/(p(i,k)-omeps*esn), one)
2638 dum = min(zero, (qtmp-qsn)/(one+cons28*qsn*tx1))
2639
2640! modify rates if needed, divide by cldmax to get local (in-precip) value
2641 prds(k) = dum*(one-dum1)/(deltat*cldmax(i,k))
2642 end if
2643 end if
2644
2645 end if
2646
2647! bergeron process - evaporation of droplets and deposition onto snow
2648 if (qniic(i,k) >= qsmall .and. qcic(i,k) >= qsmall
2649 & .and. t(i,k) < tmelt) then
2650 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
2651 qvs = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
2652 dqsidt = xxls*qvi/(rv*t(i,k)*t(i,k))
2653 abi = one + dqsidt*xxls/cpp
2654 epss = (pi+pi)*n0s(k)*rho(i,k)*dv(i,k)
2655 & * (f1s/(lams(k)*lams(k))
2656 & + f2s*sqrt(asn(i,k)*rho(i,k)/mu(i,k))
2657 & * sc(i,k)**oneb3
2658 & * cons14 / (lams(k)**((five+bs)*half)))
2659 bergs(k) = epss*(qvs-qvi)/abi
2660 else
2661 bergs(k) = zero
2662 end if
2663
2664
2665
2666!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2667! conservation to ensure no negative values of cloud water/precipitation
2668! in case microphysical process rates are large
2669
2670! make sure and use end-of-time step values for cloud water, ice, due
2671! condensation/deposition
2672
2673! note: for check on conservation, processes are multiplied by omsm
2674! to prevent problems due to round off error
2675
2676! include mixing timescale (mtime)
2677
2678 qce = max(qc(i,k)-berg(i,k)*deltat, zero)
2679 nce = nc(i,k) + npccn(i,k)*deltat*mtime
2680 qie = qi(i,k) + (cmei(i,k)+berg(i,k))*deltat
2681 nie = ni(i,k) + nnuccd(k)*deltat*mtime
2682
2683! conservation of qc
2684
2685 tx1 = lcldm(i,k) * deltat
2686 dum = (prc(k) + pra(k) + mnuccc(k) + mnucct(k)
2687 & + msacwi(k) + psacws(k) + bergs(k)) * tx1
2688
2689 if (dum > qce) then
2690
2691 ratio = qce/dum * omsm
2692
2693 prc(k) = prc(k) * ratio
2694 pra(k) = pra(k) * ratio
2695 mnuccc(k) = mnuccc(k) * ratio
2696 mnucct(k) = mnucct(k) * ratio
2697 msacwi(k) = msacwi(k) * ratio
2698 psacws(k) = psacws(k) * ratio
2699 bergs(k) = bergs(k) * ratio
2700 end if
2701
2702! conservation of nc
2703
2704
2705 dum = (nprc1(k) + npra(k) + nnuccc(k) + nnucct(k)
2706 & + npsacws(k) - nsubc(k)) * tx1
2707
2708 if (dum > nce) then
2709 ratio = nce/dum * omsm
2710
2711 nprc1(k) = nprc1(k) * ratio
2712 npra(k) = npra(k) * ratio
2713 nnuccc(k) = nnuccc(k) * ratio
2714 nnucct(k) = nnucct(k) * ratio
2715 npsacws(k) = npsacws(k) * ratio
2716 nsubc(k) = nsubc(k) * ratio
2717 end if
2718
2719! conservation of qi
2720
2721 dum = ((-mnuccc(k)-mnucct(k)-msacwi(k))*lcldm(i,k)
2722 & + (prci(k)+prai(k))*icldm(i,k)) * deltat
2723
2724 if (dum > qie) then
2725 if (prci(k)+prai(k) > zero) then
2726 ratio = (qie*dti
2727 & +(mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k))
2728 & / ((prci(k)+prai(k))*icldm(i,k))*omsm
2729 else
2730 ratio = zero
2731 end if
2732
2733 prci(k) = prci(k)*ratio
2734 prai(k) = prai(k)*ratio
2735 end if
2736
2737! conservation of ni
2738
2739 dum = ((-nnucct(k)-nsacwi(k))*lcldm(i,k)
2740 & + (nprci(k)+ nprai(k)-nsubi(k))*icldm(i,k))*deltat
2741
2742 if (dum > nie) then
2743 if ( abs(nprci(k)+nprai(k)-nsubi(k)) > zero) then
2744 ratio = (nie*dti+(nnucct(k)+nsacwi(k))*lcldm(i,k))
2745 & / ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm
2746 else
2747
2748 ratio = zero
2749 end if
2750 nprci(k) = nprci(k) * ratio
2751 nprai(k) = nprai(k) * ratio
2752 nsubi(k) = nsubi(k) * ratio
2753 end if
2754
2755! for preciptiation conservation, use logic that vertical integral
2756! of tendency from current level to top of model (i.e., qrtot) cannot be negative
2757
2758! conservation of rain mixing rat
2759
2760 if ( ((prc(k)+pra(k))*lcldm(i,k)
2761 & + (-mnuccr(k)+pre(k)-pracs(k))*cldmax(i,k))*rdz
2762 & + qrtot < zero) then
2763
2764 if (-pre(k)+pracs(k)+mnuccr(k) >= qsmall) then
2765
2766 ratio = (qrtot*rdzi + (prc(k)+pra(k))*lcldm(i,k))
2767 & / ((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm
2768
2769 pre(k) = pre(k) * ratio
2770 pracs(k) = pracs(k) * ratio
2771 mnuccr(k) = mnuccr(k) * ratio
2772 end if
2773 end if
2774
2775! conservation of nr - for now neglect evaporation of nr
2776
2777 nsubr(k) = zero
2778
2779 if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k) +
2780 & nragg(k))*cldmax(i,k))*rdz + nrtot < zero) then
2781
2782 if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k) >= qsmall)
2783 & then
2784 ratio = (nrtot*rdzi+nprc(k)*lcldm(i,k))
2785 & / ((-nsubr(k)-nragg(k)+npracs(k)
2786 & +nnuccr(k))*cldmax(i,k))*omsm
2787
2788 nsubr(k) = nsubr(k) * ratio
2789 npracs(k) = npracs(k) * ratio
2790 nnuccr(k) = nnuccr(k) * ratio
2791 nragg(k) = nragg(k) * ratio
2792 end if
2793 end if
2794
2795! conservation of snow mix ratio
2796
2797 tx1 = (bergs(k)+psacws(k))*lcldm(i,k)
2798 & + (prai(k)+prci(k))*icldm(i,k)
2799 if ((tx1+(pracs(k)+ mnuccr(k)+prds(k))*cldmax(i,k))
2800 & *rdz+qstot < zero) then
2801
2802 if (-prds(k) >= qsmall) then
2803
2804 ratio = (qstot*rdzi + tx1
2805 & + (pracs(k)+mnuccr(k))*cldmax(i,k))
2806 & / (-prds(k)*cldmax(i,k))*omsm
2807
2808 prds(k) = prds(k) * ratio
2809 end if
2810 end if
2811
2812! conservation of ns
2813
2814! calculate loss of number due to sublimation
2815! for now neglect sublimation of ns
2816 nsubs(k) = zero
2817
2818 if ((nprci(k)*icldm(i,k)
2819 & +(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))
2820 & * rdz + nstot < zero) then
2821
2822 if (-nsubs(k)-nsagg(k) >= qsmall) then
2823
2824 ratio = (nstot*rdzi
2825 & + nprci(k)*icldm(i,k)+ nnuccr(k)*cldmax(i,k))
2826 & / ((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm
2827
2828 nsubs(k) = nsubs(k) * ratio
2829 nsagg(k) = nsagg(k) * ratio
2830 end if
2831 end if
2832
2833! get tendencies due to microphysical conversion processes
2834! note: tendencies are multiplied by appropaiate cloud/precip
2835! fraction to get grid-scale values
2836! note: cmei is already grid-average values
2837
2838 qvlat(i,k) = qvlat(i,k)
2839 & - (pre(k)+prds(k))*cldmax(i,k) - cmei(i,k)
2840! if (lprint .and. k == 29) write(0,*)' qvlata=',qvlat(i,k),
2841! &' pre=',pre(k),' prds=',prds(k),' cldmax=',cldmax(i,k),cmei(i,k)
2842! &,' it=',it
2843
2844 tlat(i,k) = tlat(i,k)
2845 & + pre(k)*cldmax(i,k) * xxlv
2846 & + (prds(k)*cldmax(i,k)+cmei(i,k)) * xxls
2847 & + ((bergs(k)+psacws(k)+mnuccc(k)+
2848 & mnucct(k)+msacwi(k))*lcldm(i,k)
2849 & + (mnuccr(k)+ pracs(k))*cldmax(i,k)
2850 & + berg(i,k)) * xlf
2851
2852
2853! if(xlon<0.01.and.xlon>-0.01.and.xlat>1.346
2854! & .and.xlat<1.347.and.k==38)
2855! & write(*,*)"anning_m0",pre(k),prds(k),cmei(i,k),
2856! & bergs(k),psacws(k),
2857! & mnuccc(k),mnucct(k),msacwi(k),mnuccr(k),pracs(k),berg(i,k)
2858
2859 qctend(i,k) = qctend(i,k)
2860 & + (-pra(k)-prc(k)-mnuccc(k)-mnucct(k)
2861 & -msacwi(k)-psacws(k)-bergs(k))*lcldm(i,k)
2862 & - berg(i,k)
2863
2864 qitend(i,k) = qitend(i,k)
2865 & + (mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)
2866 & + (-prci(k)-prai(k))*icldm(i,k)
2867 & + cmei(i,k) + berg(i,k)
2868
2869 qrtend(i,k) = qrtend(i,k)
2870 & + (pra(k)+prc(k))*lcldm(i,k)
2871 & + (pre(k)-pracs(k)-mnuccr(k))*cldmax(i,k)
2872
2873 qnitend(i,k) = qnitend(i,k)
2874 & + (prai(k)+prci(k))*icldm(i,k)
2875 & + (psacws(k)+bergs(k))*lcldm(i,k)
2876 & + (prds(k)+pracs(k)+mnuccr(k))*cldmax(i,k)
2877
2878! if (lprint) write(0,*)' k=',k,' qnitend=',qnitend(i,k),
2879! & prai(k), prci(k), icldm(i,k),psacws(k),bergs(k),lcldm(i,k)
2880! &,prds(k),pracs(k),mnuccr(k),' cldmax=',cldmax(i,k)
2881
2882! add output for cmei (accumulate)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2883 cmeiout(i,k) = cmeiout(i,k) + cmei(i,k)
2884
2885!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2886
2887
2888
2889! assign variables for trop_mozart, these are grid-average
2890! evaporation/sublimation is stored here as positive term
2891
2892 evapsnow(i,k) = evapsnow(i,k) + prds(k) * cldmax(i,k)
2893 nevapr(i,k) = nevapr(i,k) + pre(k) * cldmax(i,k)
2894
2895! change to make sure prain is positive: do not remove snow from
2896! prain used for wet deposition
2897 prain(i,k) = prain(i,k)
2898 & + (pra(k)+prc(k))*lcldm(i,k)
2899 & + (-pracs(k)-mnuccr(k))*cldmax(i,k)
2900
2901 prodsnow(i,k) = prodsnow(i,k)
2902 & + (prai(k)+prci(k))*icldm(i,k)
2903 & + (psacws(k)+bergs(k))*lcldm(i,k)
2904 & + (pracs(k)+mnuccr(k))*cldmax(i,k)
2905
2906! following are used to calculate 1st order conversion rate of cloud water
2907! to rain and snow (1/s), for later use in aerosol wet removal routine
2908! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
2909! used to calculate pra, prc, ... in this routine
2910! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow }
2911! (no cloud ice or bergeron terms)
2912! qcsum_rate1ord = sum over iterations{ qc used in calculation of the transfer terms }
2913
2914 qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k)
2915 & +(pra(k)+prc(k)+psacws(k))*lcldm(i,k)
2916 qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k)
2917
2918! microphysics output, note this is grid-averaged
2919 prao(i,k) = prao(i,k) + pra(k) * lcldm(i,k)
2920 prco(i,k) = prco(i,k) + prc(k) * lcldm(i,k)
2921 mnuccco(i,k) = mnuccco(i,k) + mnuccc(k) * lcldm(i,k)
2922 mnuccto(i,k) = mnuccto(i,k) + mnucct(k) * lcldm(i,k)
2923 msacwio(i,k) = msacwio(i,k) + msacwi(k) * lcldm(i,k)
2924 psacwso(i,k) = psacwso(i,k) + psacws(k) * lcldm(i,k)
2925 bergso(i,k) = bergso(i,k) + bergs(k) * lcldm(i,k)
2926 bergo(i,k) = bergo(i,k) + berg(i,k)
2927 prcio(i,k) = prcio(i,k) + prci(k) * icldm(i,k)
2928 praio(i,k) = praio(i,k) + prai(k) * icldm(i,k)
2929 mnuccro(i,k) = mnuccro(i,k) + mnuccr(k) * cldmax(i,k)
2930 pracso(i,k) = pracso(i,k) + pracs(k) * cldmax(i,k)
2931
2932 mnuccdo(i,k) = mnuccdo(i,k) + mnuccd(k)
2933 nnuccto(i,k) = nnuccto(i,k) + nnucct(k) * lcldm(i,k)
2934
2935 nnuccdo(i,k) = nnuccdo(i,k) + nnuccd(k)
2936 nnuccco(i,k) = nnuccco(i,k) + nnuccc(k) * lcldm(i,k)
2937 nsacwio(i,k) = nsacwio(i,k) + nsacwi(k) * icldm(i,k)
2938 nsubio(i,k) = nsubio(i,k) + nsubi(k) * icldm(i,k)
2939 nprcio(i,k) = nprcio(i,k) + nprci(k) * icldm(i,k)
2940 npraio(i,k) = npraio(i,k) + nprai(k) * icldm(i,k)
2941
2942 npccno(i, k) = npccno(i,k) + npccn(i,k)*lcldm(i,k)
2943 npsacwso(i,k) = npsacwso(i,k) + npsacws(k)*lcldm(i,k)
2944 nsubco(i,k) = nsubco(i,k) + nsubc(k) * lcldm(i,k)
2945 nprao(i,k) = nprao(i,k) + npra(k) * lcldm(i,k)
2946 nprc1o(i,k) = nprc1o(i,k) + nprc1(k) * lcldm(i,k)
2947
2948
2949! multiply activation/nucleation by mtime to account for fast timescale
2950
2951 nctend(i,k) = nctend(i,k) + npccn(i,k)*mtime
2952 & + (-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k)
2953 & -npra(k)-nprc1(k))*lcldm(i,k)
2954
2955 nitend(i,k) = nitend(i,k) + nnuccd(k)*mtime
2956 & + (nnuccc(k)+nnucct(k)+nsacwi(k))*lcldm(i,k)
2957 & + (nsubi(k)-nprci(k)- nprai(k))*icldm(i,k)
2958
2959 nstend(i,k) = nstend(i,k)
2960 & + (nsubs(k)+ nsagg(k)+nnuccr(k))*cldmax(i,k)
2961 & + nprci(k)*icldm(i,k)
2962
2963 nrtend(i,k) = nrtend(i,k) + nprc(k)*lcldm(i,k)
2964 & + (nsubr(k)-npracs(k)-nnuccr(k)+nragg(k))
2965 & * cldmax(i,k)
2966
2967! make sure that nc and ni at advanced time step do not exceed
2968! maximum (existing N + source terms*dt), which is possible due to
2969! fast nucleation timescale
2970
2971 if (nctend(i,k) > zero .and.
2972 & nc(i,k)+nctend(i,k)*deltat > ncmax(i,k)) then
2973 nctend(i,k) = max(zero, (ncmax(i,k)-nc(i,k))*dti)
2974 end if
2975 if (nitend(i,k) > zero .and.
2976 & ni(i,k)+nitend(i,k)*deltat > nimax) then
2977 nitend(i,k) = max(zero, (nimax-ni(i,k))*dti)
2978 end if
2979
2980
2981! get final values for precipitation q and N, based on
2982! flux of precip from above, source/sink term, and terminal fallspeed
2983! see eq. 15-16 in MG2008
2984 if (fprcp == 0) then
2985! rain
2986
2987 if (qric(i,k) >= qsmall) then
2988 if (k == 1) then
2989 qric(i,k) = qrtend(i,k)*dz(i,k)/(cldmax(i,k)*umr(k))
2990 nric(i,k) = nrtend(i,k)*dz(i,k)/(cldmax(i,k)*unr(k))
2991 else
2992 tx1 = rho(i,km) * cldmax(i,km)
2993 tx3 = rho(i,k) * cldmax(i,k)
2994
2995 qric(i,k) = (tx1*umr(km)*qric(i,km)
2996 & + rdz*qrtend(i,k)) / (umr(k)*tx3)
2997
2998 nric(i,k) = (tx1*unr(km)*nric(i,km)
2999 & + rdz*nrtend(i,k)) / (unr(k)*tx3)
3000
3001 end if
3002 else
3003 qric(i,k) = zero
3004 nric(i,k) = zero
3005 end if
3006
3007! snow
3008
3009 if (qniic(i,k) >= qsmall) then
3010 if (k == 1) then
3011 tx1 = dz(i,k)/cldmax(i,k)
3012 qniic(i,k) = qnitend(i,k)*tx1/ums(k)
3013 nsic(i,k) = nstend(i,k)*tx1/uns(k)
3014 else
3015 tx1 = rho(i,km) * cldmax(i,km)
3016 tx3 = rho(i,k) * cldmax(i,k)
3017
3018 qniic(i,k) = (tx1*ums(km)*qniic(i,km)
3019 & + rdz*qnitend(i,k)) / (ums(k)*tx3)
3020 nsic(i,k) = (tx1*uns(km)*nsic(i,km)
3021 & + rdz*nstend(i,k)) / (uns(k)*tx3)
3022 end if
3023 else
3024 qniic(i,k) = zero
3025 nsic(i,k) = zero
3026 end if
3027
3028! calculate precipitation flux at surface
3029! divide by density of water to get units of m/s
3030
3031 tx1 = rdz/rhow
3032 prect(i) = prect(i) + (qrtend(i,k)+ qnitend(i,k))*tx1
3033 preci(i) = preci(i) + qnitend(i,k)*tx1
3034
3035! if (lprint) write(0,*)' prect=',prect(i),' preci=',preci(i)
3036! &,' qrtend=',qrtend(i,k),' qnitend=',qnitend(i,k),' rdz=',rdz
3037! &,' k=',k,' it=',it, 'rhow=',rhow,' rho=',rho(i,k),' dz=',dz(i,k)
3038
3039
3040
3041 rainrt(i,k) = qric(i,k)*rho(i,k)*umr(k)
3042 & / rhow*3600._r8*1000._r8
3043
3044! vertically-integrated precip source/sink terms (note: grid-averaged)
3045
3046 qrtot = max(qrtot + qrtend(i,k)*rdz, zero)
3047 qstot = max(qstot + qnitend(i,k)*rdz, zero)
3048 nrtot = max(nrtot + nrtend(i,k)*rdz, zero)
3049 nstot = max(nstot + nstend(i,k)*rdz, zero)
3050
3051!!!! done up to here - Moorthi
3052! calculate melting and freezing of precip
3053! melt snow at +2 C
3054 taux = 1.0
3055
3056 if (.true.) then
3057!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3058 tx1 = t(i,k) + tlat(i,k)*onebcp*deltat
3059 if (tx1 > 275.15_r8) then
3060
3061 if (qstot > zero) then
3062
3063! make sure melting snow doesn't reduce temperature below threshold
3064 dum = xlfocp*qstot*rdzi
3065 if (tx1-dum*deltat < 273.15_r8) then
3066 tx2 = (tx1-275.15_r8) * dti
3067 dum = min(one,max(zero,tx2/dum))
3068 else
3069 dum = one
3070 end if
3071
3072 qric(i,k) = qric(i,k) + dum*qniic(i,k)
3073 nric(i,k) = nric(i,k) + dum*nsic(i,k)
3074 qniic(i,k) = (one-dum)*qniic(i,k)
3075 nsic(i,k) = (one-dum)*nsic(i,k)
3076! heating tendency
3077 tmp = -xlf*dum*qstot*rdzi
3078 meltsdt(i,k) = meltsdt(i,k) + tmp
3079
3080 tlat(i,k) = tlat(i,k) + tmp
3081
3082! if(xlon<0.01.and.xlon>-0.01.and.xlat>1.346
3083! & .and.xlat<1.347.and.k==38)
3084! & write(*,*)"anning_m1",tmp
3085
3086 qrtot = qrtot + dum*qstot
3087 nrtot = nrtot + dum*nstot
3088 qstot = (one-dum)*qstot
3089 nstot = (one-dum)*nstot
3090 preci(i) = (one-dum)*preci(i)
3091 endif
3092 endif
3093
3094 end if
3095 tlataux(i, k) = tlat(i, k)
3096
3097! freeze all rain at -5C for Arctic
3098 if(.true.) then
3099!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3100 tx1 = t(i,k) + tlat(i,k)*onebcp*deltat
3101 if (tx1 < tmelt-five) then
3102
3103 if (qrtot > zero) then
3104
3105! make sure freezing rain doesn't increase temperature above threshold
3106 dum = xlfocp*qrtot*rdzi
3107 if (tx1+dum*deltat > tmelt-five) then
3108 tx2 = -(tx1 - (tmelt-five)) * dti
3109 dum = min(one,max(zero,tx2/dum))
3110 else
3111 dum = one
3112 endif
3113 qniic(i,k) = qniic(i,k) + dum*qric(i,k)
3114 nsic(i,k) = nsic(i,k) + dum*nric(i,k)
3115 qric(i,k) = (one-dum)*qric(i,k)
3116 nric(i,k) = (one-dum)*nric(i,k)
3117! heating tendency
3118 tmp = xlf*dum*qrtot*rdzi
3119 frzrdt(i,k) = frzrdt(i,k) + tmp
3120
3121 tlat(i,k) = tlat(i,k) + tmp
3122
3123! if(xlon<0.01.and.xlon>-0.01.and.xlat>1.346
3124! & .and.xlat<1.347.and.k==38)
3125! & write(*,*)"anning_m2",tmp
3126
3127 qstot = qstot + dum*qrtot
3128 qrtot = (one-dum)*qrtot
3129 nstot = nstot + dum*nrtot
3130 nrtot = (one-dum)*nrtot
3131 preci(i) = preci(i) + dum*(prect(i)-preci(i))
3132 end if
3133 end if
3134 end if
3135
3136!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3137!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3138!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3139
3140
3141
3142! if rain/snow mix ratio is zero so should number concentration
3143
3144 if (qniic(i,k) < qsmall) then
3145 qniic(i,k) = zero
3146 nsic(i,k) = zero
3147 end if
3148
3149 if (qric(i,k) < qsmall) then
3150 qric(i,k) = zero
3151 nric(i,k) = zero
3152 end if
3153
3154! make sure number concentration is a positive number to avoid
3155! taking root of negative
3156
3157 nric(i,k) = max(nric(i,k), zero)
3158 nsic(i,k) = max(nsic(i,k), zero)
3159
3160!.......................................................................
3161! get size distribution parameters for fallspeed calculations
3162!......................................................................
3163! rain
3164
3165 if (qric(i,k) >= qsmall) then
3166 lamr(k) = (pirhow*nric(i,k)/qric(i,k))**oneb3
3167 n0r(k) = nric(i,k)*lamr(k)
3168
3169! check for slope
3170! change lammax and lammin for rain and snow
3171! adjust vars
3172
3173 if (lamr(k) < lamminr) then
3174 lamr(k) = lamminr
3175 tx1 = lamr(k) * lamr(k)
3176 n0r(k) = tx1*tx1*qric(i,k)/pirhow
3177 nric(i,k) = n0r(k)/lamr(k)
3178 else if (lamr(k) > lammaxr) then
3179 lamr(k) = lammaxr
3180 tx1 = lamr(k) * lamr(k)
3181 n0r(k) = tx1*tx1*qric(i,k)/pirhow
3182 nric(i,k) = n0r(k)/lamr(k)
3183 end if
3184
3185
3186! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
3187
3188 tx1 = arn(i,k) / lamr(k)**br
3189 tx2 = 9.1_r8*rhof(i,k)
3190 unr(k) = min(tx1*cons4, tx2)
3191 umr(k) = min(tx1*(cons5/6._r8),tx2)
3192
3193 else
3194 lamr(k) = zero
3195 n0r(k) = zero
3196 umr(k) = zero
3197 unr(k) = zero
3198 end if
3199
3200!calculate mean size of combined rain and snow
3201
3202 if (lamr(k) > zero) then
3203 artmp = n0r(k) * (0.5*pi) / (lamr(k)*lamr(k)*lamr(k))
3204 else
3205 artmp = zero
3206 endif
3207
3208 if (lamc(k) > zero) then
3209 actmp = cdist1(k) * pi * gamma(pgam(k)+three)
3210 & / (four * lamc(k)*lamc(k))
3211 else
3212 actmp = zero
3213 endif
3214
3215 if (actmp > zero .or.artmp > zero) then
3216 rercld(i,k) = rercld(i,k)+three*(qric(i,k)+qcic(i,k))
3217 & /(four*rhow*(actmp+artmp))
3218 arcld(i,k) = arcld(i,k) + one
3219 endif
3220
3221!......................................................................
3222! snow
3223
3224 if (qniic(i,k) >= qsmall) then
3225 lams(k) = (cons6*cs*nsic(i,k) / qniic(i,k))**(one/ds)
3226 n0s(k) = nsic(i,k)*lams(k)
3227
3228! check for slope
3229! adjust vars
3230
3231 if (lams(k) < lammins) then
3232 lams(k) = lammins
3233 n0s(k) = lams(k)**(ds+one)*qniic(i,k)/(cs*cons6)
3234 nsic(i,k) = n0s(k)/lams(k)
3235
3236 else if (lams(k) > lammaxs) then
3237 lams(k) = lammaxs
3238 n0s(k) = lams(k)**(ds+one)*qniic(i,k)/(cs*cons6)
3239 nsic(i,k) = n0s(k)/lams(k)
3240 end if
3241
3242! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
3243
3244 tx1 = asn(i,k) / lams(k)**bs
3245 tx2 = 1.2_r8*rhof(i,k)
3246 ums(k) = min(tx1*(cons8/6._r8), tx2)
3247 uns(k) = min(tx1*cons7, tx2)
3248
3249 else
3250 lams(k) = zero
3251 n0s(k) = zero
3252 ums(k) = zero
3253 uns(k) = zero
3254 end if
3255
3256!c........................................................................
3257! sum over sub-step for average process rates
3258
3259! convert rain/snow q and N for output to history, note,
3260! output is for gridbox average
3261
3262 qrout(i,k) = qrout(i,k) + qric(i,k)*cldmax(i,k)
3263 qsout(i,k) = qsout(i,k) + qniic(i,k)*cldmax(i,k)
3264 tx1 = rho(i,k)*cldmax(i,k)
3265 nrout(i,k) = nrout(i,k) + nric(i,k)*tx1
3266 nsout(i,k) = nsout(i,k) + nsic(i,k)*tx1
3267
3268 end if !fprcp Anning Cheng 9/16/2016
3269 tlat1(i,k) = tlat1(i,k) + tlat(i,k)
3270 tlat1_aux(i,k) = tlat1_aux(i,k) + tlataux(i,k)
3271
3272 qvlat1(i,k) = qvlat1(i,k) + qvlat(i,k)
3273 qctend1(i,k) = qctend1(i,k) + qctend(i,k)
3274 qitend1(i,k) = qitend1(i,k) + qitend(i,k)
3275 nctend1(i,k) = nctend1(i,k) + nctend(i,k)
3276 nitend1(i,k) = nitend1(i,k) + nitend(i,k)
3277
3278 t(i,k) = t(i,k) + tlat(i,k) * deltat/cpp
3279 q(i,k) = q(i,k) + qvlat(i,k) * deltat
3280 qc(i,k) = qc(i,k) + qctend(i,k) * deltat
3281 qi(i,k) = qi(i,k) + qitend(i,k) * deltat
3282 nc(i,k) = nc(i,k) + nctend(i,k) * deltat
3283 ni(i,k) = ni(i,k) + nitend(i,k) * deltat
3284
3285 rainrt1(i,k) = rainrt1(i,k) + rainrt(i,k)
3286
3287!divide rain radius over substeps for average
3288 if (arcld(i,k) > zero) then
3289 rercld(i,k) = rercld(i,k) / arcld(i,k)
3290 end if
3291!calculate precip fluxes and adding them to summing sub-stepping variables
3292
3293 rflx(i,1) = zero
3294 sflx(i,1) = zero
3295
3296
3297 rflx(i,k+1) = qrout(i,k)*rho(i,k)*umr(k)
3298 sflx(i,k+1) = qsout(i,k)*rho(i,k)*ums(k)
3299
3300
3301 rflx1(i,k+1) = rflx1(i,k+1) + rflx(i,k+1)
3302 sflx1(i,k+1) = sflx1(i,k+1) + sflx(i,k+1)
3303
3304!c........................................................................
3305
3306 enddo ! big k loop
3307
3308 prect1(i) = prect1(i) + prect(i)
3309 preci1(i) = preci1(i) + preci(i)
3310
3311! if (lprint) write(0,*)' prect1=',prect1(i),' prect=',
3312! &prect(i),' iter=',iter,' it=',it
3313
3314 enddo !end of big iter loop
3315
3316 do k = 1, pver
3317 rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)
3318
3319 & / max(qcsum_rate1ord(k),1.0e-30_r8)
3320 enddo
3321
3322 endif ! end of if (ltrue(i) == 0) then
3323 enddo ! end of big i loop2
3324
3325! convert dt from sub-step back to full time step
3326
3327 deltat = deltat*real(iter)
3328 dti = one / deltat
3329
3330
3331!.............................................................................
3332
3333 do i=1,ncol !big i loop3
3334
3335 if (ltrue(i) == 0) then ! skip all calculations if no cloud water
3336
3337 do k=1,pver
3338 ! assign default values for effective radius
3339 effc(i,k) = 10._r8
3340 effi(i,k) = 25._r8
3341 effc_fn(i,k) = 10._r8
3342 lamcrad(i,k) = zero
3343 pgamrad(i,k) = zero
3344 deffi(i,k) = zero
3345 end do
3346
3347 else
3348!
3349 nstep = 1 ! initialize nstep for sedimentation sub-steps
3350
3351! divide precip rate by number of sub-steps to get average over time step
3352
3353 prect(i) = prect1(i) * riter
3354 preci(i) = preci1(i) * riter
3355
3356! if (lprint) write(0,*)' prect=',prect(i),' prect1=',prect1(i)
3357! &,' riter=',riter
3358
3359 do k=1,pver
3360
3361! assign variables back to start-of-timestep values before updating after sub-steps
3362
3363 t(i,k) = t1(i,k)
3364 q(i,k) = q1(i,k)
3365 qc(i,k) = qc1(i,k)
3366 qi(i,k) = qi1(i,k)
3367 nc(i,k) = nc1(i,k)
3368 ni(i,k) = ni1(i,k)
3369
3370! divide microphysical tendencies by number of sub-steps to get average over time step
3371
3372 tlat(i,k) = tlat1(i,k) * riter
3373 qvlat(i,k) = qvlat1(i,k) * riter
3374 qctend(i,k) = qctend1(i,k) * riter
3375 qitend(i,k) = qitend1(i,k) * riter
3376 nctend(i,k) = nctend1(i,k) * riter
3377 nitend(i,k) = nitend1(i,k) * riter
3378
3379 tlataux(i,k) = tlat1_aux(i,k) * riter
3380
3381
3382 rainrt(i,k) = rainrt1(i,k) * riter
3383
3384! divide by number of sub-steps to find final values
3385 rflx(i,k+1) = rflx1(i,k+1) * riter
3386 sflx(i,k+1) = sflx1(i,k+1) * riter
3387
3388! divide output precip q and N by number of sub-steps to get average over time step
3389
3390 qrout(i,k) = qrout(i,k) * riter
3391 qsout(i,k) = qsout(i,k) * riter
3392 nrout(i,k) = nrout(i,k) * riter
3393 nsout(i,k) = nsout(i,k) * riter
3394
3395! divide trop_mozart variables by number of sub-steps to get average over time step
3396
3397 nevapr(i,k) = nevapr(i,k) * riter
3398 evapsnow(i,k) = evapsnow(i,k) * riter
3399 prain(i,k) = prain(i,k) * riter
3400 prodsnow(i,k) = prodsnow(i,k) * riter
3401 cmeout(i,k) = cmeout(i,k) * riter
3402
3403 cmeiout(i,k) = cmeiout(i,k) * riter
3404 meltsdt(i,k) = meltsdt(i,k) * riter
3405 frzrdt(i,k) = frzrdt(i,k) * riter
3406
3407
3408! microphysics output
3409 prao(i,k) = prao(i,k) * riter
3410 prco(i,k) = prco(i,k) * riter
3411 mnuccco(i,k) = mnuccco(i,k) * riter
3412 mnuccto(i,k) = mnuccto(i,k) * riter
3413 msacwio(i,k) = msacwio(i,k) * riter
3414 psacwso(i,k) = psacwso(i,k) * riter
3415 bergso(i,k) = bergso(i,k) * riter
3416 bergo(i,k) = bergo(i,k) * riter
3417 prcio(i,k) = prcio(i,k) * riter
3418 praio(i,k) = praio(i,k) * riter
3419
3420 mnuccdo(i,k) = mnuccdo(i,k) * riter
3421 mnuccto(i,k) = mnuccto(i,k) * riter
3422
3423 mnuccro(i,k) = mnuccro(i,k) * riter
3424 pracso(i,k) = pracso(i,k) * riter
3425
3426!!!!DONIFFFF==========================
3427 nnuccdo(i,k) = nnuccdo(i,k) * riter
3428 nnuccco(i,k) = nnuccco(i,k) * riter
3429 nsacwio(i,k) = nsacwio(i,k) * riter
3430 nsubio(i,k) = nsubio(i,k) * riter
3431 nprcio(i,k) = nprcio(i,k) * riter
3432 npraio(i,k) = npraio(i,k) * riter
3433
3434 npccno(i,k) = npccno(i,k) * riter
3435 npsacwso(i,k) = npsacwso(i,k) * riter
3436 nsubco(i,k) = nsubco(i,k) * riter
3437 nprao(i,k) = nprao(i,k) * riter
3438 nprc1o(i,k) = nprc1o(i,k) * riter
3439
3440!=====================================
3441
3442
3443! modify to include snow. in prain & evap (diagnostic here: for wet dep)
3444
3445 prain(i,k) = prain(i,k) + prodsnow(i,k)
3446
3447!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3448! calculate sedimentation for cloud water and ice
3449!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3450
3451! update in-cloud cloud mixing ratio and number concentration
3452! with microphysical tendencies to calculate sedimentation, assign to dummy vars
3453! note: these are in-cloud values***, hence we divide by cloud fraction
3454
3455 tx1 = one / lcldm(i,k)
3456 tx2 = one / icldm(i,k)
3457 dumc(i,k) = (qc(i,k) + qctend(i,k)*deltat) * tx1
3458 dumi(i,k) = (qi(i,k) + qitend(i,k)*deltat) * tx2
3459 dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)*tx1, zero)
3460 dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)*tx2, zero)
3461
3462 if (nccons) then
3463 dumnc(i,k) = ncnst*irho(i,k)
3464 end if
3465
3466 if (nicons) then
3467 dumni(i,k) = ninst*irho(i,k)
3468 end if
3469
3470! obtain new slope parameter to avoid possible singularity
3471
3472 if (dumi(i,k) >= qsmall) then
3473! add upper limit to in-cloud number concentration to prevent numerical error
3474 dumni(i,k) = min(dumni(i,k),dumi(i,k)*1.e20_r8)
3475 lami(k) = (cons1*ci* dumni(i,k)/dumi(i,k))**oneodi
3476
3477! miu_ice(k) = mui_hemp_l(lami(k)) Anning changed here
3478 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
3479 & 10.0_r8), 0.1_r8)
3480 tx1 = one + miu_ice(k)
3481 tx2 = one / gamma(tx1)
3482 aux = (gamma(tx1+di)*tx2) ** oneodi
3483 lami(k) = aux*lami(k)
3484
3485 n0i(k) = niic(i,k)*lami(k)**tx1 * tx2
3486
3487
3488 if (lami(k) < lammini*aux) then
3489 lami(k) = lammini
3490 miu_ice(k) = zero
3491 niic(i,k) = n0i(k)/lami(k)
3492 end if
3493 if (lami(k) > lammaxi*aux) then
3494 lami(k) = lammaxi
3495 miu_ice(k) = zero
3496 niic(i,k) = n0i(k)/lami(k)
3497 end if
3498
3499 else
3500 lami(k) = zero
3501 end if
3502
3503 if (dumc(i,k) >= qsmall) then
3504! add upper limit to in-cloud number concentration to prevent numerical error
3505 dumnc(i,k) = min(dumnc(i,k),dumc(i,k)*1.e20_r8)
3506! add lower limit to in-cloud number concentration
3507 dumnc(i,k) = max(dumnc(i,k),cdnl*irho(i,k))
3508 tx1 = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
3509 pgam(k) = max(two, min(15._r8, one/(tx1*tx1)-one))
3510
3511 lamc(k) = (cr*dumnc(i,k)*gamma(pgam(k)+four)
3512 & / (dumc(i,k)*gamma(pgam(k)+one)))**oneb3
3513 lammin = (pgam(k)+one) / 50.e-6_r8
3514 lammax = (pgam(k)+one) / 2.e-6_r8
3515 lamc(k) = min(lammax, max(lamc(k),lammin))
3516 else
3517 lamc(k) = zero
3518 end if
3519
3520! calculate number and mass weighted fall velocity for droplets
3521! include effects of sub-grid distribution of cloud water
3522
3523
3524 if (dumc(i,k) >= qsmall) then
3525 tx1 = lamc(k) ** bc
3526 unc = acn(i,k)*gamma(one+bc+pgam(k))
3527 & / (tx1*gamma(pgam(k)+one))
3528 umc = acn(i,k)*gamma(four+bc+pgam(k))
3529 & / (tx1*gamma(pgam(k)+four))
3530! fallspeed for output
3531 vtrmc(i,k) = umc
3532 else
3533 umc = zero
3534 unc = zero
3535 end if
3536
3537
3538
3539! calculate number and mass weighted fall velocity for cloud ice
3540
3541 if (dumi(i,k) >= qsmall) then
3542 cons16 = gamma(one+bi+miu_ice(k))/gamma(one+miu_ice(k))
3543 cons17 = gamma(four+bi+miu_ice(k))/gamma(four+miu_ice(k))
3544
3545 tx1 = ain(i,k) / lami(k)**bi
3546 tx2 = 1.2_r8*rhof(i,k)
3547 uni = min(tx1*cons16, tx2)
3548 umi = min(tx1*cons17, tx2)
3549! fallspeed
3550 vtrmi(i,k) = umi
3551 else
3552 umi = zero
3553 uni = zero
3554 end if
3555
3556!DONIFFFF tune up sedimentation vel!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3557 umi = umi*ui_scale
3558 uni = uni*ui_scale
3559!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3560
3561 tx1 = g*rho(i,k)
3562 fi(k) = tx1*umi
3563 fni(k) = tx1*uni
3564 fc(k) = tx1*umc
3565 fnc(k) = tx1*unc
3566
3567! calculate number of split time steps to ensure courant stability criteria
3568! for sedimentation calculations
3569
3570 rgvm = max(fi(k),fc(k),fni(k),fnc(k))
3571 nstep = max(int(rgvm*deltat/pdel(i,k)+one),nstep)
3572
3573! redefine dummy variables - sedimentation is calculated over grid-scale
3574! quantities to ensure conservation
3575
3576 dumc(i,k) = qc(i,k) + qctend(i,k)*deltat
3577 dumi(i,k) = qi(i,k) + qitend(i,k)*deltat
3578 dumnc(i,k) = max(nc(i,k) + nctend(i,k)*deltat, zero)
3579 dumni(i,k) = max(ni(i,k) + nitend(i,k)*deltat, zero)
3580
3581 if (dumc(i,k) < qsmall) dumnc(i,k) = zero
3582 if (dumi(i,k) < qsmall) dumni(i,k) = zero
3583
3584 end do ! end of k loop
3585
3586 do n = 1,nstep
3587
3588 do k = 1,pver
3589 falouti(k) = max(fi(k) * dumi(i,k), zero)
3590 faloutni(k) = max(fni(k) * dumni(i,k), zero)
3591 faloutc(k) = max(fc(k) * dumc(i,k), zero)
3592 faloutnc(k) = max(fnc(k) * dumnc(i,k), zero)
3593 end do
3594
3595! top of model
3596
3597 k = 1
3598 tx1 = one / pdel(i,k)
3599 faltndi = falouti(k) * tx1
3600 faltndni = faloutni(k) * tx1
3601 faltndc = faloutc(k) * tx1
3602 faltndnc = faloutnc(k) * tx1
3603
3604! add fallout terms to microphysical tendencies
3605
3606 tx2 = one / float(nstep)
3607 tx3 = deltat*tx2
3608
3609 qitend(i,k) = qitend(i,k) - faltndi * tx2
3610 nitend(i,k) = nitend(i,k) - faltndni * tx2
3611 qctend(i,k) = qctend(i,k) - faltndc * tx2
3612 nctend(i,k) = nctend(i,k) - faltndnc * tx2
3613
3614! sedimentation tendencies for output
3615 qcsedten(i,k) = qcsedten(i,k) - faltndc * tx2
3616 qisedten(i,k) = qisedten(i,k) - faltndi * tx2
3617
3618 dumi(i,k) = dumi(i,k) - faltndi *tx3
3619 dumni(i,k) = dumni(i,k) - faltndni*tx3
3620 dumc(i,k) = dumc(i,k) - faltndc *tx3
3621 dumnc(i,k) = dumnc(i,k) - faltndnc*tx3
3622
3623 do k = 2,pver
3624
3625! for cloud liquid and ice, if cloud fraction increases with height
3626! then add flux from above to both vapor and cloud water of current level
3627! this means that flux entering clear portion of cell from above evaporates
3628! instantly
3629 tx4 = tx1
3630 tx1 = one / pdel(i,k)
3631
3632 if (lcldm(i,k-1) > zero) then
3633 dum = min(one, lcldm(i,k)/lcldm(i,k-1))
3634 else
3635 dum = min(one, lcldm(i,k))
3636 endif
3637 if (icldm(i,k-1) > zero) then
3638 dum1 = min(one, icldm(i,k)/icldm(i,k-1))
3639 else
3640 dum1 = min(one, icldm(i,k))
3641 endif
3642
3643 faltndqie = (falouti(k) - falouti(k-1)) * tx1
3644 faltndi = (falouti(k) - dum1*falouti(k-1)) * tx1
3645 faltndni = (faloutni(k) - dum1*faloutni(k-1)) * tx1
3646
3647 faltndqce = (faloutc(k) - faloutc(k-1)) * tx1
3648 faltndc = (faloutc(k) - dum*faloutc(k-1)) * tx1
3649 faltndnc = (faloutnc(k)- dum*faloutnc(k-1)) * tx1
3650
3651! add fallout terms to eulerian tendencies
3652
3653 qitend(i,k) = qitend(i,k) - faltndi * tx2
3654 nitend(i,k) = nitend(i,k) - faltndni * tx2
3655 qctend(i,k) = qctend(i,k) - faltndc * tx2
3656 nctend(i,k) = nctend(i,k) - faltndnc * tx2
3657
3658
3659! sedimentation tendencies for output
3660 qcsedten(i,k) = qcsedten(i,k) - faltndc * tx2
3661 qisedten(i,k) = qisedten(i,k) - faltndi * tx2
3662
3663! add terms to to evap/sub of cloud water
3664
3665 qvlat(i,k) = qvlat(i,k) - (faltndqie-faltndi) * tx2
3666
3667! if (lprint .and. k == 29) write(0,*)' qvlatb=',qvlat(i,k),
3668! &' tx2=',tx2,' faltndqie=',faltndqie,' faltndi=',faltndi
3669! for output
3670 qisevap(i,k) = qisevap(i,k) + (faltndqie-faltndi) * tx2
3671
3672
3673 qvlat(i,k) = qvlat(i,k) - (faltndqce-faltndc) * tx2
3674! if (lprint .and. k == 29) write(0,*)' qvlatc=',qvlat(i,k),
3675! &' tx2=',tx2,' faltndqce=',faltndqce,' faltndc=',faltndc
3676
3677 qcsevap(i,k) = qcsevap(i,k) + (faltndqce-faltndc) * tx2
3678
3679
3680 tlat(i,k) = tlat(i,k) + ((faltndqie-faltndi)*xxls
3681 & + (faltndqce-faltndc)*xxlv) * tx2
3682
3683
3684 dumi(i,k) = dumi(i,k) - faltndi * tx3
3685 dumni(i,k) = dumni(i,k) - faltndni * tx3
3686 dumc(i,k) = dumc(i,k) - faltndc * tx3
3687 dumnc(i,k) = dumnc(i,k) - faltndnc * tx3
3688
3689 fni(k) = max(fni(k)*tx1, fni(k-1)*tx4) * pdel(i,k)
3690 fi(k) = max(fi(k)*tx1, fi(k-1)*tx4) * pdel(i,k)
3691 fnc(k) = max(fnc(k)*tx1, fnc(k-1)*tx4) * pdel(i,k)
3692 fc(k) = max(fc(k)*tx1, fc(k-1)*tx4) * pdel(i,k)
3693
3694 end do ! end of k loop
3695
3696! units below are m/s
3697! cloud water/ice sedimentation flux at surface
3698! is added to precip flux at surface to get total precip (cloud + precip water)
3699! rate
3700
3701! if (lprint) write(0,*)' befend prect=',prect(i),' preci=',preci(i)
3702 tx3 = tx2 / (g*1000._r8)
3703 prect(i) = prect(i) + (faloutc(pver)+falouti(pver)) * tx3
3704 preci(i) = preci(i) + falouti(pver) * tx3
3705
3706! if (lprint) write(0,*)' end prect=',prect(i),' preci=',preci(i)
3707
3708 end do ! end of n loop
3709
3710
3711! end sedimentation
3712!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3713! calculate sedimentation for rain and snow
3714! Anning Cheng 9/19/2016, forecast rain and snow
3715! reuse dummy variable for cloud water and ice
3716! iter =1 for fprcp >= 1
3717!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3718! initialize nstep for sedimentation sub-steps
3719! reuse dumc, dumi, dumnc, and dumni
3720 if (fprcp == 1) then
3721 nstep = 1
3722
3723 do k=1,pver
3724
3725! tx1 = deltat / lcldm(i,k)
3726! tx2 = deltat / icldm(i,k)
3727! dumc(i,k) = max(qrn(i,k)+qrtend(i,k)*tx1, zero)
3728! dumi(i,k) = max(qsnw(i,k)+qnitend(i,k)*tx2,zero)
3729! dumnc(i,k) = max(nrn(i,k)+nrtend(i,k)*tx1, zero)
3730! dumni(i,k) = max(nsnw(i,k)+nstend(i,k)*tx2, zero)
3731
3732 dumc(i,k) = max(qrn(i,k)+qrtend(i,k)*deltat, zero)
3733 dumi(i,k) = max(qsnw(i,k)+qnitend(i,k)*deltat,zero)
3734 dumnc(i,k) = max(nrn(i,k)+nrtend(i,k)*deltat, zero)
3735 dumni(i,k) = max(nsnw(i,k)+nstend(i,k)*deltat, zero)
3736
3737! if rain/snow mix ratio is zero so should number concentration
3738
3739 if (dumi(i,k) < qsmall) then
3740 dumi(i,k) = zero
3741 dumni(i,k) = zero
3742 endif
3743
3744 if (dumc(i,k) < qsmall) then
3745 dumc(i,k) = zero
3746 dumnc(i,k) = zero
3747 endif
3748
3749! make sure number concentration is a positive number to avoid
3750! taking root of negative
3751
3752 dumnc(i,k) = max(dumnc(i,k),zero)
3753 dumni(i,k) = max(dumni(i,k),zero)
3754
3755!.......................................................................
3756! get size distribution parameters for fallspeed calculations
3757!......................................................................
3758! rain
3759
3760 if (dumc(i,k) >= qsmall) then
3761 lamr(k) = (pi*rhow*dumnc(i,k)/dumc(i,k))**oneb3
3762 n0r(k) = dumnc(i,k)*lamr(k)
3763
3764! check for slope
3765! change lammax and lammin for rain and snow
3766! adjust vars
3767
3768 if (lamr(k) < lamminr) then
3769
3770 lamr(k) = lamminr
3771
3772 n0r(k) = lamr(k)**4*dumc(i,k)/(pi*rhow)
3773 dumnc(i,k) = n0r(k)/lamr(k)
3774 else if (lamr(k) > lammaxr) then
3775 lamr(k) = lammaxr
3776 n0r(k) = lamr(k)**4*dumc(i,k)/(pi*rhow)
3777 dumnc(i,k) = n0r(k)/lamr(k)
3778 end if
3779
3780
3781! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
3782
3783 tx1 = arn(i,k) / lamr(k)**br
3784 tx2 = 9.1_r8*rhof(i,k)
3785 unr(k) = min(tx1*cons4, tx2)
3786 umr(k) = min(tx1*(cons5/6._r8),tx2)
3787
3788 else
3789 lamr(k) = zero
3790 n0r(k) = zero
3791 umr(k) = zero
3792 unr(k) = zero
3793 end if
3794
3795!......................................................................
3796! snow
3797
3798 if (dumi(i,k) >= qsmall) then
3799 lams(k) = (cons6*cs*dumni(i,k)/ dumi(i,k))**(one/ds)
3800 n0s(k) = dumni(i,k)*lams(k)
3801
3802! check for slope
3803! adjust vars
3804
3805 if (lams(k) < lammins) then
3806 lams(k) = lammins
3807 n0s(k) = lams(k)**(ds+one)*dumi(i,k)/(cs*cons6)
3808 dumni(i,k) = n0s(k)/lams(k)
3809 else if (lams(k) > lammaxs) then
3810 lams(k) = lammaxs
3811 n0s(k) = lams(k)**(ds+one)*dumi(i,k)/(cs*cons6)
3812 dumni(i,k) = n0s(k)/lams(k)
3813 end if
3814
3815! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
3816
3817 tx1 = asn(i,k) / lams(k)**bs
3818 tx2 = 1.2_r8*rhof(i,k)
3819 ums(k) = min(tx1*(cons8/6._r8), tx2)
3820 uns(k) = min(tx1*cons7, tx2)
3821
3822 else
3823 lams(k) = zero
3824 n0s(k) = zero
3825 ums(k) = zero
3826 uns(k) = zero
3827 end if
3828
3829 tx1 = g*rho(i,k)
3830 fi(k) = tx1*ums(k)
3831 fni(k) = tx1*uns(k)
3832 fc(k) = tx1*umr(k)
3833 fnc(k) = tx1*unr(k)
3834
3835! calculate number of split time steps to ensure courant stability criteria
3836! for sedimentation calculations
3837
3838 rgvm = max(fi(k),fc(k),fni(k),fnc(k))
3839 nstep = max(int(rgvm*deltat/pdel(i,k)+one),nstep)
3840
3841! redefine dummy variables - sedimentation is calculated over grid-scale
3842! quantities to ensure conservation
3843
3844 qrn(i,k) = (qrn(i,k) + qrtend(i,k)*deltat)
3845 qsnw(i,k) = (qsnw(i,k) + qnitend(i,k)*deltat)
3846 nrn(i,k) = max((nrn(i,k) + nrtend(i,k)*deltat),zero)
3847 nsnw(i,k) = max((nsnw(i,k) + nstend(i,k)*deltat),zero)
3848
3849 if (qrn(i,k) < qsmall) nrn(i,k) = zero
3850 if (qsnw(i,k) < qsmall) nsnw(i,k) = zero
3851
3852 enddo ! end of k loop
3853
3854 tx2 = one / float(nstep)
3855 tx3 = deltat*tx2
3856 do n = 1,nstep
3857
3858 do k = 1,pver
3859 falouti(k) = max(fi(k) * qsnw(i,k), zero)
3860 faloutni(k) = max(fni(k) * nsnw(i,k), zero)
3861 faloutc(k) = max(fc(k) * qrn(i,k), zero)
3862 faloutnc(k) = max(fnc(k) * nrn(i,k), zero)
3863 end do
3864
3865! top of model
3866
3867 k = 1
3868 tx1 = one / pdel(i,k)
3869 faltndi = falouti(k) * tx1
3870 faltndni = faloutni(k) * tx1
3871 faltndc = faloutc(k) * tx1
3872 faltndnc = faloutnc(k) * tx1
3873
3874! add fallout terms to microphysical tendencies
3875
3876! qnitend(i,k) = qnitend(i,k) - faltndi * tx2
3877! nstend(i,k) = nstend(i,k) - faltndni * tx2
3878! qrtend(i,k) = qrtend(i,k) - faltndc * tx2
3879! nrtend(i,k) = nrtend(i,k) - faltndnc * tx2
3880
3881! sedimentation tendencies for output
3882
3883 qsnw(i,k) = qsnw(i,k) - faltndi * tx3
3884 nsnw(i,k) = nsnw(i,k) - faltndni * tx3
3885 qrn(i,k) = qrn(i,k) - faltndc * tx3
3886 nrn(i,k) = nrn(i,k) - faltndnc * tx3
3887
3888 do k = 2,pver
3889
3890! for rain and snow
3891 tx1 = one / pdel(i,k)
3892! dum = min(one, lcldm(i,k)/lcldm(i,k-1))
3893! dum1 = min(one, icldm(i,k)/icldm(i,k-1))
3894
3895 faltndc = (faloutc(k) - faloutc(k-1)) * tx1
3896 faltndnc = (faloutnc(k) - faloutnc(k-1)) * tx1
3897
3898 faltndi = (falouti(k) - falouti(k-1)) * tx1
3899 faltndni = (faloutni(k) - faloutni(k-1)) * tx1
3900
3901 qsnw(i,k) = qsnw(i,k) - faltndi * tx3
3902 nsnw(i,k) = nsnw(i,k) - faltndni * tx3
3903 qrn(i,k) = qrn(i,k) - faltndc * tx3
3904 nrn(i,k) = nrn(i,k) - faltndnc * tx3
3905
3906 end do
3907
3908
3909 tx5 = tx2 / (g*1000._r8)
3910 prect(i) = prect(i) + (faloutc(pver)+falouti(pver)) * tx5
3911 preci(i) = preci(i) + (falouti(pver)) * tx5
3912
3913 end do ! end of n loop
3914 end if !fprcp ==1
3915! end sedimentation for rain and snow
3916!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3917! get new update for variables that includes sedimentation tendency
3918! note : here dum variables are grid-average, NOT in-cloud
3919
3920! DONE UP TO HERE
3921 do k=1,pver
3922 if (fprcp == 1) then
3923! calculate melting and freezing of precip
3924! melt snow at +2 C
3925
3926!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3927 tx1 = t(i,k) + tlat(i,k)*onebcp*deltat
3928 if (tx1 > 275.15_r8) then
3929 if (qsnw(i,k) > zero) then
3930
3931! make sure melting snow doesn't reduce temperature below threshold
3932 dum = xlfocp*qsnw(i,k)
3933
3934 if (tx1-dum < 273.15_r8) then
3935 tx2 = tx1-275.15_r8
3936 dum = min(one,max(zero,tx2/dum))
3937 else
3938 dum = one
3939 end if
3940
3941 qrn(i,k) = qrn(i,k) + dum*qsnw(i,k)
3942 nrn(i,k) = nrn(i,k) + dum*nsnw(i,k)
3943 qsnw(i,k) = (one-dum)*qsnw(i,k)
3944 nsnw(i,k) = (one-dum)*nsnw(i,k)
3945! heating tendency
3946 tmp = -xlf*dum*qsnw(i,k)/deltat
3947
3948 tlat(i,k) = tlat(i,k) + tmp
3949
3950 preci(i) = (one-dum)*preci(i)
3951 end if
3952 end if
3953
3954! freeze all rain at -5C for Arctic
3955!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3956 tx1 = t(i,k) + tlat(i,k)*onebcp*deltat
3957 if (tx1 < tmelt-five) then
3958
3959 if (qrn(i,k) > zero) then
3960
3961! make sure freezing rain doesn't increase temperature above threshold
3962 dum = xlfocp*qrn(i,k)
3963 if (tx1+dum > tmelt-five) then
3964 tx2 = -(tx1-(tmelt-5._r8))
3965 dum = min(one,max(zero,tx2/dum))
3966 else
3967 dum = one
3968 end if
3969
3970 qsnw(i,k) = qsnw(i,k) + dum*qrn(i,k)
3971 nsnw(i,k) = nsnw(i,k) + dum*nrn(i,k)
3972 qrn(i,k) = (one-dum)*qrn(i,k)
3973 nrn(i,k) = (one-dum)*nrn(i,k)
3974! heating tendency
3975 tmp = xlf*dum*qrn(i,k)/deltat
3976
3977 tlat(i,k) = tlat(i,k) + tmp
3978
3979 preci(i) = preci(i) + dum*(prect(i)-preci(i))
3980 end if
3981 end if
3982
3983! if rain/snow mix ratio is zero so should number concentration
3984
3985 if (qsnw(i,k) < qsmall) then
3986 nsnw(i,k) = zero
3987 end if
3988
3989 if (qrn(i,k) < qsmall) then
3990 nrn(i,k) = zero
3991 end if
3992
3993! make sure number concentration is a positive number to avoid
3994! taking root of negative
3995
3996 nrn(i,k) = max(nrn(i,k),zero)
3997 nsnw(i,k) = max(nsnw(i,k),zero)
3998
3999 qrout(i,k) = qrout(i,k) + qrn(i,k)
4000 qsout(i,k) = qsout(i,k) + qsnw(i,k)
4001 nrout(i,k) = nrout(i,k) + nrn(i,k)*rho(i,k)
4002 nsout(i,k) = nsout(i,k) + nsnw(i,k)*rho(i,k)
4003!.......................................................................
4004 end if !fprcp ==1
4005
4006
4007
4008 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero)
4009 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero)
4010 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero)
4011 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero)
4012
4013
4014 if (nccons) then
4015 dumnc(i,k) = ncnst * irho(i,k) * lcldm(i,k)
4016 endif
4017
4018 if (nicons) then
4019 dumni(i,k) = ninst * irho(i,k) * icldm(i,k)
4020 endif
4021
4022 if (dumc(i,k) < qsmall) dumnc(i,k) = zero
4023 if (dumi(i,k) < qsmall) dumni(i,k) = zero
4024
4025! calculate instantaneous processes (melting, homogeneous freezing)
4026
4027 tx1 = t(i,k) + tlat(i,k)*onebcp*deltat
4028 if (tx1 > tmelt) then
4029 if (dumi(i,k) > zero) then
4030
4031! limit so that melting does not push temperature below freezing
4032 dum = dumi(i,k)*xlfocp
4033 if (tx1-dum < tmelt) then
4034 dum = (tx1-tmelt)*cpoxlf / dumi(i,k)
4035 dum = max(zero, min(one, dum))
4036 else
4037 dum = one
4038 end if
4039
4040 tx2 = dum*dumi(i,k)*dti
4041 qctend(i,k) = qctend(i,k) + tx2
4042! for output
4043 melto(i,k) = tx2
4044
4045! assume melting ice produces droplet
4046! mean volume radius of 8 micron
4047
4048 nctend(i,k) = nctend(i,k) + three*tx2
4049 & / (four*pi*5.12e-16_r8*rhow)
4050
4051 qitend(i,k) = ((one-dum)*dumi(i,k)-qi(i,k)) * dti
4052 nitend(i,k) = ((one-dum)*dumni(i,k)-ni(i,k)) * dti
4053 tlat(i,k) = tlat(i,k) - xlf*tx2
4054
4055 endif
4056 endif
4057
4058! homogeneously freeze droplets between -35 C and -40 C
4059
4060 tx1 = t(i,k) + tlat(i,k)*onebcp*deltat
4061 if (tx1 < 233.15_r8) then
4062
4063 if (dumc(i,k) > zero .and. qc(i,k) > zero) then
4064
4065! limit so that freezing does not push temperature above threshold
4066 dum = dumc(i,k)*xlfocp
4067 if (tx1+dum > 233.15_r8) then
4068 dum = -(tx1-233.15_r8)*cpoxlf / dumc(i,k)
4069 dum = max(zero, min(one, dum))
4070 else
4071 dum = one
4072 end if
4073
4074 tx2 = dum*dumc(i,k)*dti
4075 qitend(i,k) = qitend(i,k) + tx2
4076! for output
4077 homoo(i,k) = tx2
4078
4079! assume 25 micron mean volume radius of homogeneously frozen droplets
4080! consistent with size of detrained ice in stratiform.F90
4081
4082 nitend(i,k) = nitend(i,k) + dum*three*dumc(i,k)
4083 & / (four*3.14_r8*1.563e-14_r8* 500._r8) * dti
4084 qctend(i,k) = ((one-dum)*dumc(i,k)-qc(i,k)) * dti
4085 nctend(i,k) = ((one-dum)*dumnc(i,k)-nc(i,k)) * dti
4086 tlat(i,k) = tlat(i,k) + xlf*tx2
4087
4088 endif
4089 endif
4090
4091! remove any excess over-saturation, which is possible due to non-linearity when adding
4092! together all microphysical processes
4093! follow code similar to old CAM scheme
4094!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4095!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4096! qsn = min(epsqs*esn/(p(i,k)-omeps*esn),one)
4097! if (qtmp > qsn .and. qsn > 0) then
4098! expression below is approximate since there may be ice deposition
4099! dum = (qtmp-qsn)/(one+cons27*qsn/(cpp*rv*ttmp**2))/deltat
4100! add to output cme
4101! now add to tendencies, partition between liquid and ice based on temperature
4102! now add to tendencies, partition between liquid and ice based on te
4103! dum = (qtmp-qsn)/(one+(xxls*dum1+xxlv*(one-dum1))**2 &
4104! *qsn/(cpp*rv*ttmp**2))/deltat
4105!
4106! qctend(i,k)=qctend(i,k)+dum*(one-dum1)
4107! for output
4108! qcreso(i,k)=dum*(one-dum1)
4109! qitend(i,k)=qitend(i,k)+dum*dum1
4110! qireso(i,k)=dum*dum1
4111! qvlat(i,k)=qvlat(i,k)-dum
4112! for output
4113! qvres(i,k)=-dum
4114! tlat(i,k)=tlat(i,k)+dum*(one-dum1)*xxlv+dum*dum1*xxls
4115! end if
4116!end if
4117
4118!!!!!!!!!!!!!!!!!!!!!!!!!!
4119!!!!!!!!!!!!!!!!!!!!!!!!!!!1
4120
4121
4122!...............................................................................
4123! calculate effective radius for pass to radiation code
4124! if no cloud water, default value is 10 micron for droplets,
4125! 25 micron for cloud ice
4126
4127! update cloud variables after instantaneous processes to get effective radius
4128! variables are in-cloud to calculate size dist parameters
4129
4130 tx1 = one / lcldm(i,k)
4131 tx2 = one / icldm(i,k)
4132 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero) * tx1
4133 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero) * tx2
4134 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero) * tx1
4135 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero) * tx2
4136
4137! if ((dumc(i, k)*1e6 .gt. 1.0) .and. dumnc(i, k) .lt. 1e-20) then!
4138
4139! print *, 'dumnc', dumnc(i,k)*1e-6
4140! print *, 'dumc', dumc(i, k)*1e6
4141! print *, i, k
4142
4143! end if
4144
4145 if (nccons) then
4146 dumnc(i,k) = ncnst * irho(i,k)
4147 end if
4148
4149 if (nicons) then
4150 dumni(i,k) = ninst * irho(i,k)
4151 end if
4152
4153
4154! limit in-cloud mixing ratio of water and ice to reasonable value of 5 g kg-1
4155
4156 dumc(i,k) = min(dumc(i,k),5.e-3_r8)
4157 dumi(i,k) = min(dumi(i,k),5.e-3_r8)
4158
4159!...................
4160! cloud ice effective radius
4161
4162 if (dumi(i,k) >= qsmall) then
4163! add upper limit to in-cloud number concentration to prevent numerical error
4164 dumni(i,k) = min(dumni(i,k),dumi(i,k)*1.e20_r8)
4165 lami(k) = (cons1*ci* dumni(i,k)/dumi(i,k))**oneodi
4166
4167! miu_ice(k) = mui_hemp_l(lami(k))
4168 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
4169 & 10.0_r8), 0.1_r8)
4170 tx1 = one + miu_ice(k)
4171 tx2 = one / gamma(tx1)
4172 aux = (gamma(tx1+di)*tx2) ** oneodi
4173 lami(k) = aux*lami(k)
4174 n0i(k) = niic(i,k) * lami(k)**tx1 * tx2
4175
4176 if (lami(k) < lammini*aux) then
4177 miu_ice(k) = zero
4178 lami(k) = lammini
4179 n0i(k) = lami(k)**(di+one)*dumi(i,k)/(ci*cons1)
4180 niic(i,k) = n0i(k)/lami(k)
4181! adjust number conc if needed to keep mean size in reasonable range
4182 nitend(i,k) = (niic(i,k)*icldm(i,k)-ni(i,k))/deltat
4183 else if (lami(k) > lammaxi*aux) then
4184 miu_ice(k) = zero
4185 lami(k) = lammaxi
4186 n0i(k) = lami(k)**(di+one)*dumi(i,k)/(ci*cons1)
4187 niic(i,k) = n0i(k)/lami(k)
4188 nitend(i,k) = (niic(i,k)*icldm(i,k)-ni(i,k))/deltat
4189 end if
4190
4191 effi(i,k) = 1.e6_r8*gamma(four+miu_ice(k))
4192 & / (two*lami(k)*gamma(three+miu_ice(k)))
4193 else
4194 effi(i,k) = 25._r8
4195 end if
4196
4197!...................
4198! cloud droplet effective radius
4199
4200
4201 if (dumc(i,k) >= qsmall) then
4202
4203! add upper limit to in-cloud number concentration to prevent numerical error
4204 dumnc(i,k) = min(dumnc(i,k),dumc(i,k)*1.e20_r8)
4205
4206!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4207! set tendency to ensure minimum droplet concentration
4208! after update by microphysics, except when lambda exceeds bounds on mean drop
4209! size or if there is no cloud water
4210 if (dumnc(i,k) < cdnl*irho(i,k)) then
4211 nctend(i,k) = (cdnl*irho(i,k)*cldm(i,k)-nc(i,k))*dti
4212 end if
4213 dumnc(i,k) = max(dumnc(i,k),cdnl*irho(i,k))
4214
4215 if (nccons) then
4216! make sure nc is consistence with the constant N by adjusting tendency, need
4217! to multiply by cloud fraction
4218! note that nctend may be further adjusted below if mean droplet size is
4219! out of bounds
4220
4221 nctend(i,k) = (ncnst*irho(i,k)*lcldm(i,k)-nc(i,k))
4222 & * dti
4223 end if
4224
4225
4226!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4227 tx1 = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
4228 pgam(k) = max(two, min(15._r8, one/(tx1*tx1)-one))
4229
4230 tx1 = gamma(pgam(k)+one)
4231 tx2 = gamma(pgam(k)+four)
4232 lamc(k) = (cr*dumnc(i,k)*tx2
4233 & / (dumc(i,k)*tx1))**oneb3
4234 lammin = (pgam(k)+one) / 50.e-6_r8
4235 lammax = (pgam(k)+one) / 2.e-6_r8
4236 if (lamc(k) < lammin) then
4237 lamc(k) = lammin
4238 ncic(i,k) = 6._r8*lamc(k)*lamc(k)*lamc(k)*dumc(i,k)
4239 & * tx1 / (pirhow*tx2)
4240! adjust number conc if needed to keep mean size in reasonable range
4241 nctend(i,k) = (ncic(i,k)*lcldm(i,k)-nc(i,k))*dti
4242
4243 else if (lamc(k) > lammax) then
4244 lamc(k) = lammax
4245 ncic(i,k) = 6._r8*lamc(k)*lamc(k)*lamc(k)*dumc(i,k)
4246 & * tx1 / (pirhow*tx2)
4247! adjust number conc if needed to keep mean size in reasonable range
4248 nctend(i,k) = (ncic(i,k)*lcldm(i,k)-nc(i,k))*dti
4249 end if
4250
4251 effc(i,k) = tx2 / (gamma(pgam(k)+three)*lamc(k)*2.e6_r8)
4252!assign output fields for shape here
4253 lamcrad(i,k) = lamc(k)
4254 pgamrad(i,k) = pgam(k)
4255
4256 else
4257 effc(i,k) = 10._r8
4258 lamcrad(i,k) = zero
4259 pgamrad(i,k) = zero
4260 end if
4261
4262
4263! ice effective diameter for david mitchell's optics
4264 deffi(i,k) = effi(i,k) * (rhoi / 917._r8*two)
4265
4266
4267!!! recalculate effective radius for constant number, in order to separate
4268! first and second indirect effects
4269! assume constant number of 10^8 kg-1
4270
4271 dumnc(i,k) = 1.e8
4272
4273 if (dumc(i,k) >= qsmall) then
4274 tx1 = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
4275 pgam(k) = max(two, min(15._r8, one/(tx1*tx1)-one))
4276 tx2 = gamma(pgam(k)+four)
4277 lamc(k) = (cr*dumnc(i,k)*tx2
4278 & / (dumc(i,k)*gamma(pgam(k)+one)))**oneb3
4279 lammin = (pgam(k)+one) / 50.e-6_r8
4280 lammax = (pgam(k)+one) / 2.e-6_r8
4281 if (lamc(k) < lammin) then
4282 lamc(k) = lammin
4283 else if (lamc(k) > lammax) then
4284 lamc(k) = lammax
4285 end if
4286 effc_fn(i,k) = tx2/(gamma(pgam(k)+three)*lamc(k)*2.e6_r8)
4287
4288 else
4289 effc_fn(i,k) = 10._r8
4290 end if
4291
4292
4293!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
4294
4295 end do ! end of k loop after n loop
4296
4297
4298 endif ! end of if (ltrue(i) == 0) loop
4299
4300! convert dt from sub-step back to full time step
4301
4302 deltat = deltat*real(iter)
4303 dti = one / deltat
4304
4305
4306 do k=1,pver
4307! if updated q (after microphysics) is zero, then ensure updated n is also zero
4308
4309 if (qc(i,k)+qctend(i,k)*deltat < qsmall)
4310 & nctend(i,k) = -nc(i,k) * dti
4311 if (qi(i,k)+qitend(i,k)*deltat < qsmall)
4312 & nitend(i,k) = -ni(i,k) * dti
4313 end do
4314
4315
4316 end do !end big i loop3
4317
4318
4319!print*, 'nctend', 1.0e-6*nctend*deltat
4320
4321! hm add rain/snow mixing ratio and number concentration as diagnostic
4322
4323#ifdef CAM
4324 call outfld('QRAIN',qrout, pcols, lchnk)
4325 call outfld('QSNOW',qsout, pcols, lchnk)
4326 call outfld('NRAIN',nrout, pcols, lchnk)
4327 call outfld('NSNOW',nsout, pcols, lchnk)
4328#endif
4329
4330! add snow output
4331 do k=1,pver
4332 do i = 1,ncol
4333 if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero) then
4334 dsout(i,k) = (pirhosn * nsout(i,k)/qsout(i,k)) **(-oneb3)
4335 endif
4336 end do
4337 end do
4338
4339#ifdef CAM
4340 call outfld('DSNOW',dsout, pcols, lchnk)
4341
4342! add ip fluxes as output fields
4343 call outfld('MGFLXPRC',rflx, pcols, lchnk)
4344 call outfld('MGFLXSNW',sflx, pcols, lchnk)
4345#endif
4346
4347 do k=1,pver
4348 do i = 1,ncol
4349 if ((qc(i,k)+qctend(i,k)*deltat >= qsmall) .and.
4350 & (cldmax(i,k) > zero) .and. (lcldm(i, k) > zero) .and.
4351 & (nc(i,k)+nctend(i,k)*deltat > 0.0)) then
4352
4353 tx1 = rho(i,k) / lcldm(i,k)
4354 tx2 = (qc(i,k)+qctend(i,k)*deltat)*tx1*1000._r8
4355 dum = tx2 * tx2 * lcldm(i,k)
4356 & / (0.109_r8*(nc(i,k)+nctend(i,k)*deltat)
4357 & *tx1*1.e-6_r8*cldmax(i,k))
4358 else
4359 dum = zero
4360 end if
4361 if (qi(i,k)+qitend(i,k)*deltat >= qsmall) then
4362 dum1 = ((qi(i,k)+qitend(i,k)*deltat)*rho(i,k)/icldm(i,k)
4363 & * 1000._r8/0.1_r8)**(one/0.63_r8)
4364 & * icldm(i,k)/cldmax(i,k)
4365 else
4366 dum1 = zero
4367 end if
4368
4369 if (qsout(i,k) >= qsmall) then
4370 dum1 = dum1 + (qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)
4371 & **(one/0.63_r8)
4372 end if
4373
4374 refl(i,k) = dum + dum1
4375
4376 if (rainrt(i,k) >= 0.001) then
4377 dum = log10(rainrt(i,k)**6._r8) + 16._r8
4378 dum = 10._r8**(dum/10._r8)
4379 else
4380 dum = 0.
4381 end if
4382
4383 refl(i,k) = refl(i,k) + dum
4384
4385 areflz(i,k) = refl(i,k)
4386
4387 if (refl(i,k) > minrefl) then
4388 refl(i,k) = 10._r8*log10(refl(i,k))
4389 else
4390 refl(i,k) = -9999._r8
4391 end if
4392
4393 if (refl(i,k) > mindbz) then
4394 arefl(i,k) = refl(i,k)
4395 frefl(i,k) = one
4396 else
4397 arefl(i,k) = zero
4398 areflz(i,k) = zero
4399 frefl(i,k) = zero
4400 end if
4401
4402 csrfl(i,k) = min(csmax,refl(i,k))
4403
4404 if (csrfl(i,k) > csmin) then
4405 acsrfl(i,k) = refl(i,k)
4406 fcsrfl(i,k) = one
4407 else
4408 acsrfl(i,k) = zero
4409 fcsrfl(i,k) = zero
4410 end if
4411
4412 end do
4413 end do
4414
4415#ifdef CAM
4416 call outfld('REFL',refl, pcols, lchnk)
4417 call outfld('AREFL',arefl, pcols, lchnk)
4418 call outfld('AREFLZ',areflz, pcols, lchnk)
4419 call outfld('FREFL',frefl, pcols, lchnk)
4420 call outfld('CSRFL',csrfl, pcols, lchnk)
4421 call outfld('ACSRFL',acsrfl, pcols, lchnk)
4422 call outfld('FCSRFL',fcsrfl, pcols, lchnk)
4423
4424 call outfld('RERCLD',rercld, pcols, lchnk)
4425#endif
4426
4427! averaging for snow and rain number and diameter
4428
4429 do k=1,pver
4430 do i = 1,ncol
4431 qrout2(i,k) = zero
4432 nrout2(i,k) = zero
4433 drout2(i,k) = zero
4434 freqr(i,k) = zero
4435
4436 qsout2(i,k) = zero
4437 nsout2(i,k) = zero
4438 dsout2(i,k) = zero
4439 freqs(i,k) = zero
4440
4441 if (qrout(i,k) > 1.e-7_r8 .and. nrout(i,k) > zero) then
4442 qrout2(i,k) = qrout(i,k)
4443 nrout2(i,k) = nrout(i,k)
4444 drout2(i,k) = (pirhow*nrout(i,k) / qrout(i,k))**(-oneb3)
4445 freqr(i,k) = one
4446 endif
4447 if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero) then
4448 qsout2(i,k) = qsout(i,k)
4449 nsout2(i,k) = nsout(i,k)
4450 dsout2(i,k) = (pirhosn*nsout(i,k) / qsout(i,k))**(-oneb3)
4451 freqs(i,k) = one
4452 endif
4453
4454! output activated liquid and ice (convert from #/kg -> #/m3)
4455 ncai(i,k) = dum2i(i,k)*rho(i,k)
4456 ncal(i,k) = dum2l(i,k)*rho(i,k)
4457 end do
4458 end do
4459
4460
4461#ifdef CAM
4462 call outfld('NCAL',ncal, pcols,lchnk)
4463 call outfld('NCAI',ncai, pcols,lchnk)
4464
4465!add averaged output fields.
4466 call outfld('AQRAIN',qrout2, pcols,lchnk)
4467 call outfld('AQSNOW',qsout2, pcols,lchnk)
4468 call outfld('ANRAIN',nrout2, pcols,lchnk)
4469 call outfld('ANSNOW',nsout2, pcols,lchnk)
4470 call outfld('ADRAIN',drout2, pcols,lchnk)
4471 call outfld('ADSNOW',dsout2, pcols,lchnk)
4472 call outfld('FREQR',freqr, pcols,lchnk)
4473 call outfld('FREQS',freqs, pcols,lchnk)
4474#endif
4475
4476!redefine fice here....
4477 do k=1,pver
4478 do i=1,ncol
4479 nfice(i,k) = zero
4480 tx1 = qc(i,k) + qctend(i,k)*deltat ! water
4481 tx2 = qi(i,k) + qitend(i,k)*deltat ! ice
4482 dumfice = qsout(i,k) + qrout(i,k) + tx1 + tx2
4483 if (dumfice > zero) then
4484 nfice(i,k) = (qsout(i,k) + tx2) / dumfice
4485 endif
4486 enddo
4487 enddo
4488
4489#ifdef CAM
4490 call outfld('FICE',nfice, pcols, lchnk)
4491#endif
4492
4493 return
4494 end subroutine mmicro_pcond
4495
4496#ifdef CAM
4497
4498!##############################################################################
4502 subroutine findsp1 (lchnk, ncol, q, t, p, tsp, qsp)
4503!-----------------------------------------------------------------------
4504!
4505! Purpose:
4506! find the wet bulb temperature for a given t and q
4507! in a longitude height section
4508! wet bulb temp is the temperature and spec humidity that is
4509! just saturated and has the same enthalpy
4510! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
4511! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
4512!
4513! Method:
4514! a Newton method is used
4515! first guess uses an algorithm provided by John Petch from the UKMO
4516! we exclude points where the physical situation is unrealistic
4517! e.g. where the temperature is outside the range of validity for the
4518! saturation vapor pressure, or where the water vapor pressure
4519! exceeds the ambient pressure, or the saturation specific humidity is
4520! unrealistic
4521!
4522! Author: P. Rasch
4523!
4524!-----------------------------------------------------------------------
4525!
4526! input arguments
4527!
4528 integer, intent(in) :: lchnk
4529 integer, intent(in) :: ncol
4530
4531 real(r8), intent(in) :: q(pcols,pver)
4532 real(r8), intent(in) :: t(pcols,pver)
4533 real(r8), intent(in) :: p(pcols,pver)
4534!
4535! output arguments
4536!
4537 real(r8), intent(out) :: tsp(pcols,pver)
4538 real(r8), intent(out) :: qsp(pcols,pver)
4539!
4540! local variables
4541!
4542 integer i
4543 integer k
4544 logical lflg
4545 integer iter
4546 integer l
4547 logical :: error_found
4548
4549 real(r8) omeps
4550 real(r8) trinv
4551 real(r8) es
4552 real(r8) desdt
4553! real(r8) desdp ! change in sat vap pressure wrt pressure
4554 real(r8) dqsdt
4555 real(r8) dgdt
4556 real(r8) g
4557 real(r8) weight(pcols)
4558 real(r8) hlatsb
4559 real(r8) hlatvp
4560 real(r8) hltalt(pcols,pver)
4561 real(r8) tterm
4562 real(r8) qs
4563 real(r8) tc
4564
4565! work variables
4566 real(r8) t1, q1, dt, dq
4567 real(r8) dtm, dqm
4568 real(r8) qvd, a1, tmp
4569 real(r8) rair
4570 real(r8) r1b, c1, c2, c3
4571 real(r8) denom
4572 real(r8) dttol
4573 real(r8) dqtol
4574 integer doit(pcols)
4575 real(r8) enin(pcols), enout(pcols)
4576 real(r8) tlim(pcols)
4577
4578 omeps = one - epsqs
4579 trinv = one / ttrice
4580 a1 = 7.5_r8*log(10._r8)
4581 rair = 287.04_r8
4582 c3 = rair*a1/cp
4583 dtm = zero
4584 dqm = zero
4585 dttol = 1.e-4_r8
4586 dqtol = 1.e-4_r8
4587! tmin = 173.16 ! the coldest temperature we can deal with
4588!
4589! max number of times to iterate the calculation
4590 iter = 10
4591!
4592 do k = 1,pver
4593
4594!
4595! first guess on the wet bulb temperature
4596!
4597 do i = 1,ncol
4598
4599#ifdef DEBUG
4600 if ( (lchnk == lchnklook(nlook) ) .and. (i ==
4601 & icollook(nlook) ) ) then
4602 write (iulog,*) ' '
4603 write (iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
4604 endif
4605#endif
4606
4607! limit the temperature range to that relevant to the sat vap pres tables
4608#if ( defined WACCM_PHYS )
4609 tlim(i) = min(max(t(i,k),128._r8),373._r8)
4610#else
4611 tlim(i) = min(max(t(i,k),173._r8),373._r8)
4612#endif
4613 es = estblf(tlim(i))
4614 denom = p(i,k) - omeps*es
4615 qs = epsqs*es/denom
4616 doit(i) = 0
4617 enout(i) = one
4618! make sure a meaningful calculation is possible
4619 if (p(i,k) > 5._r8*es .and. qs > zero .and. qs < 0.5_r8) then
4620!
4621! Saturation specific humidity
4622!
4623 qs = min(epsqs*es/denom,one)
4624!
4625! "generalized" analytic expression for t derivative of es
4626! accurate to within 1 percent for 173.16 < t < 373.16
4627!
4628! Weighting of hlat accounts for transition from water to ice
4629! polynomial expression approximates difference between es over
4630! water and es over ice from 0 to -ttrice (C) (min of ttrice is
4631! -40): required for accurate estimate of es derivative in transition
4632! range from ice to water also accounting for change of hlatv with t
4633! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
4634!
4635 tc = tlim(i) - tt0
4636 lflg = (tc >= -ttrice .and. tc < zero)
4637 weight(i) = min(-tc*trinv,one)
4638 hlatsb = hlatv + weight(i)*hlatf
4639 hlatvp = hlatv - 2369.0_r8*tc
4640 if (tlim(i) < tt0) then
4641 hltalt(i,k) = hlatsb
4642 else
4643 hltalt(i,k) = hlatvp
4644 end if
4645 enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)
4646
4647! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
4648 tmp = q(i,k) - qs
4649 c1 = hltalt(i,k)*c3
4650 c2 = (tlim(i) + 36._r8)**2
4651 r1b = c2/(c2 + c1*qs)
4652 qvd = r1b*tmp
4653 tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
4654#ifdef DEBUG
4655 if ( (lchnk == lchnklook(nlook) ) .and. (i ==
4656 & icollook(nlook) ) ) then
4657 write (iulog,*) ' relative humidity ', q(i,k)/qs
4658 write (iulog,*) ' first guess ', tsp(i,k)
4659 endif
4660#endif
4661 es = estblf(tsp(i,k))
4662 qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),one)
4663 else
4664 doit(i) = 1
4665 tsp(i,k) = tlim(i)
4666 qsp(i,k) = q(i,k)
4667 enin(i) = one
4668 endif
4669 end do
4670!
4671! now iterate on first guess
4672!
4673 do l = 1, iter
4674 dtm = 0
4675 dqm = 0
4676 do i = 1,ncol
4677 if (doit(i) == 0) then
4678 es = estblf(tsp(i,k))
4679!
4680! Saturation specific humidity
4681!
4682 qs = min(epsqs*es/(p(i,k) - omeps*es),one)
4683!
4684! "generalized" analytic expression for t derivative of es
4685! accurate to within 1 percent for 173.16 < t < 373.16
4686!
4687! Weighting of hlat accounts for transition from water to ice
4688! polynomial expression approximates difference between es over
4689! water and es over ice from 0 to -ttrice (C) (min of ttrice is
4690! -40): required for accurate estimate of es derivative in transition
4691! range from ice to water also accounting for change of hlatv with t
4692! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
4693!
4694 tc = tsp(i,k) - tt0
4695 lflg = (tc >= -ttrice .and. tc < zero)
4696 weight(i) = min(-tc*trinv,one)
4697 hlatsb = hlatv + weight(i)*hlatf
4698 hlatvp = hlatv - 2369.0_r8*tc
4699 if (tsp(i,k) < tt0) then
4700 hltalt(i,k) = hlatsb
4701 else
4702 hltalt(i,k) = hlatvp
4703 end if
4704 if (lflg) then
4705 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*
4706 &pcf(5))))
4707 else
4708 tterm = zero
4709 end if
4710 desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
4711 dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
4712! g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
4713 g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
4714 dgdt = -(cp + hltalt(i,k)*dqsdt)
4715 t1 = tsp(i,k) - g/dgdt
4716 dt = abs(t1 - tsp(i,k))/t1
4717 tsp(i,k) = max(t1,tmin)
4718 es = estblf(tsp(i,k))
4719 q1 = min(epsqs*es/(p(i,k) - omeps*es),one)
4720 dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
4721 qsp(i,k) = q1
4722#ifdef DEBUG
4723 if ( (lchnk == lchnklook(nlook) ) .and. (i ==
4724 & icollook(nlook) ) ) then
4725 write (iulog,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g
4726 endif
4727#endif
4728 dtm = max(dtm,dt)
4729 dqm = max(dqm,dq)
4730! if converged at this point, exclude it from more iterations
4731 if (dt < dttol .and. dq < dqtol) then
4732 doit(i) = 2
4733 endif
4734 enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
4735! bail out if we are too near the end of temp range
4736#if ( defined WACCM_PHYS )
4737 if (tsp(i,k) < 130.16_r8) then
4738#else
4739 if (tsp(i,k) < 174.16_r8) then
4740#endif
4741 doit(i) = 4
4742 endif
4743 else
4744 endif
4745 end do
4746
4747 if (dtm < dttol .and. dqm < dqtol) then
4748 go to 10
4749 endif
4750
4751 end do
4752 10 continue
4753
4754 error_found = .false.
4755 if (dtm > dttol .or. dqm > dqtol) then
4756 do i = 1,ncol
4757 if (doit(i) == 0) error_found = .true.
4758 end do
4759 if (error_found) then
4760 do i = 1,ncol
4761 if (doit(i) == 0) then
4762 write (iulog,*) ' findsp not converging at point i, k ', i, k
4763 write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k),
4764 & enin(i)
4765 write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k),
4766 & enout(i)
4767 call endrun ('FINDSP')
4768 endif
4769 end do
4770 endif
4771 endif
4772 do i = 1,ncol
4773 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+
4774 &enout(i))) > 1.e-4_r8) then
4775 error_found = .true.
4776 endif
4777 end do
4778 if (error_found) then
4779 do i = 1,ncol
4780 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+
4781 &enout(i))) > 1.e-4_r8) then
4782 write (iulog,*) ' the enthalpy is not conserved for point ', i,
4783 & k, enin(i), enout(i)
4784 write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k),
4785 & enin(i)
4786 write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k),
4787 & enout(i)
4788 call endrun ('FINDSP')
4789 endif
4790 end do
4791 endif
4792
4793 end do
4794
4795 return
4796 end subroutine findsp1
4797
4801 subroutine findsp1_water (lchnk, ncol, q, t, p, tsp, qsp)
4802!-----------------------------------------------------------------------
4803!
4804! Purpose:
4805! find the wet bulb temperature for a given t and q
4806! in a longitude height section
4807! wet bulb temp is the temperature and spec humidity that is
4808! just saturated and has the same enthalpy
4809! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
4810! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
4811!
4812! Method:
4813! a Newton method is used
4814! first guess uses an algorithm provided by John Petch from the UKMO
4815! we exclude points where the physical situation is unrealistic
4816! e.g. where the temperature is outside the range of validity for the
4817! saturation vapor pressure, or where the water vapor pressure
4818! exceeds the ambient pressure, or the saturation specific humidity is
4819! unrealistic
4820!
4821! Author: P. Rasch
4822!
4823!-----------------------------------------------------------------------
4824!
4825! input arguments
4826!
4827 integer, intent(in) :: lchnk
4828 integer, intent(in) :: ncol
4829
4830 real(r8), intent(in) :: q(pcols,pver)
4831 real(r8), intent(in) :: t(pcols,pver)
4832 real(r8), intent(in) :: p(pcols,pver)
4833!
4834! output arguments
4835!
4836 real(r8), intent(out) :: tsp(pcols,pver)
4837 real(r8), intent(out) :: qsp(pcols,pver)
4838!
4839! local variables
4840!
4841 integer i
4842 integer k
4843 logical lflg
4844 integer iter
4845 integer l
4846 logical :: error_found
4847
4848 real(r8) omeps
4849 real(r8) trinv
4850 real(r8) es
4851 real(r8) desdt
4852! real(r8) desdp ! change in sat vap pressure wrt pressure
4853 real(r8) dqsdt
4854 real(r8) dgdt
4855 real(r8) g
4856 real(r8) weight(pcols)
4857 real(r8) hlatsb
4858 real(r8) hlatvp
4859 real(r8) hltalt(pcols,pver)
4860 real(r8) tterm
4861 real(r8) qs
4862 real(r8) tc
4863
4864! work variables
4865 real(r8) t1, q1, dt, dq
4866 real(r8) dtm, dqm
4867 real(r8) qvd, a1, tmp
4868 real(r8) rair
4869 real(r8) r1b, c1, c2, c3
4870 real(r8) denom
4871 real(r8) dttol
4872 real(r8) dqtol
4873 integer doit(pcols)
4874 real(r8) enin(pcols), enout(pcols)
4875 real(r8) tlim(pcols)
4876
4877 omeps = one - epsqs
4878 a1 = 7.5_r8*log(10._r8)
4879 rair = 287.04_r8
4880 c3 = rair*a1/cp
4881 dtm = zero
4882 dqm = zero
4883 dttol = 1.e-4_r8
4884 dqtol = 1.e-4_r8
4885!
4886! max number of times to iterate the calculation
4887 iter = 8
4888!
4889 do k = 1,pver
4890
4891!
4892! first guess on the wet bulb temperature
4893!
4894 do i = 1,ncol
4895
4896#ifdef DEBUG
4897 if ( (lchnk == lchnklook(nlook) ) .and. (i ==
4898 & icollook(nlook) ) ) then
4899 write (iulog,*) ' '
4900 write (iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
4901 endif
4902#endif
4903! limit the temperature range to that relevant to the sat vap pres tables
4904#if ( defined WACCM_PHYS )
4905 tlim(i) = min(max(t(i,k),128._r8),373._r8)
4906#else
4907 tlim(i) = min(max(t(i,k),173._r8),373._r8)
4908#endif
4909
4910 es = polysvp(tlim(i),0)
4911 denom = p(i,k) - omeps*es
4912 qs = epsqs*es/denom
4913 doit(i) = 0
4914 enout(i) = one
4915! make sure a meaningful calculation is possible
4916 if (p(i,k) > 5._r8*es .and. qs > zero .and. qs < 0.5_r8) then
4917!
4918! Saturation specific humidity
4919!
4920 qs = min(epsqs*es/denom,one)
4921!
4922! "generalized" analytic expression for t derivative of es
4923! accurate to within 1 percent for 173.16 < t < 373.16
4924
4925!
4926! No icephs or water to ice transition
4927!
4928 hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
4929 hlatsb = hlatv
4930 if (tlim(i) < tt0) then
4931 hltalt(i,k) = hlatsb
4932 else
4933 hltalt(i,k) = hlatvp
4934 end if
4935!--xl
4936 enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)
4937
4938! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
4939 tmp = q(i,k) - qs
4940 c1 = hltalt(i,k)*c3
4941 c2 = (tlim(i) + 36._r8)**2
4942 r1b = c2/(c2 + c1*qs)
4943 qvd = r1b*tmp
4944 tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
4945#ifdef DEBUG
4946 if ( (lchnk == lchnklook(nlook) ) .and. (i ==
4947 & icollook(nlook) ) ) then
4948 write (iulog,*) ' relative humidity ', q(i,k)/qs
4949 write (iulog,*) ' first guess ', tsp(i,k)
4950 endif
4951#endif
4952 es = polysvp(tsp(i,k),0)
4953 qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),one)
4954 else
4955 doit(i) = 1
4956 tsp(i,k) = tlim(i)
4957 qsp(i,k) = q(i,k)
4958 enin(i) = one
4959 endif
4960 end do
4961!
4962! now iterate on first guess
4963!
4964 do l = 1, iter
4965 dtm = 0
4966 dqm = 0
4967 do i = 1,ncol
4968 if (doit(i) == 0) then
4969 es = polysvp(tsp(i,k),0)
4970!
4971! Saturation specific humidity
4972!
4973 qs = min(epsqs*es/(p(i,k) - omeps*es),one)
4974!
4975! "generalized" analytic expression for t derivative of es
4976! accurate to within 1 percent for 173.16 < t < 373.16
4977!
4978!
4979! No icephs or water to ice transition
4980!
4981 hlatvp = hlatv - 2369.0*(tsp(i,k)-tt0)
4982 hlatsb = hlatv
4983 if (tsp(i,k) < tt0) then
4984 hltalt(i,k) = hlatsb
4985 else
4986 hltalt(i,k) = hlatvp
4987 end if
4988 desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k))
4989 dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
4990 g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
4991 dgdt = -(cp + hltalt(i,k)*dqsdt)
4992 t1 = tsp(i,k) - g/dgdt
4993 dt = abs(t1 - tsp(i,k))/t1
4994 tsp(i,k) = max(t1,tmin)
4995
4996 es = polysvp(tsp(i,k),0)
4997 q1 = min(epsqs*es/(p(i,k) - omeps*es),one)
4998 dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
4999 qsp(i,k) = q1
5000#ifdef DEBUG
5001 if ( (lchnk == lchnklook(nlook) ) .and. (i ==
5002 & icollook(nlook) ) ) then
5003 write (iulog,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g
5004 endif
5005#endif
5006 dtm = max(dtm,dt)
5007 dqm = max(dqm,dq)
5008! if converged at this point, exclude it from more iterations
5009 if (dt < dttol .and. dq < dqtol) then
5010 doit(i) = 2
5011 endif
5012 enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
5013! bail out if we are too near the end of temp range
5014#if ( defined WACCM_PHYS )
5015 if (tsp(i,k) < 130.16_r8) then
5016#else
5017 if (tsp(i,k) < 174.16_r8) then
5018#endif
5019 doit(i) = 4
5020 endif
5021 else
5022 endif
5023 end do
5024
5025 if (dtm < dttol .and. dqm < dqtol) then
5026 go to 10
5027 endif
5028
5029 end do
5030 10 continue
5031
5032 error_found = .false.
5033 if (dtm > dttol .or. dqm > dqtol) then
5034 do i = 1,ncol
5035 if (doit(i) == 0) error_found = .true.
5036 end do
5037 if (error_found) then
5038 do i = 1,ncol
5039 if (doit(i) == 0) then
5040 write (iulog,*) ' findsp not converging at point i, k ', i, k
5041 write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k),
5042 & enin(i)
5043 write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k),
5044 & enout(i)
5045 call endrun ('FINDSP')
5046 endif
5047 end do
5048 endif
5049 endif
5050 do i = 1,ncol
5051 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+
5052 &enout(i))) > 1.e-4_r8) then
5053 error_found = .true.
5054 endif
5055 end do
5056 if (error_found) then
5057 do i = 1,ncol
5058 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+
5059 &enout(i))) > 1.e-4_r8) then
5060 write (iulog,*) ' the enthalpy is not conserved for point ', i,
5061 & k, enin(i), enout(i)
5062 write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k),
5063 & enin(i)
5064 write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k),
5065 & enout(i)
5066 call endrun ('FINDSP')
5067 endif
5068 end do
5069 endif
5070
5071 end do
5072
5073 return
5074 end subroutine findsp1_water
5075#endif
5076!--jtb : end of GEOS5 exclusion (begins at top of findsp)
5077
5078!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
5084
5085 function derf(x)
5086 implicit real (a - h, o - z)
5087 real(r8) a,b,x
5088 dimension a(0 : 64), b(0 : 64)
5089 integer i,k
5090 data (a(i), i = 0, 12) / 0.00000000005958930743d0, -
5091 &0.00000000113739022964d0, 0.00000001466005199839d0, -
5092 &0.00000016350354461960d0, 0.00000164610044809620d0, -
5093 &0.00001492559551950604d0, 0.00012055331122299265d0, -
5094 &0.00085483269811296660d0, 0.00522397762482322257d0, -
5095 &0.02686617064507733420d0, 0.11283791670954881569d0, -
5096 &0.37612638903183748117d0, 1.12837916709551257377d0 /
5097 data (a(i), i = 13, 25) / 0.00000000002372510631d0, -
5098 &0.00000000045493253732d0, 0.00000000590362766598d0, -
5099 &0.00000006642090827576d0, 0.00000067595634268133d0, -
5100 &0.00000621188515924000d0, 0.00005103883009709690d0, -
5101 &0.00037015410692956173d0, 0.00233307631218880978d0, -
5102 &0.01254988477182192210d0, 0.05657061146827041994d0, -
5103 &0.21379664776456006580d0, 0.84270079294971486929d0 /
5104 data (a(i), i = 26, 38) / 0.00000000000949905026d0, -
5105 &0.00000000018310229805d0, 0.00000000239463074000d0, -
5106 &0.00000002721444369609d0, 0.00000028045522331686d0, -
5107 &0.00000261830022482897d0, 0.00002195455056768781d0, -
5108 &0.00016358986921372656d0, 0.00107052153564110318d0, -
5109 &0.00608284718113590151d0, 0.02986978465246258244d0, -
5110 &0.13055593046562267625d0, 0.67493323603965504676d0 /
5111 data (a(i), i = 39, 51) / 0.00000000000382722073d0, -
5112 &0.00000000007421598602d0, 0.00000000097930574080d0, -
5113 &0.00000001126008898854d0, 0.00000011775134830784d0, -
5114 &0.00000111992758382650d0, 0.00000962023443095201d0, -
5115 &0.00007404402135070773d0, 0.00050689993654144881d0, -
5116 &0.00307553051439272889d0, 0.01668977892553165586d0, -
5117 &0.08548534594781312114d0, 0.56909076642393639985d0 /
5118 data (a(i), i = 52, 64) / 0.00000000000155296588d0, -
5119 &0.00000000003032205868d0, 0.00000000040424830707d0, -
5120 &0.00000000471135111493d0, 0.00000005011915876293d0, -
5121 &0.00000048722516178974d0, 0.00000430683284629395d0, -
5122 &0.00003445026145385764d0, 0.00024879276133931664d0, -
5123 &0.00162940941748079288d0, 0.00988786373932350462d0, -
5124 &0.05962426839442303805d0, 0.49766113250947636708d0 /
5125 data (b(i), i = 0, 12) / -0.00000000029734388465d0,
5126 & 0.00000000269776334046d0, -0.00000000640788827665d0, -
5127 &0.00000001667820132100d0, -0.00000021854388148686d0,
5128 & 0.00000266246030457984d0, 0.00001612722157047886d0, -
5129 &0.00025616361025506629d0, 0.00015380842432375365d0,
5130 & 0.00815533022524927908d0, -0.01402283663896319337d0, -
5131 &0.19746892495383021487d0, 0.71511720328842845913d0 /
5132 data (b(i), i = 13, 25) / -0.00000000001951073787d0, -
5133 &0.00000000032302692214d0, 0.00000000522461866919d0,
5134 & 0.00000000342940918551d0, -0.00000035772874310272d0,
5135 & 0.00000019999935792654d0, 0.00002687044575042908d0, -
5136 &0.00011843240273775776d0, -0.00080991728956032271d0,
5137 & 0.00661062970502241174d0, 0.00909530922354827295d0, -
5138 &0.20160072778491013140d0, 0.51169696718727644908d0 /
5139 data (b(i), i = 26, 38) / 0.00000000003147682272d0, -
5140 &0.00000000048465972408d0, 0.00000000063675740242d0,
5141 & 0.00000003377623323271d0, -0.00000015451139637086d0, -
5142 &0.00000203340624738438d0, 0.00001947204525295057d0,
5143 & 0.00002854147231653228d0, -0.00101565063152200272d0,
5144 & 0.00271187003520095655d0, 0.02328095035422810727d0, -
5145 &0.16725021123116877197d0, 0.32490054966649436974d0 /
5146 data (b(i), i = 39, 51) / 0.00000000002319363370d0, -
5147 &0.00000000006303206648d0, -0.00000000264888267434d0,
5148 & 0.00000002050708040581d0, 0.00000011371857327578d0, -
5149 &0.00000211211337219663d0, 0.00000368797328322935d0,
5150 & 0.00009823686253424796d0, -0.00065860243990455368d0, -
5151 &0.00075285814895230877d0, 0.02585434424202960464d0, -
5152 &0.11637092784486193258d0, 0.18267336775296612024d0 /
5153 data (b(i), i = 52, 64) / -0.00000000000367789363d0,
5154 & 0.00000000020876046746d0, -0.00000000193319027226d0, -
5155 &0.00000000435953392472d0, 0.00000018006992266137d0, -
5156 &0.00000078441223763969d0, -0.00000675407647949153d0,
5157 & 0.00008428418334440096d0, -0.00017604388937031815d0, -
5158 &0.00239729611435071610d0, 0.02064129023876022970d0, -
5159 &0.06905562880005864105d0, 0.09084526782065478489d0 /
5160 w = abs(x)
5161 if (w .lt. 2.2d0) then
5162 t = w * w
5163 k = int(t)
5164 t = t - k
5165 k = k * 13
5166 y = ((((((((((((a(k) * t + a(k + 1)) * t + a(k + 2)) * t + a(k +
5167 & 3)) * t + a(k + 4)) * t + a(k + 5)) * t + a(k + 6)) * t + a(k +
5168 & 7)) * t + a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + a(k +
5169 & 11)) * t + a(k + 12)) * w
5170 else if (w .lt. 6.9d0) then
5171 k = int(w)
5172 t = w - k
5173 k = 13 * (k - 2)
5174 y = (((((((((((b(k) * t + b(k + 1)) * t + b(k + 2)) * t + b(k +
5175 & 3)) * t + b(k + 4)) * t + b(k + 5)) * t + b(k + 6)) * t + b(k +
5176 & 7)) * t + b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + b(k +
5177 & 11)) * t + b(k + 12)
5178 y = y * y
5179 y = y * y
5180 y = y * y
5181 y = 1 - y * y
5182 else
5183 y = 1
5184 end if
5185 if (x .lt. 0) y = -y
5186 derf = y
5187 end function derf
5188!
5189
5190
5191!cccccccccccccccccccccDONIFccccccccccccccccccccccccccccccccccccccccccccccccc
5192
5193
5194
5195!**********************************
5198 FUNCTION mui_hemp(T)
5199
5200
5201 real(r8) :: mui_hemp
5202 REAL(r8), intent(in) :: t
5203 REAL(r8) :: tc, mui, lambdai
5204 tc=t-273.15_r8
5205
5206 tc=min(max(tc, -70.0), -15.0)
5207
5208 if (tc > -27.0) then
5209 lambdai = 6.8_r8*exp(-0.096_r8*tc)
5210 else
5211 lambdai = 24.8_r8*exp(-0.049_r8*tc)
5212 end if
5213
5214 mui=(0.13_r8*(lambdai**0.64_r8))-two
5215 mui=max(mui, 1.5_r8)
5216 mui_hemp=mui
5217
5218
5219 END FUNCTION mui_hemp
5220
5221
5222!cccccccccccccccccccccDONIFccccccccccccccccccccccccccccccccccccccccccccccccc
5223
5224
5225
5226!**********************************
5229 FUNCTION mui_hemp_l(lambda)
5230
5231
5232 real(r8) :: mui_hemp_l
5233 REAL(r8), intent(in) :: lambda
5234 REAL(r8) :: tc, mui, lx
5235 lx = lambda*0.01
5236
5237 mui=(0.008_r8*(lx**0.87_r8))
5238 mui=max(min(mui, 10.0_r8), 0.1_r8)
5239 mui_hemp_l=mui !Anning for multithread to work
5240
5241
5242 END FUNCTION mui_hemp_l
5243
5244
5245
5248 FUNCTION gamma_incomp(muice, x)
5249
5250
5251
5252 real(r8) :: gamma_incomp
5253 REAL(r8), intent(in) :: muice, x
5254 REAL(r8) :: xog, kg, alfa, auxx
5255 alfa = min(max(muice+1._r8, 1._r8), 20._r8)
5256
5257 xog = log(alfa -0.3068_r8)
5258 kg = 1.44818*(alfa**0.5357_r8)
5259 auxx = max(min(kg*(log(x)-xog), 30._r8), -30._r8)
5260 gamma_incomp= one/(one +exp(-auxx))
5261 gamma_incomp = max(gamma_incomp, 1.0e-20)
5262
5263 END FUNCTION gamma_incomp
5264
5265
5266!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
5268 FUNCTION gamma(X)
5269
5270!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
5271
5272!D DOUBLE PRECISION FUNCTION DGAMMA(X)
5273!----------------------------------------------------------------------
5274!
5275! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
5276! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
5277! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
5278! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
5279! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
5280! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
5281! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
5282! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
5283! MACHINE-DEPENDENT CONSTANTS.
5284!
5285!
5286!*******************************************************************
5287!*******************************************************************
5288!
5289! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
5290!
5291! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
5292! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
5293! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
5294! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
5295! GAMMA(XBIG) = BETA**MAXEXP
5296! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
5297! APPROXIMATELY BETA**MAXEXP
5298! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
5299! 1.0+EPS .GT. 1.0
5300! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
5301! 1/XMININ IS MACHINE REPRESENTABLE
5302!
5303! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
5304!
5305! BETA MAXEXP XBIG
5306!
5307! CRAY-1 (S.P.) 2 8191 966.961
5308! CYBER 180/855
5309! UNDER NOS (S.P.) 2 1070 177.803
5310! IEEE (IBM/XT,
5311! SUN, ETC.) (S.P.) 2 128 35.040
5312! IEEE (IBM/XT,
5313! SUN, ETC.) (D.P.) 2 1024 171.624
5314! IBM 3033 (D.P.) 16 63 57.574
5315! VAX D-FORMAT (D.P.) 2 127 34.844
5316! VAX G-FORMAT (D.P.) 2 1023 171.489
5317!
5318! XINF EPS XMININ
5319!
5320! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
5321! CYBER 180/855
5322! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
5323! IEEE (IBM/XT,
5324! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
5325! IEEE (IBM/XT,
5326! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
5327! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
5328! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
5329! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
5330!
5331!*******************************************************************
5332!*******************************************************************
5333!
5334! ERROR RETURNS
5335!
5336! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
5337! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
5338! TO BE FREE OF UNDERFLOW AND OVERFLOW.
5339!
5340!
5341! INTRINSIC FUNCTIONS REQUIRED ARE:
5342!
5343! INT, DBLE, EXP, LOG, REAL, SIN
5344!
5345!
5346! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
5347! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
5348! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
5349! (ED.), SPRINGER VERLAG, BERLIN, 1976.
5350!
5351! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
5352! SONS, NEW YORK, 1968.
5353!
5354! LATEST MODIFICATION: OCTOBER 12, 1989
5355!
5356! AUTHORS: W. J. CODY AND L. STOLTZ
5357! APPLIED MATHEMATICS DIVISION
5358! ARGONNE NATIONAL LABORATORY
5359! ARGONNE, IL 60439
5360!
5361!----------------------------------------------------------------------
5362 INTEGER i,n
5363 LOGICAL parity
5364
5365 REAL(r8) :: gamma,conv,fact,res,sum,x,xden,xnum,y,y1,ysq,z
5366 real(r8) :: c(7),p(8),q(8)
5367!----------------------------------------------------------------------
5368! MATHEMATICAL CONSTANTS
5369!----------------------------------------------------------------------
5370 real(r8), parameter :: one=1.0e0_r8, half=0.5e0_r8, &
5371 & twelve=12.0e0_r8, two=2.0e0_r8, &
5372 & zero=0.0e0_r8, &
5373 & pi=3.1415926535897932384626434e0_r8, &
5374 & sqrtpi=0.9189385332046727417803297e0_r8
5375
5376!D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
5377!D 1 SQRTPI/0.9189385332046727417803297D0/,
5378!D 2 PI/3.1415926535897932384626434D0/
5379!----------------------------------------------------------------------
5380! MACHINE DEPENDENT PARAMETERS
5381!----------------------------------------------------------------------
5382 real(r8), parameter :: xbig=35.040e0_r8, xminin=1.18e-38_r8, &
5383 & eps=1.19e-7_r8, xinf=3.4e38_r8
5384
5385!D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
5386!D 1 XINF/1.79D308/
5387!----------------------------------------------------------------------
5388! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
5389! APPROXIMATION OVER (1,2).
5390!----------------------------------------------------------------------
5391 DATA p/-1.71618513886549492533811e+0_r8,
5392 &2.47656508055759199108314e+1_r8, -3.79804256470945635097577e+2_r8,
5393 &6.29331155312818442661052e+2_r8, 8.66966202790413211295064e+2_r8,-
5394 &3.14512729688483675254357e+4_r8, -3.61444134186911729807069e+4_r8,
5395 &6.64561438202405440627855e+4_r8/
5396 DATA q/-3.08402300119738975254353e+1_r8,
5397 &3.15350626979604161529144e+2_r8, -1.01515636749021914166146e+3_r8,
5398 &-3.10777167157231109440444e+3_r8, 2.25381184209801510330112e+4_r8,
5399 &4.75584627752788110767815e+3_r8, -1.34659959864969306392456e+5_r8,
5400 &-1.15132259675553483497211e+5_r8/
5401!D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
5402!D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
5403!D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
5404!D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
5405!D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
5406!D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
5407!D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3,
5408!D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
5409!----------------------------------------------------------------------
5410! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
5411!----------------------------------------------------------------------
5412 DATA c/-1.910444077728e-03_r8,8.4171387781295e-04_r8, -
5413 &5.952379913043012e-04_r8,7.93650793500350248e-04_r8, -
5414 &2.777777777777681622553e-03_r8,8.333333333333333331554247e-02_r8,
5415 & 5.7083835261e-03_r8/
5416!D DATA C/-1.910444077728D-03,8.4171387781295D-04,
5417!D 1 -5.952379913043012D-04,7.93650793500350248D-04,
5418!D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
5419!D 3 5.7083835261D-03/
5420!----------------------------------------------------------------------
5421! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
5422!----------------------------------------------------------------------
5423 conv(i) = real(i,r8)
5424!D CONV(I) = DBLE(I)
5425 parity = .false.
5426 fact = one
5427 n = 0
5428 y = x
5429 IF(y <= zero) THEN
5430!----------------------------------------------------------------------
5431! ARGUMENT IS NEGATIVE
5432!----------------------------------------------------------------------
5433 y = -x
5434 y1 = aint(y)
5435 res = y - y1
5436 IF(res.NE.zero)THEN
5437 IF(y1.NE.aint(y1*half)*two) parity = .true.
5438 fact = -pi/sin(pi*res)
5439 y = y + one
5440 ELSE
5441 res=xinf
5442 GOTO 900
5443 ENDIF
5444 ENDIF
5445!----------------------------------------------------------------------
5446! ARGUMENT IS POSITIVE
5447!----------------------------------------------------------------------
5448 IF(y.LT.eps)THEN
5449!----------------------------------------------------------------------
5450! ARGUMENT .LT. EPS
5451!----------------------------------------------------------------------
5452 IF(y.GE.xminin)THEN
5453 res = one/y
5454 ELSE
5455 res = xinf
5456 GOTO 900
5457 ENDIF
5458 ELSEIF(y.LT.twelve)THEN
5459 y1 = y
5460 IF(y.LT.one)THEN
5461!----------------------------------------------------------------------
5462! 0.0 .LT. ARGUMENT .LT. 1.0
5463!----------------------------------------------------------------------
5464 z = y
5465 y = y + one
5466 ELSE
5467!----------------------------------------------------------------------
5468! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
5469!----------------------------------------------------------------------
5470 n = int(y) - 1
5471 y = y - conv(n)
5472 z = y - one
5473 ENDIF
5474!----------------------------------------------------------------------
5475! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
5476!----------------------------------------------------------------------
5477 xnum = zero
5478 xden = one
5479 DO 260 i=1,8
5480 xnum = (xnum+p(i))*z
5481 xden = xden*z + q(i)
5482260 CONTINUE
5483 res = xnum/xden + one
5484 IF(y1.LT.y)THEN
5485!----------------------------------------------------------------------
5486! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
5487!----------------------------------------------------------------------
5488 res = res/y1
5489 ELSEIF(y1.GT.y)THEN
5490!----------------------------------------------------------------------
5491! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
5492!----------------------------------------------------------------------
5493 DO 290 i=1,n
5494 res = res*y
5495 y = y + one
5496290 CONTINUE
5497 ENDIF
5498 ELSE
5499!----------------------------------------------------------------------
5500! EVALUATE FOR ARGUMENT .GE. 12.0,
5501!----------------------------------------------------------------------
5502 IF(y.LE.xbig)THEN
5503 ysq = y*y
5504 sum = c(7)
5505 DO 350 i=1,6
5506 sum = sum / ysq + c(i)
5507350 CONTINUE
5508 sum = sum / y - y + sqrtpi
5509 sum = sum + (y-half)*log(y)
5510 res = exp(sum)
5511 ELSE
5512 res = xinf
5513 GOTO 900
5514 ENDIF
5515 ENDIF
5516!----------------------------------------------------------------------
5517! FINAL ADJUSTMENTS AND RETURN
5518!----------------------------------------------------------------------
5519 IF(parity) res = -res
5520 IF(fact.NE.one) res = fact/res
5521900 gamma = res
5522!D900 DGAMMA = RES
5523 RETURN
5524! ---------- LAST LINE OF GAMMA ----------
5525 END function gamma
5526
5527
5528
5529 end module cldwat2m_micro
subroutine findsp1_water(lchnk, ncol, q, t, p, tsp, qsp)
This subroutine finds the wet bulb temperature for a given t and q in a longitude height section.
real(r8) function mui_hemp_l(lambda)
This function.
real(r8), private ci
real(r8) function mui_hemp(t)
This subroutine.
real(r8), private bs
real(r8), private rhosn
real(r8), parameter rh2o
real(r8), parameter three
real(r8), parameter tmelt
real(r8), private minrefl
real(r8), private di
subroutine, public mmicro_pcond(lchnk, ncol, deltatin, tn, ttend, pcols, pver, qn, qtend, cwtend, qc, qi, nc, ni, fprcp, qrn, qsnw, nrn, nsnw, p, pdel, cldn, liqcldf, icecldf, cldo, pint, rpdel, zm, rate1ord_cw2pr_st, naai, npccnin, rndst, nacon, rhdfda, rhu00, fice, tlat, qvlat, qctend, qitend, nctend, nitend, effc, effc_fn, effi, prect, preci, nevapr, evapsnow, prain, prodsnow, cmeout, deffi, pgamrad, lamcrad, qsout2, dsout2, qrout2, drout2, qcsevap, qisevap, qvres, cmeiout, vtrmc, vtrmi, qcsedten, qisedten, prao, prco, mnuccco, mnuccto, msacwio, psacwso, bergso, bergo, melto, homoo, qcreso, prcio, praio, qireso, mnuccro, pracso, meltsdt, frzrdt, ncal, ncai, mnuccdo, nnuccto, nsout2, nrout2, ncnst, ninst, nimm, miu_disp, nsoot, rnsoot, ui_scale, dcrit, nnuccdo, nnuccco, nsacwio, nsubio, nprcio, npraio, npccno, npsacwso, nsubco, nprao, nprc1o, tlataux, nbincontactdust, lprint, xlat, xlon, rhc)
This subroutine is the microphysics routine for each timestep goes here...
real(r8), private ac
real(r8), private xlf
real(r8), private oneodi
real(r8), private mindbz
real(r8), parameter oneb3
real(r8), parameter cpair
real(r8), parameter latice
real(r8), parameter gravit
real(r8), private cr
real(r8), private ai
real(r8), private ds
real(r8), parameter half
real(r8), parameter one
real(r8), private cs
real(r8), parameter epsilon
integer, parameter iulog
real(r8), private cons3
real(r8), private ar
real(r8), private lammaxi
real(r8), private pirhosn
real(r8), private csmin
real(r8), private f1r
real(r8), private qsmall
logical, public liu_in
real(r8), private ginv
real(r8), private eii
function, public derf(x)
error function in single precision
real(r8), private lammini
real(r8), parameter mwh2o
real(r8), private rhosu
real(r8), private lamminr
real(r8), private lammaxr
real(r8), private cpp
real(r8), parameter latvap
real(r8), private f2s
real(r8), private rv
real(r8), private pirhoi
real(r8), private cons1
real(r8), private f2r
real(r8), private xlfocp
real(r8), private rhoi
real(r8), private cpoxlf
real(r8), private ecc
real(r8) function gamma_incomp(muice, x)
This function.
real(r8), parameter rair
real(r8), private mi0
real(r8), private pirhow
real(r8), private qvsmall
real(r8), private xxls
real(r8), private r
real(r8), private tt0
real(r8), parameter four
real(r8), private ts_auto_ice
real(r8), private bimm
real(r8), parameter rhmaxi
real(r8), private aimm
subroutine, public ini_micro(dcs_, qcvar_, ts_auto_ice_)
This subroutine initializes constants for MG microphysics.
real(r8), private qcvar
real(r8), private br
subroutine findsp1(lchnk, ncol, q, t, p, tsp, qsp)
This subroutine finds the wet bulb temperature for a given t and q in a longitude hight section.
real(r8), parameter rhoh2o
real(r8), private bi
real(r8), private dr
real(r8), private dcs
real(r8), private g
real(r8), private bc
real(r8), private as
real(r8), parameter, private tmax_fsnow
parameters for snow/rain fraction for convective clouds
real(r8), private csmax
real(r8), parameter onebcp
real(r8), private rhow
real(r8), private f1s
real(r8), private xxlv
real(r8), private cons2
real(r8), parameter five
real(r8), private ecr
real(r8), private cons4
real(r8), parameter two
real(r8) function, public gamma(x)
real(r8), parameter r_universal
real(r8), private rin
real(r8), parameter rhmini
real(r8), private cons5
subroutine, public vqsatd(t, p, es, qs, gam, len)
This subroutine is the utility procedure to look up and return saturation vapor pressure from precomp...
subroutine, public vqsatd_water(t, p, es, qs, gam, len)
subroutine, public vqsatd2(t, p, es, qs, dqsdt, len)
real(kp) function, public polysvp(t, typ)
subroutine, public gestbl(tmn, tmx, trice, ip, epsil, latvap, latice, rh2o, cpair, tmeltx)
This subroutine builds saturation vapor pressure table for later procedure.
subroutine, public vqsatd2_single(t, p, es, qs, dqsdt)
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36
Definition sflx.f:3