FV3 Bundle
convection_tl.F90
Go to the documentation of this file.
1 module convection_tl
2 
3 USE convection
4 
5 IMPLICIT NONE
6 
7 PRIVATE
8 PUBLIC :: rase_d, rase0_d, rase_tracer_d
9 
10 CONTAINS
11 
12 ! Generated by TAPENADE (INRIA, Tropics team)
13 ! Tapenade 3.9 (r5096) - 24 Feb 2014 16:53
14 !
15 ! Differentiation of rase in forward (tangent) mode:
16 ! variations of useful results: clw cnv_prc3 tho qho vho cnv_updfrc
17 ! uho flxd
18 ! with respect to varying inputs: tho qho vho uho
19 ! RW status of diff variables: clw:out cnv_prc3:out tho:in-out
20 ! qho:in-out vho:in-out cnv_updfrc:out uho:in-out
21 ! flxd:out
22 SUBROUTINE rase_d(idim, irun, k0, icmin, dt, cons_cp, cons_alhl, &
23 & cons_grav, cons_rgas, cons_h2omw, cons_airmw, cons_vireps, seedras, &
24 & sige, kcbl, wgt0, wgt1, frland, ts, tho, thod, qho, qhod, uho, uhod, &
25 & vho, vhod, co_auto, ple, clw, clwd, flxd, flxdd, cnv_prc3, cnv_prc3d, &
26 & cnv_updfrc, cnv_updfrcd, rasparams, estblx)
27  IMPLICIT NONE
28 !INPUTS
29  INTEGER, INTENT(IN) :: idim, irun, k0, icmin
30  REAL*8, DIMENSION(idim, k0 + 1), INTENT(IN) :: ple
31  REAL*8, DIMENSION(k0 + 1), INTENT(IN) :: sige
32  REAL*8, INTENT(IN) :: dt, cons_cp, cons_alhl, cons_grav, cons_rgas
33  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_vireps
34  INTEGER, DIMENSION(idim), INTENT(IN) :: seedras
35  INTEGER, DIMENSION(idim), INTENT(IN) :: kcbl
36  REAL*8, DIMENSION(idim), INTENT(IN) :: ts, frland
37  REAL*8, DIMENSION(idim), INTENT(IN) :: co_auto
38  REAL*8, DIMENSION(idim, k0), INTENT(IN) :: wgt0, wgt1
39  REAL*8, DIMENSION(:), INTENT(IN) :: rasparams
40  REAL*8, DIMENSION(:), INTENT(IN) :: estblx
41 !OUTPUTS
42  REAL*8, DIMENSION(idim, k0), INTENT(OUT) :: clw, flxd
43  REAL*8, DIMENSION(idim, k0), INTENT(OUT) :: clwd, flxdd
44  REAL*8, DIMENSION(idim, k0), INTENT(OUT) :: cnv_prc3
45  REAL*8, DIMENSION(idim, k0), INTENT(OUT) :: cnv_prc3d
46  REAL*8, DIMENSION(idim, k0), INTENT(OUT) :: cnv_updfrc
47  REAL*8, DIMENSION(idim, k0), INTENT(OUT) :: cnv_updfrcd
48 !PROGNOSTIC
49  REAL*8, DIMENSION(idim, k0), INTENT(INOUT) :: tho, qho, uho, vho
50  REAL*8, DIMENSION(idim, k0), INTENT(INOUT) :: thod, qhod, uhod, vhod
51 !LOCALS
52  INTEGER :: i, ic, l, kk, k
53 !Parameters
54  REAL*8, PARAMETER :: onepkap=1.+2./7., daylen=86400.0
55  REAL*8, PARAMETER :: rhmax=0.9999
56  REAL*8, PARAMETER :: cbl_qpert=0.0, cbl_tpert=1.0
57  REAL*8, PARAMETER :: cbl_tpert_mxocn=2.0, cbl_tpert_mxlnd=4.0
58 !Constants
59  REAL*8 :: grav, cp, alhl, cpbg, alhi, cpi, gravi, ddt, lbcp
60 !Rasparams
61  REAL*8 :: fricfac, cli_crit, rasal1, rasal2
62  REAL*8 :: friclambda
63  REAL*8 :: sdqv2, sdqv3, sdqvt1
64  REAL*8 :: acritfac, pblfrac, autorampb
65  REAL*8 :: maxdallowed, rhmn, rhmx
66  REAL*8 :: mxdiam
67  REAL*8 :: tx2, tx3, akm, acr, alm, tth, qqh, dqx
68  REAL*8 :: tx2d, tx3d, akmd, almd
69  REAL*8 :: wfn, tem, trg, trgexp, evp, wlq, qcc
70  REAL*8 :: wfnd, temd, trgd, wlqd, qccd
71  REAL*8 :: cli, te_a, c00_x, cli_crit_x, toki
72  REAL*8 :: clid, te_ad, c00_xd, cli_crit_xd, tokid
73  REAL*8 :: dt_lyr, rate, cvw_x, closs, f2, f3, f4
74  REAL*8 :: dt_lyrd, rated, cvw_xd, clossd, f2d
75  REAL*8 :: wght0, prcbl, rndu
76  REAL*8 :: lambda_min, lambda_max
77  REAL*8 :: tpert, qpert
78  REAL*8 :: tpertd
79  REAL*8 :: uht, vht
80  REAL*8 :: uhtd, vhtd
81  REAL*8, DIMENSION(k0) :: poi_sv, qoi_sv, uoi_sv, voi_sv
82  REAL*8, DIMENSION(k0) :: poi_svd, qoi_svd, uoi_svd, voi_svd
83  REAL*8, DIMENSION(k0) :: poi, qoi, uoi, voi, dqq, bet, gam, cll
84  REAL*8, DIMENSION(k0) :: poid, qoid, uoid, void0, dqqd, betd, gamd, &
85 & clld
86  REAL*8, DIMENSION(k0) :: poi_c, qoi_c
87  REAL*8, DIMENSION(k0) :: poi_cd, qoi_cd
88  REAL*8, DIMENSION(k0) :: prh, pri, ght, dpt, dpb, pki
89  REAL*8, DIMENSION(k0) :: prhd, prid, ghtd, dptd, dpbd, pkid
90  REAL*8, DIMENSION(k0) :: ucu, vcu
91  REAL*8, DIMENSION(k0) :: ucud, vcud
92  REAL*8, DIMENSION(k0) :: cln, rns, pol
93  REAL*8, DIMENSION(k0) :: rnsd, pold
94  REAL*8, DIMENSION(k0) :: qst, ssl, rmf, rnn, rn1, rmfc, rmfp
95  REAL*8, DIMENSION(k0) :: qstd, ssld, rnnd, rmfpd
96  REAL*8, DIMENSION(k0) :: gms, eta, gmh, eht, gm1, hcc, rmfd
97  REAL*8, DIMENSION(k0) :: gmsd, etad, gmhd, ehtd, gm1d, hccd, rmfdd
98  REAL*8, DIMENSION(k0) :: hol, hst, qol, zol, hcld, cll0, cllx, clli
99  REAL*8, DIMENSION(k0) :: hold, hstd, qold, zold, hcldd, cll0d
100  REAL*8, DIMENSION(k0) :: bke, cvw, updfrc
101  REAL*8, DIMENSION(k0) :: cvwd, updfrcd
102  REAL*8, DIMENSION(k0) :: rasal, updfrp, bk2, dll0, dllx
103  REAL*8, DIMENSION(k0) :: rasald, updfrpd, bk2d
104  REAL*8, DIMENSION(k0) :: wght, massf
105  REAL*8, DIMENSION(k0) :: wghtd
106  REAL*8, DIMENSION(k0) :: qss, dqs, pf, pk, tempf, zlo
107  REAL*8, DIMENSION(k0) :: qssd, dqsd, tempfd, zlod
108  REAL*8, DIMENSION(k0 + 1) :: prj, prs, qht, sht, zet, zle, pke
109  REAL*8, DIMENSION(k0+1) :: prjd, prsd, qhtd, shtd, zetd, zled
110  INTRINSIC max
111  INTRINSIC min
112  INTRINSIC sqrt
113  INTRINSIC exp
114  INTRINSIC sum
115  REAL*8, DIMENSION(k0+1) :: pwx1
116  REAL*8 :: pwy1
117  REAL*8, DIMENSION(k0) :: pwx10
118  REAL*8 :: pwx11
119  REAL*8 :: pwr1
120  REAL*8 :: arg1
121  REAL*8 :: arg1d
122  REAL*8 :: max2d
123  REAL*8 :: x1
124  REAL*8 :: max1d
125  REAL*8 :: x1d
126  REAL*8 :: max2
127  REAL*8 :: max1
128  REAL*8 :: y1
129 !Initialize Local Arrays
130  poi = 0.0
131  qoi = 0.0
132  uoi = 0.0
133  voi = 0.0
134  dqq = 0.0
135  bet = 0.0
136  gam = 0.0
137  poi_c = 0.0
138  qoi_c = 0.0
139  prh = 0.0
140  pri = 0.0
141  ght = 0.0
142  dpt = 0.0
143  dpb = 0.0
144  pki = 0.0
145  ucu = 0.0
146  vcu = 0.0
147  cln = 0.0
148  pol = 0.0
149  qst = 0.0
150  ssl = 0.0
151  rmf = 0.0
152  rnn = 0.0
153  rn1 = 0.0
154  gms = 0.0
155  eta = 0.0
156  gmh = 0.0
157  eht = 0.0
158  gm1 = 0.0
159  hcc = 0.0
160  hol = 0.0
161  hst = 0.0
162  qol = 0.0
163  zol = 0.0
164  hcld = 0.0
165  bke = 0.0
166  cvw = 0.0
167  updfrc = 0.0
168  rasal = 0.0
169  updfrp = 0.0
170  bk2 = 0.0
171  wght = 0.0
172  massf = 0.0
173  qss = 0.0
174  dqs = 0.0
175  pf = 0.0
176  pk = 0.0
177  tempf = 0.0
178  zlo = 0.0
179  prj = 0.0
180  prs = 0.0
181  qht = 0.0
182  sht = 0.0
183  zet = 0.0
184  zle = 0.0
185  pke = 0.0
186 !Initialize Outputs
187  cnv_prc3 = 0.
188  cnv_updfrc = 0.
189  flxd = 0.
190  clw = 0.
191 ! --- 1
192  fricfac = rasparams(1)
193 ! --- 4
194  cli_crit = rasparams(4)
195 ! --- 5
196  rasal1 = rasparams(5)
197 ! --- 6
198  rasal2 = rasparams(6)
199 ! --- 11
200  friclambda = rasparams(11)
201 ! --- 14
202  sdqv2 = rasparams(14)
203 ! --- 15
204  sdqv3 = rasparams(15)
205 ! --- 16
206  sdqvt1 = rasparams(16)
207 ! --- 17
208  acritfac = rasparams(17)
209 ! --- 20
210  pblfrac = rasparams(20)
211 ! --- 21
212  autorampb = rasparams(21)
213 ! --- 24
214  rhmn = rasparams(24)
215 ! --- 24
216  maxdallowed = rasparams(23)
217 ! --- 25
218  rhmx = rasparams(25)
219  grav = cons_grav
220  alhl = cons_alhl
221  cp = cons_cp
222  cpi = 1.0/cp
223  alhi = 1.0/alhl
224  gravi = 1.0/grav
225  cpbg = cp*gravi
226  ddt = daylen/dt
227  lbcp = alhl*cpi
228  i = 1
229 !CALL FINDBASE
230  k = kcbl(i)
231  IF (k .GT. 0) THEN
232 !Get saturation specific humidity and gradient wrt to T
233  pwx1 = ple(i, :)/1000.
234  pwy1 = cons_rgas/cons_cp
235  pke = pwx1**pwy1
236  pf = 0.5*(ple(i, 1:k0)+ple(i, 2:k0+1))
237  pwx10 = pf/1000.
238  pwy1 = cons_rgas/cons_cp
239  pk = pwx10**pwy1
240  tempfd = pk*thod(i, :)
241  tempf = tho(i, :)*pk
242  zle = 0.0
243  zlo = 0.0
244  zled(k0+1) = 0.0_8
245  zle(k0+1) = 0.
246  zlod = 0.0_8
247  zled = 0.0_8
248  DO l=k0,1,-1
249  zled(l) = thod(i, l)*(1.+cons_vireps*qho(i, l)) + tho(i, l)*&
250 & cons_vireps*qhod(i, l)
251  zle(l) = tho(i, l)*(1.+cons_vireps*qho(i, l))
252  zlod(l) = zled(l+1) + cons_cp*(pke(l+1)-pk(l))*zled(l)/cons_grav
253  zlo(l) = zle(l+1) + cons_cp/cons_grav*(pke(l+1)-pk(l))*zle(l)
254  zled(l) = zlod(l) + cons_cp*(pk(l)-pke(l))*zled(l)/cons_grav
255  zle(l) = zlo(l) + cons_cp/cons_grav*(pk(l)-pke(l))*zle(l)
256  END DO
257  tpertd = cbl_tpert*(-tempfd(k0)-cons_grav*zlod(k0)/cons_cp)
258  tpert = cbl_tpert*(ts(i)-(tempf(k0)+cons_grav*zlo(k0)/cons_cp))
259 !* ( QSSFC - Q(:,:,K0) ) [CBL_QPERT = 0.0]
260  qpert = cbl_qpert
261  IF (tpert .LT. 0.0) THEN
262  tpert = 0.0
263  tpertd = 0.0_8
264  ELSE
265  tpert = tpert
266  END IF
267  IF (qpert .LT. 0.0) THEN
268  qpert = 0.0
269  ELSE
270  qpert = qpert
271  END IF
272  IF (frland(i) .LT. 0.1) THEN
273  IF (tpert .GT. cbl_tpert_mxocn) THEN
274  tpert = cbl_tpert_mxocn
275  tpertd = 0.0_8
276  ELSE
277  tpert = tpert
278  END IF
279  ELSE IF (tpert .GT. cbl_tpert_mxlnd) THEN
280  tpert = cbl_tpert_mxlnd
281  tpertd = 0.0_8
282  ELSE
283  tpert = tpert
284  END IF
285  CALL dqsat_ras_d(dqs, dqsd, qss, qssd, tempf, tempfd, pf, k0, estblx&
286 & , cons_h2omw, cons_airmw)
287  DO kk=icmin,k+1
288  prjd(kk) = 0.0_8
289  prj(kk) = pke(kk)
290  END DO
291  prsd(icmin:k0+1) = 0.0_8
292  prs(icmin:k0+1) = ple(i, icmin:k0+1)
293  poid = 0.0_8
294  poid(icmin:k) = thod(i, icmin:k)
295  poi(icmin:k) = tho(i, icmin:k)
296  qoid = 0.0_8
297  qoid(icmin:k) = qhod(i, icmin:k)
298  qoi(icmin:k) = qho(i, icmin:k)
299  uoid = 0.0_8
300  uoid(icmin:k) = uhod(i, icmin:k)
301  uoi(icmin:k) = uho(i, icmin:k)
302  void0 = 0.0_8
303  void0(icmin:k) = vhod(i, icmin:k)
304  voi(icmin:k) = vho(i, icmin:k)
305  qstd = 0.0_8
306  qstd(icmin:k) = qssd(icmin:k)
307  qst(icmin:k) = qss(icmin:k)
308  dqqd = 0.0_8
309  dqqd(icmin:k) = dqsd(icmin:k)
310  dqq(icmin:k) = dqs(icmin:k)
311 !Mass fraction of each layer below cloud base
312  massf(:) = wgt0(i, :)
313 !RESET PRESSURE at bottom edge of CBL
314  prcbl = prs(k)
315  DO l=k,k0
316  prcbl = prcbl + massf(l)*(prs(l+1)-prs(l))
317  END DO
318  prsd(k+1) = 0.0_8
319  prs(k+1) = prcbl
320  pwx11 = prs(k+1)/1000.
321  pwy1 = cons_rgas/cons_cp
322  prjd(k+1) = 0.0_8
323  prj(k+1) = pwx11**pwy1
324  DO l=k,icmin,-1
325  pold(l) = 0.0_8
326  pol(l) = 0.5*(prs(l)+prs(l+1))
327  prhd(l) = 0.0_8
328  prh(l) = (prs(l+1)*prj(l+1)-prs(l)*prj(l))/(onepkap*(prs(l+1)-prs(&
329 & l)))
330  pkid(l) = 0.0_8
331  pki(l) = 1.0/prh(l)
332  dptd(l) = 0.0_8
333  dpt(l) = prh(l) - prj(l)
334  dpbd(l) = 0.0_8
335  dpb(l) = prj(l+1) - prh(l)
336  prid(l) = 0.0_8
337  pri(l) = .01/(prs(l+1)-prs(l))
338  END DO
339 !RECALCULATE PROFILE QUAN. IN LOWEST STRAPPED LAYER
340  IF (k .LE. k0) THEN
341  poid(k) = 0.0_8
342  poi(k) = 0.
343  qoid(k) = 0.0_8
344  qoi(k) = 0.
345  uoid(k) = 0.0_8
346  uoi(k) = 0.
347  void0(k) = 0.0_8
348  voi(k) = 0.
349 !SPECIFY WEIGHTS GIVEN TO EACH LAYER WITHIN SUBCLOUD "SUPERLAYER"
350  wght = 0.
351  DO l=k,k0
352  wghtd(l) = 0.0_8
353  wght(l) = massf(l)*(ple(i, l+1)-ple(i, l))/(prs(k+1)-prs(k))
354  END DO
355  DO l=k,k0
356  poid(k) = poid(k) + wght(l)*thod(i, l)
357  poi(k) = poi(k) + wght(l)*tho(i, l)
358  qoid(k) = qoid(k) + wght(l)*qhod(i, l)
359  qoi(k) = qoi(k) + wght(l)*qho(i, l)
360  uoid(k) = uoid(k) + wght(l)*uhod(i, l)
361  uoi(k) = uoi(k) + wght(l)*uho(i, l)
362  void0(k) = void0(k) + wght(l)*vhod(i, l)
363  voi(k) = voi(k) + wght(l)*vho(i, l)
364  END DO
365  CALL dqsats_ras_d(dqq(k), dqqd(k), qst(k), qstd(k), poi(k)*prh(k)&
366 & , prh(k)*poid(k), pol(k), estblx, cons_h2omw, &
367 & cons_airmw)
368  END IF
369  IF (seedras(i)/1000000. .LT. 1e-6) THEN
370  rndu = 1e-6
371  ELSE
372  rndu = seedras(i)/1000000.
373  END IF
374  pwr1 = rndu**(-(1./2.))
375  mxdiam = maxdallowed*pwr1
376  gm1d = 0.0_8
377  ghtd = 0.0_8
378  betd = 0.0_8
379  gamd = 0.0_8
380  DO l=k,icmin,-1
381 !*
382  betd(l) = pki(l)*dqqd(l)
383  bet(l) = dqq(l)*pki(l)
384 !*
385  gamd(l) = -(pki(l)*lbcp*dqqd(l)/(1.0+lbcp*dqq(l))**2)
386  gam(l) = pki(l)/(1.0+lbcp*dqq(l))
387  IF (l .LT. k) THEN
388  ghtd(l+1) = dpb(l)*gamd(l) + dpt(l+1)*gamd(l+1)
389  ght(l+1) = gam(l)*dpb(l) + gam(l+1)*dpt(l+1)
390  gm1d(l+1) = 0.5*lbcp*((dqqd(l)*alhl*(1.0+lbcp*dqq(l))-dqq(l)*&
391 & alhl*lbcp*dqqd(l))/(alhl*(1.0+lbcp*dqq(l)))**2+(dqqd(l+1)*alhl&
392 & *(1.0+lbcp*dqq(l+1))-dqq(l+1)*alhl*lbcp*dqqd(l+1))/(alhl*(1.0+&
393 & lbcp*dqq(l+1)))**2)
394  gm1(l+1) = 0.5*lbcp*(dqq(l)/(alhl*(1.0+lbcp*dqq(l)))+dqq(l+1)/(&
395 & alhl*(1.0+lbcp*dqq(l+1))))
396  END IF
397  END DO
398  rns = 0.
399  cll = 0.
400  rmf = 0.
401  rmfd = 0.
402  rmfc = 0.
403  rmfp = 0.
404  cll0 = 0.
405  dll0 = 0.
406  cllx = 0.
407  dllx = 0.
408  clli = 0.
409  poi_svd = poid
410  poi_sv = poi
411  qoi_svd = qoid
412  qoi_sv = qoi
413  uoi_svd = uoid
414  uoi_sv = uoi
415  voi_svd = void0
416  voi_sv = voi
417  cvw = 0.0
418  updfrc = 0.0
419  updfrp = 0.0
420 ! HOL initialized here in order not to confuse Valgrind debugger
421  hol = 0.
422  zetd(k+1) = 0.0_8
423  zet(k+1) = 0
424  shtd = 0.0_8
425  shtd(k+1) = cp*prj(k+1)*poid(k)
426  sht(k+1) = cp*poi(k)*prj(k+1)
427  hstd = 0.0_8
428  qold = 0.0_8
429  ssld = 0.0_8
430  zetd = 0.0_8
431  zold = 0.0_8
432  hold = 0.0_8
433  DO l=k,icmin,-1
434  IF (qst(l)*rhmax .GT. qoi(l)) THEN
435  qold(l) = qoid(l)
436  qol(l) = qoi(l)
437  ELSE
438  qold(l) = rhmax*qstd(l)
439  qol(l) = qst(l)*rhmax
440  END IF
441  IF (0.000 .LT. qol(l)) THEN
442  qol(l) = qol(l)
443  ELSE
444  qold(l) = 0.0_8
445  qol(l) = 0.000
446  END IF
447  ssld(l) = cp*prj(l+1)*poid(l) + grav*zetd(l+1)
448  ssl(l) = cp*prj(l+1)*poi(l) + grav*zet(l+1)
449  hold(l) = ssld(l) + alhl*qold(l)
450  hol(l) = ssl(l) + qol(l)*alhl
451  hstd(l) = ssld(l) + alhl*qstd(l)
452  hst(l) = ssl(l) + qst(l)*alhl
453  temd = (prj(l+1)-prj(l))*cpbg*poid(l)
454  tem = poi(l)*(prj(l+1)-prj(l))*cpbg
455  zetd(l) = zetd(l+1) + temd
456  zet(l) = zet(l+1) + tem
457  zold(l) = zetd(l+1) + (prj(l+1)-prh(l))*cpbg*poid(l)
458  zol(l) = zet(l+1) + (prj(l+1)-prh(l))*poi(l)*cpbg
459  END DO
460  ehtd = 0.0_8
461  hccd = 0.0_8
462  f2d = 0.0_8
463  cvwd = 0.0_8
464  ucud = 0.0_8
465  qhtd = 0.0_8
466  bk2d = 0.0_8
467  rasald = 0.0_8
468  hcldd = 0.0_8
469  etad = 0.0_8
470  gmhd = 0.0_8
471  rnnd = 0.0_8
472  gmsd = 0.0_8
473  updfrcd = 0.0_8
474  rnsd = 0.0_8
475  vcud = 0.0_8
476  rmfdd = 0.0_8
477  updfrpd = 0.0_8
478  cll0d = 0.0_8
479  clld = 0.0_8
480  rmfpd = 0.0_8
481  DO ic=k,icmin+1,-1
482  ucud(icmin:) = 0.0_8
483  ucu(icmin:) = 0.
484  vcud(icmin:) = 0.0_8
485  vcu(icmin:) = 0.
486  alm = 0.
487  IF (1. .GT. (qoi(k)/qst(k)-rhmn)/(rhmx-rhmn)) THEN
488  trgd = (qoid(k)*qst(k)-qoi(k)*qstd(k))/qst(k)**2/(rhmx-rhmn)
489  trg = (qoi(k)/qst(k)-rhmn)/(rhmx-rhmn)
490  ELSE
491  trg = 1.
492  trgd = 0.0_8
493  END IF
494  IF (0.0 .LT. (autorampb-sige(ic))/0.2) THEN
495  y1 = (autorampb-sige(ic))/0.2
496  ELSE
497  y1 = 0.0
498  END IF
499  IF (1.0 .GT. y1) THEN
500  f4 = y1
501  ELSE
502  f4 = 1.0
503  END IF
504  IF (trg .GT. 1.0e-5) THEN
505 !================>>
506 !RECOMPUTE SOUNDING UP TO DETRAINMENT LEVEL
507  poi_cd = poid
508  poi_c = poi
509  qoi_cd = qoid
510  qoi_c = qoi
511  poi_cd(k) = poi_cd(k) + tpertd
512  poi_c(k) = poi_c(k) + tpert
513  qoi_c(k) = qoi_c(k) + qpert
514  zetd(k+1) = 0.0_8
515  zet(k+1) = 0.
516  shtd(k+1) = cp*prj(k+1)*poi_cd(k)
517  sht(k+1) = cp*poi_c(k)*prj(k+1)
518  DO l=k,ic,-1
519  IF (qst(l)*rhmax .GT. qoi_c(l)) THEN
520  qold(l) = qoi_cd(l)
521  qol(l) = qoi_c(l)
522  ELSE
523  qold(l) = rhmax*qstd(l)
524  qol(l) = qst(l)*rhmax
525  END IF
526  IF (0.000 .LT. qol(l)) THEN
527  qol(l) = qol(l)
528  ELSE
529  qold(l) = 0.0_8
530  qol(l) = 0.000
531  END IF
532  ssld(l) = cp*prj(l+1)*poi_cd(l) + grav*zetd(l+1)
533  ssl(l) = cp*prj(l+1)*poi_c(l) + grav*zet(l+1)
534  hold(l) = ssld(l) + alhl*qold(l)
535  hol(l) = ssl(l) + qol(l)*alhl
536  hstd(l) = ssld(l) + alhl*qstd(l)
537  hst(l) = ssl(l) + qst(l)*alhl
538  temd = (prj(l+1)-prj(l))*cpbg*poi_cd(l)
539  tem = poi_c(l)*(prj(l+1)-prj(l))*cpbg
540  zetd(l) = zetd(l+1) + temd
541  zet(l) = zet(l+1) + tem
542  zold(l) = zetd(l+1) + (prj(l+1)-prh(l))*cpbg*poi_cd(l)
543  zol(l) = zet(l+1) + (prj(l+1)-prh(l))*poi_c(l)*cpbg
544  END DO
545  DO l=ic+1,k
546  tem = (prj(l)-prh(l-1))/(prh(l)-prh(l-1))
547  shtd(l) = ssld(l-1) + tem*(ssld(l)-ssld(l-1))
548  sht(l) = ssl(l-1) + tem*(ssl(l)-ssl(l-1))
549  qhtd(l) = .5*(qold(l)+qold(l-1))
550  qht(l) = .5*(qol(l)+qol(l-1))
551  END DO
552 !CALCULATE LAMBDA, ETA, AND WORKFUNCTION
553  lambda_min = .2/mxdiam
554  lambda_max = .2/200.
555  IF (hol(k) .GT. hst(ic)) THEN
556 !================>>
557 !LAMBDA CALCULATION: MS-A18
558  temd = (hstd(ic)-hold(ic))*(zol(ic)-zet(ic+1)) + (hst(ic)-hol(&
559 & ic))*(zold(ic)-zetd(ic+1))
560  tem = (hst(ic)-hol(ic))*(zol(ic)-zet(ic+1))
561  DO l=ic+1,k-1
562  temd = temd + (hstd(ic)-hold(l))*(zet(l)-zet(l+1)) + (hst(ic&
563 & )-hol(l))*(zetd(l)-zetd(l+1))
564  tem = tem + (hst(ic)-hol(l))*(zet(l)-zet(l+1))
565  END DO
566  IF (tem .GT. 0.0) THEN
567 !================>>
568  almd = ((hold(k)-hstd(ic))*tem-(hol(k)-hst(ic))*temd)/tem**2
569  alm = (hol(k)-hst(ic))/tem
570  IF (alm .LE. lambda_max) THEN
571 !================>>
572  toki = 1.0
573  IF (alm .LT. lambda_min) THEN
574  tokid = 2*alm*almd/lambda_min**2
575  toki = (alm/lambda_min)**2
576  ELSE
577  tokid = 0.0_8
578  END IF
579 !ETA CALCULATION: MS-A2
580  DO l=ic+1,k
581  etad(l) = almd*(zet(l)-zet(k)) + alm*(zetd(l)-zetd(k))
582  eta(l) = 1.0 + alm*(zet(l)-zet(k))
583  END DO
584  etad(ic) = almd*(zol(ic)-zet(k)) + alm*(zold(ic)-zetd(k))
585  eta(ic) = 1.0 + alm*(zol(ic)-zet(k))
586 !WORKFUNCTION CALCULATION: MS-A22
587  wfn = 0.0
588  hccd(k) = hold(k)
589  hcc(k) = hol(k)
590  wfnd = 0.0_8
591  DO l=k-1,ic+1,-1
592  hccd(l) = hccd(l+1) + (etad(l)-etad(l+1))*hol(l) + (eta(&
593 & l)-eta(l+1))*hold(l)
594  hcc(l) = hcc(l+1) + (eta(l)-eta(l+1))*hol(l)
595  temd = dpb(l)*hccd(l+1) + dpt(l)*hccd(l)
596  tem = hcc(l+1)*dpb(l) + hcc(l)*dpt(l)
597  ehtd(l) = dpb(l)*etad(l+1) + dpt(l)*etad(l)
598  eht(l) = eta(l+1)*dpb(l) + eta(l)*dpt(l)
599  wfnd = wfnd + (temd-ehtd(l)*hst(l)-eht(l)*hstd(l))*gam(l&
600 & ) + (tem-eht(l)*hst(l))*gamd(l)
601  wfn = wfn + (tem-eht(l)*hst(l))*gam(l)
602  END DO
603  hccd(ic) = hstd(ic)*eta(ic) + hst(ic)*etad(ic)
604  hcc(ic) = hst(ic)*eta(ic)
605  wfnd = wfnd + dpb(ic)*((hccd(ic+1)-hstd(ic)*eta(ic+1)-hst(&
606 & ic)*etad(ic+1))*gam(ic)+(hcc(ic+1)-hst(ic)*eta(ic+1))*&
607 & gamd(ic))
608  wfn = wfn + (hcc(ic+1)-hst(ic)*eta(ic+1))*gam(ic)*dpb(ic)
609 !VERTICAL VELOCITY/KE CALCULATION (ADDED 12/2001 JTB)
610  bk2d(k) = 0.0_8
611  bk2(k) = 0.0
612  bke(k) = 0.0
613  hcldd(k) = hold(k)
614  hcld(k) = hol(k)
615  DO l=k-1,ic,-1
616  hcldd(l) = ((etad(l+1)*hcld(l+1)+eta(l+1)*hcldd(l+1)+(&
617 & etad(l)-etad(l+1))*hol(l)+(eta(l)-eta(l+1))*hold(l))*&
618 & eta(l)-(eta(l+1)*hcld(l+1)+(eta(l)-eta(l+1))*hol(l))*&
619 & etad(l))/eta(l)**2
620  hcld(l) = (eta(l+1)*hcld(l+1)+(eta(l)-eta(l+1))*hol(l))/&
621 & eta(l)
622  temd = (((hcldd(l)-hstd(l))*(zet(l)-zet(l+1))+(hcld(l)-&
623 & hst(l))*(zetd(l)-zetd(l+1)))*(1.0+lbcp*dqq(l))-(hcld(l&
624 & )-hst(l))*(zet(l)-zet(l+1))*lbcp*dqqd(l))/(1.0+lbcp*&
625 & dqq(l))**2
626  tem = (hcld(l)-hst(l))*(zet(l)-zet(l+1))/(1.0+lbcp*dqq(l&
627 & ))
628  bke(l) = bke(l+1) + grav*tem/(cp*prj(l+1)*poi(l))
629  IF (tem .LT. 0.0) THEN
630  max1 = 0.0
631  max1d = 0.0_8
632  ELSE
633  max1d = temd
634  max1 = tem
635  END IF
636  bk2d(l) = bk2d(l+1) + (grav*max1d*cp*prj(l+1)*poi(l)-&
637 & grav*max1*cp*prj(l+1)*poid(l))/(cp*prj(l+1)*poi(l))**2
638  bk2(l) = bk2(l+1) + grav*max1/(cp*prj(l+1)*poi(l))
639  IF (bk2(l) .LT. 0.0) THEN
640  max2 = 0.0
641  max2d = 0.0_8
642  ELSE
643  max2d = bk2d(l)
644  max2 = bk2(l)
645  END IF
646  IF (2.0*max2 .EQ. 0.0) THEN
647  cvwd(l) = 0.0_8
648  ELSE
649  cvwd(l) = max2d/sqrt(2.0*max2)
650  END IF
651  cvw(l) = sqrt(2.0*max2)
652  END DO
653 !ALPHA CALCULATION
654  IF (zet(ic) .LT. 2000.) THEN
655  rasald(ic) = 0.0_8
656  rasal(ic) = rasal1
657  END IF
658  IF (zet(ic) .GE. 2000.) THEN
659  rasald(ic) = (rasal2-rasal1)*zetd(ic)/8000.
660  rasal(ic) = rasal1 + (rasal2-rasal1)*(zet(ic)-2000.)/&
661 & 8000.
662  END IF
663  IF (rasal(ic) .GT. 1.0e5) THEN
664  rasald(ic) = 0.0_8
665  rasal(ic) = 1.0e5
666  ELSE
667  rasal(ic) = rasal(ic)
668  END IF
669  rasald(ic) = -(dt*rasald(ic)/rasal(ic)**2)
670  rasal(ic) = dt/rasal(ic)
671  DO l=ic,k
672  IF (cvw(l) .LT. 1.00) THEN
673  cvwd(l) = 0.0_8
674  cvw(l) = 1.00
675  ELSE
676  cvw(l) = cvw(l)
677  END IF
678  END DO
679  CALL acritn(pol(ic), prs(k), acr, acritfac)
680  IF (wfn .GT. acr) THEN
681 !================>>
682  wlqd = qold(k)
683  wlq = qol(k)
684  uhtd = uoid(k)
685  uht = uoi(k)
686  vhtd = void0(k)
687  vht = voi(k)
688  rnnd(k) = 0.0_8
689  rnn(k) = 0.
690  cll0d(k) = 0.0_8
691  cll0(k) = 0.
692  DO l=k-1,ic,-1
693  temd = etad(l) - etad(l+1)
694  tem = eta(l) - eta(l+1)
695  wlqd = wlqd + temd*qol(l) + tem*qold(l)
696  wlq = wlq + tem*qol(l)
697  uhtd = uhtd + temd*uoi(l) + tem*uoid(l)
698  uht = uht + tem*uoi(l)
699  vhtd = vhtd + temd*voi(l) + tem*void0(l)
700  vht = vht + tem*voi(l)
701  IF (l .GT. ic) THEN
702  tx2d = 0.5*((qstd(l)+qstd(l-1))*eta(l)+(qst(l)+qst(l&
703 & -1))*etad(l))
704  tx2 = 0.5*(qst(l)+qst(l-1))*eta(l)
705  tx3d = 0.5*((hstd(l)+hstd(l-1))*eta(l)+(hst(l)+hst(l&
706 & -1))*etad(l))
707  tx3 = 0.5*(hst(l)+hst(l-1))*eta(l)
708  qccd = tx2d + gm1d(l)*(hcc(l)-tx3) + gm1(l)*(hccd(l)&
709 & -tx3d)
710  qcc = tx2 + gm1(l)*(hcc(l)-tx3)
711  cll0d(l) = wlqd - qccd
712  cll0(l) = wlq - qcc
713  ELSE
714  cll0d(l) = wlqd - qstd(ic)*eta(ic) - qst(ic)*etad(ic&
715 & )
716  cll0(l) = wlq - qst(ic)*eta(ic)
717  END IF
718  IF (cll0(l) .LT. 0.00) THEN
719  cll0d(l) = 0.0_8
720  cll0(l) = 0.00
721  ELSE
722  cll0(l) = cll0(l)
723  END IF
724  clid = (cll0d(l)*eta(l)-cll0(l)*etad(l))/eta(l)**2
725  cli = cll0(l)/eta(l)
726  te_ad = prh(l)*poid(l)
727  te_a = poi(l)*prh(l)
728  CALL sundq3_ice_d(te_a, te_ad, sdqv2, sdqv3, sdqvt1, &
729 & f2, f2d, f3)
730  c00_xd = co_auto(i)*f3*f4*f2d
731  c00_x = co_auto(i)*f2*f3*f4
732  cli_crit_xd = -(cli_crit*f3*f2d/(f2*f3)**2)
733  cli_crit_x = cli_crit/(f2*f3)
734  arg1d = -((2*cli*clid*cli_crit_x**2-cli**2*2*&
735 & cli_crit_x*cli_crit_xd)/(cli_crit_x**2)**2)
736  arg1 = -(cli**2/cli_crit_x**2)
737  rated = c00_xd*(1.0-exp(arg1)) - c00_x*arg1d*exp(arg1)
738  rate = c00_x*(1.0-exp(arg1))
739  IF (cvw(l) .LT. 1.00) THEN
740  cvw_x = 1.00
741  cvw_xd = 0.0_8
742  ELSE
743  cvw_xd = cvwd(l)
744  cvw_x = cvw(l)
745  END IF
746  dt_lyrd = ((zetd(l)-zetd(l+1))*cvw_x-(zet(l)-zet(l+1))&
747 & *cvw_xd)/cvw_x**2
748  dt_lyr = (zet(l)-zet(l+1))/cvw_x
749  clossd = cll0d(l)*rate*dt_lyr + cll0(l)*(rated*dt_lyr+&
750 & rate*dt_lyrd)
751  closs = cll0(l)*rate*dt_lyr
752  IF (closs .GT. cll0(l)) THEN
753  clossd = cll0d(l)
754  closs = cll0(l)
755  ELSE
756  closs = closs
757  END IF
758  cll0d(l) = cll0d(l) - clossd
759  cll0(l) = cll0(l) - closs
760  dll0(l) = closs
761  IF (closs .GT. 0.) THEN
762  wlqd = wlqd - clossd
763  wlq = wlq - closs
764  rnnd(l) = clossd
765  rnn(l) = closs
766  ELSE
767  rnnd(l) = 0.0_8
768  rnn(l) = 0.
769  END IF
770  END DO
771  wlqd = wlqd - qstd(ic)*eta(ic) - qst(ic)*etad(ic)
772  wlq = wlq - qst(ic)*eta(ic)
773 !CALCULATE GAMMAS AND KERNEL
774  gmsd(k) = pri(k)*(shtd(k)-ssld(k))
775  gms(k) = (sht(k)-ssl(k))*pri(k)
776  gmhd(k) = gmsd(k) + pri(k)*alhl*(qhtd(k)-qold(k))
777  gmh(k) = gms(k) + (qht(k)-qol(k))*pri(k)*alhl
778  akmd = dpb(k-1)*(gmhd(k)*gam(k-1)+gmh(k)*gamd(k-1))
779  akm = gmh(k)*gam(k-1)*dpb(k-1)
780  tx2d = gmhd(k)
781  tx2 = gmh(k)
782  DO l=k-1,ic+1,-1
783  gmsd(l) = pri(l)*(etad(l)*(sht(l)-ssl(l))+eta(l)*(shtd&
784 & (l)-ssld(l))+etad(l+1)*(ssl(l)-sht(l+1))+eta(l+1)*(&
785 & ssld(l)-shtd(l+1)))
786  gms(l) = (eta(l)*(sht(l)-ssl(l))+eta(l+1)*(ssl(l)-sht(&
787 & l+1)))*pri(l)
788  gmhd(l) = gmsd(l) + alhl*pri(l)*(etad(l)*(qht(l)-qol(l&
789 & ))+eta(l)*(qhtd(l)-qold(l))+etad(l+1)*(qol(l)-qht(l+&
790 & 1))+eta(l+1)*(qold(l)-qhtd(l+1)))
791  gmh(l) = gms(l) + (eta(l)*(qht(l)-qol(l))+eta(l+1)*(&
792 & qol(l)-qht(l+1)))*alhl*pri(l)
793  tx2d = tx2d + (etad(l)-etad(l+1))*gmh(l) + (eta(l)-eta&
794 & (l+1))*gmhd(l)
795  tx2 = tx2 + (eta(l)-eta(l+1))*gmh(l)
796  akmd = akmd - pki(l)*(gmsd(l)*eht(l)+gms(l)*ehtd(l)) +&
797 & tx2d*ght(l) + tx2*ghtd(l)
798  akm = akm - gms(l)*eht(l)*pki(l) + tx2*ght(l)
799  END DO
800  gmsd(ic) = pri(ic)*(etad(ic+1)*(ssl(ic)-sht(ic+1))+eta(&
801 & ic+1)*(ssld(ic)-shtd(ic+1)))
802  gms(ic) = eta(ic+1)*(ssl(ic)-sht(ic+1))*pri(ic)
803  akmd = akmd - dpb(ic)*pki(ic)*(gmsd(ic)*eta(ic+1)+gms(ic&
804 & )*etad(ic+1))
805  akm = akm - gms(ic)*eta(ic+1)*dpb(ic)*pki(ic)
806  gmhd(ic) = gmsd(ic) + pri(ic)*(alhl*(etad(ic+1)*(qol(ic)&
807 & -qht(ic+1))+eta(ic+1)*(qold(ic)-qhtd(ic+1)))+etad(ic)*&
808 & (hst(ic)-hol(ic))+eta(ic)*(hstd(ic)-hold(ic)))
809  gmh(ic) = gms(ic) + (eta(ic+1)*(qol(ic)-qht(ic+1))*alhl+&
810 & eta(ic)*(hst(ic)-hol(ic)))*pri(ic)
811 !CLOUD BASE MASS FLUX
812  IF (.NOT.(akm .GE. 0.0 .OR. wlq .LT. 0.0)) THEN
813 !================>>
814  wfnd = -((wfnd*akm-(wfn-acr)*akmd)/akm**2)
815  wfn = -((wfn-acr)/akm)
816  x1d = (rasald(ic)*wfn+rasal(ic)*wfnd)*trg*toki + rasal&
817 & (ic)*wfn*(trgd*toki+trg*tokid)
818  x1 = rasal(ic)*trg*toki*wfn
819  IF (x1 .GT. (prs(k+1)-prs(k))*(100.*pblfrac)) THEN
820  wfn = (prs(k+1)-prs(k))*(100.*pblfrac)
821  wfnd = 0.0_8
822  ELSE
823  wfnd = x1d
824  wfn = x1
825  END IF
826 !CUMULATIVE PRECIP AND CLOUD-BASE MASS FLUX FOR OUTPUT
827  temd = gravi*wfnd
828  tem = wfn*gravi
829  clld(ic) = clld(ic) + wlqd*tem + wlq*temd
830  cll(ic) = cll(ic) + wlq*tem
831  rmf(ic) = rmf(ic) + tem
832  rmfdd(ic) = rmfdd(ic) + temd*eta(ic) + tem*etad(ic)
833  rmfd(ic) = rmfd(ic) + tem*eta(ic)
834  DO l=ic+1,k
835  rmfpd(l) = temd*eta(l) + tem*etad(l)
836  rmfp(l) = tem*eta(l)
837  rmfc(l) = rmfc(l) + rmfp(l)
838  dllx(l) = dllx(l) + tem*dll0(l)
839  IF (cvw(l) .GT. 0.0) THEN
840  updfrpd(l) = (ddt*1000.*rmfpd(l)*cvw(l)*prs(l)/&
841 & daylen-rmfp(l)*ddt*1000.*prs(l)*cvwd(l)/daylen)/&
842 & (cvw(l)*prs(l))**2
843  updfrp(l) = rmfp(l)*(ddt/daylen)*1000./(cvw(l)*prs&
844 & (l))
845  ELSE
846  updfrpd(l) = 0.0_8
847  updfrp(l) = 0.0
848  END IF
849  clli(l) = cll0(l)/eta(l)
850  updfrcd(l) = updfrcd(l) + updfrpd(l)
851  updfrc(l) = updfrc(l) + updfrp(l)
852  END DO
853 !THETA AND Q CHANGE DUE TO CLOUD TYPE IC
854  DO l=ic,k
855  rnsd(l) = rnsd(l) + rnnd(l)*tem + rnn(l)*temd
856  rns(l) = rns(l) + rnn(l)*tem
857  gmhd(l) = gmhd(l)*wfn + gmh(l)*wfnd
858  gmh(l) = gmh(l)*wfn
859  gmsd(l) = gmsd(l)*wfn + gms(l)*wfnd
860  gms(l) = gms(l)*wfn
861  qoid(l) = qoid(l) + alhi*(gmhd(l)-gmsd(l))
862  qoi(l) = qoi(l) + (gmh(l)-gms(l))*alhi
863  poid(l) = poid(l) + pki(l)*cpi*gmsd(l)
864  poi(l) = poi(l) + gms(l)*pki(l)*cpi
865  qstd(l) = qstd(l) + cpi*(gmsd(l)*bet(l)+gms(l)*betd(&
866 & l))
867  qst(l) = qst(l) + gms(l)*bet(l)*cpi
868  END DO
869 !*FRICFAC*0.5
870  wfnd = 0.5*wfnd
871  wfn = wfn*0.5*1.0
872 !CUMULUS FRICTION
873  IF (fricfac .GT. 0.0) THEN
874 !================>>
875  wfnd = fricfac*(wfnd*exp(-(alm/friclambda))-wfn*almd&
876 & *exp(-(alm/friclambda))/friclambda)
877  wfn = wfn*fricfac*exp(-(alm/friclambda))
878  temd = pri(k)*wfnd
879  tem = wfn*pri(k)
880  ucud(k) = ucud(k) + temd*(uoi(k-1)-uoi(k)) + tem*(&
881 & uoid(k-1)-uoid(k))
882  ucu(k) = ucu(k) + tem*(uoi(k-1)-uoi(k))
883  vcud(k) = vcud(k) + temd*(voi(k-1)-voi(k)) + tem*(&
884 & void0(k-1)-void0(k))
885  vcu(k) = vcu(k) + tem*(voi(k-1)-voi(k))
886  DO l=k-1,ic+1,-1
887  temd = pri(l)*wfnd
888  tem = wfn*pri(l)
889  ucud(l) = ucud(l) + temd*((uoi(l-1)-uoi(l))*eta(l)&
890 & +(uoi(l)-uoi(l+1))*eta(l+1)) + tem*((uoid(l-1)-&
891 & uoid(l))*eta(l)+(uoi(l-1)-uoi(l))*etad(l)+(uoid(&
892 & l)-uoid(l+1))*eta(l+1)+(uoi(l)-uoi(l+1))*etad(l+&
893 & 1))
894  ucu(l) = ucu(l) + tem*((uoi(l-1)-uoi(l))*eta(l)+(&
895 & uoi(l)-uoi(l+1))*eta(l+1))
896  vcud(l) = vcud(l) + temd*((voi(l-1)-voi(l))*eta(l)&
897 & +(voi(l)-voi(l+1))*eta(l+1)) + tem*((void0(l-1)-&
898 & void0(l))*eta(l)+(voi(l-1)-voi(l))*etad(l)+(&
899 & void0(l)-void0(l+1))*eta(l+1)+(voi(l)-voi(l+1))*&
900 & etad(l+1))
901  vcu(l) = vcu(l) + tem*((voi(l-1)-voi(l))*eta(l)+(&
902 & voi(l)-voi(l+1))*eta(l+1))
903  END DO
904  temd = pri(ic)*wfnd
905  tem = wfn*pri(ic)
906  ucud(ic) = ucud(ic) + (2.*(uhtd-uoid(ic)*(eta(ic)-&
907 & eta(ic+1))-uoi(ic)*(etad(ic)-etad(ic+1)))-(uoid(ic&
908 & )+uoid(ic+1))*eta(ic+1)-(uoi(ic)+uoi(ic+1))*etad(&
909 & ic+1))*tem + (2.*(uht-uoi(ic)*(eta(ic)-eta(ic+1)))&
910 & -(uoi(ic)+uoi(ic+1))*eta(ic+1))*temd
911  ucu(ic) = ucu(ic) + (2.*(uht-uoi(ic)*(eta(ic)-eta(ic&
912 & +1)))-(uoi(ic)+uoi(ic+1))*eta(ic+1))*tem
913  vcud(ic) = vcud(ic) + (2.*(vhtd-void0(ic)*(eta(ic)-&
914 & eta(ic+1))-voi(ic)*(etad(ic)-etad(ic+1)))-(void0(&
915 & ic)+void0(ic+1))*eta(ic+1)-(voi(ic)+voi(ic+1))*&
916 & etad(ic+1))*tem + (2.*(vht-voi(ic)*(eta(ic)-eta(ic&
917 & +1)))-(voi(ic)+voi(ic+1))*eta(ic+1))*temd
918  vcu(ic) = vcu(ic) + (2.*(vht-voi(ic)*(eta(ic)-eta(ic&
919 & +1)))-(voi(ic)+voi(ic+1))*eta(ic+1))*tem
920  DO l=ic,k
921  uoid(l) = uoid(l) + ucud(l)
922  uoi(l) = uoi(l) + ucu(l)
923  void0(l) = void0(l) + vcud(l)
924  voi(l) = voi(l) + vcu(l)
925  END DO
926  END IF
927  END IF
928  END IF
929  END IF
930  END IF
931  END IF
932  END IF
933  END DO
934 !CLOUD LOOP
935  IF (sum(rmf(icmin:k)) .GT. 0.0) THEN
936  cnv_prc3d = 0.0_8
937  DO l=icmin,k
938  tem = pri(l)*grav
939  cnv_prc3d(i, l) = tem*rnsd(l)
940  cnv_prc3(i, l) = rns(l)*tem
941  END DO
942  thod(i, icmin:k-1) = poid(icmin:k-1)
943  tho(i, icmin:k-1) = poi(icmin:k-1)
944  qhod(i, icmin:k-1) = qoid(icmin:k-1)
945  qho(i, icmin:k-1) = qoi(icmin:k-1)
946  uhod(i, icmin:k-1) = uoid(icmin:k-1)
947  uho(i, icmin:k-1) = uoi(icmin:k-1)
948  vhod(i, icmin:k-1) = void0(icmin:k-1)
949  vho(i, icmin:k-1) = voi(icmin:k-1)
950  cnv_updfrcd = 0.0_8
951  cnv_updfrcd(i, icmin:k-1) = updfrcd(icmin:k-1)
952  cnv_updfrc(i, icmin:k-1) = updfrc(icmin:k-1)
953 !De-strap tendencies from RAS
954  wght = wgt1(i, :)
955 !Scale properly by layer masses
956  wght0 = 0.
957  DO l=k,k0
958  wght0 = wght0 + wght(l)*(ple(i, l+1)-ple(i, l))
959  END DO
960  wght0 = (prs(k+1)-prs(k))/wght0
961  wght = wght0*wght
962  DO l=k,k0
963  thod(i, l) = thod(i, l) + wght(l)*(poid(k)-poi_svd(k))
964  tho(i, l) = tho(i, l) + wght(l)*(poi(k)-poi_sv(k))
965  qhod(i, l) = qhod(i, l) + wght(l)*(qoid(k)-qoi_svd(k))
966  qho(i, l) = qho(i, l) + wght(l)*(qoi(k)-qoi_sv(k))
967  uhod(i, l) = uhod(i, l) + wght(l)*(uoid(k)-uoi_svd(k))
968  uho(i, l) = uho(i, l) + wght(l)*(uoi(k)-uoi_sv(k))
969  vhod(i, l) = vhod(i, l) + wght(l)*(void0(k)-voi_svd(k))
970  vho(i, l) = vho(i, l) + wght(l)*(voi(k)-voi_sv(k))
971  END DO
972  flxdd = 0.0_8
973  flxdd(i, icmin:k) = ddt*rmfdd(icmin:k)/daylen
974  flxd(i, icmin:k) = rmfd(icmin:k)*ddt/daylen
975  clwd = 0.0_8
976  clwd(i, icmin:k) = ddt*clld(icmin:k)/daylen
977  clw(i, icmin:k) = cll(icmin:k)*ddt/daylen
978  flxdd(i, 1:icmin-1) = 0.0_8
979  flxd(i, 1:icmin-1) = 0.
980  clwd(i, 1:icmin-1) = 0.0_8
981  clw(i, 1:icmin-1) = 0.
982  IF (k .LT. k0) THEN
983  flxdd(i, k:k0) = 0.0_8
984  flxd(i, k:k0) = 0.
985  clwd(i, k:k0) = 0.0_8
986  clw(i, k:k0) = 0.
987  END IF
988  ELSE
989  flxdd(i, :) = 0.0_8
990  flxd(i, :) = 0.
991  clwd(i, :) = 0.0_8
992  clw(i, :) = 0.
993  clwd = 0.0_8
994  cnv_prc3d = 0.0_8
995  cnv_updfrcd = 0.0_8
996  flxdd = 0.0_8
997  END IF
998  ELSE
999  flxdd(i, :) = 0.0_8
1000  flxd(i, :) = 0.
1001  clwd(i, :) = 0.0_8
1002  clw(i, :) = 0.
1003  clwd = 0.0_8
1004  cnv_prc3d = 0.0_8
1005  cnv_updfrcd = 0.0_8
1006  flxdd = 0.0_8
1007  END IF
1008 END SUBROUTINE rase_d
1009 
1010 ! Differentiation of sundq3_ice in forward (tangent) mode:
1011 ! variations of useful results: f2
1012 ! with respect to varying inputs: temp f2
1013 SUBROUTINE sundq3_ice_d(temp, tempd, rate2, rate3, te1, f2, f2d, f3)
1014  IMPLICIT NONE
1015  REAL*8, INTENT(IN) :: temp, rate2, rate3, te1
1016  REAL*8, INTENT(IN) :: tempd
1017  REAL*8, INTENT(OUT) :: f2, f3
1018  REAL*8, INTENT(OUT) :: f2d
1019 !,RATE2,RATE3,TE1
1020  REAL*8 :: xx, yy, te0, te2, jump1
1021  te0 = 273.
1022  te2 = 200.
1023  jump1 = (rate2-1.0)/(te0-te1)**0.333
1024 ! Ice - phase treatment
1025  IF (temp .GE. te0) THEN
1026  f2 = 1.0
1027  f3 = 1.0
1028  f2d = 0.0_8
1029  END IF
1030  IF (temp .GE. te1 .AND. temp .LT. te0) THEN
1031  f2d = -(jump1*0.3333*(te0-temp)**(-0.6667)*tempd)
1032  f2 = 1.0 + jump1*(te0-temp)**0.3333
1033  f3 = 1.0
1034  END IF
1035  IF (temp .LT. te1) THEN
1036  f2d = (-((rate3-rate2)*tempd))/(te1-te2)
1037  f2 = rate2 + (rate3-rate2)*(te1-temp)/(te1-te2)
1038  f3 = 1.0
1039  END IF
1040  IF (f2 .GT. 27.0) THEN
1041  f2 = 27.0
1042  f2d = 0.0_8
1043  END IF
1044 END SUBROUTINE sundq3_ice_d
1045 
1046 ! Differentiation of dqsat_ras in forward (tangent) mode:
1047 ! variations of useful results: dqsi qssi
1048 ! with respect to varying inputs: temp
1049 SUBROUTINE dqsat_ras_d(dqsi, dqsid, qssi, qssid, temp, tempd, plo, lm, &
1050 & estblx, cons_h2omw, cons_airmw)
1051  IMPLICIT NONE
1052 !Inputs
1053  INTEGER :: lm
1054  REAL*8, DIMENSION(lm) :: temp, plo
1055  REAL*8, DIMENSION(lm) :: tempd
1056  REAL*8 :: estblx(:)
1057  REAL*8 :: cons_h2omw, cons_airmw
1058 !Outputs
1059  REAL*8, DIMENSION(lm) :: dqsi, qssi
1060  REAL*8, DIMENSION(lm) :: dqsid, qssid
1061 !Locals
1062  REAL*8, PARAMETER :: max_mixing_ratio=1.0
1063  REAL*8 :: esfac
1064  INTEGER :: k
1065  REAL*8 :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
1066  REAL*8 :: tld, ttd, tid, dqsatd, qsatd, qqd, ddd
1067  INTEGER :: it
1068  INTEGER, PARAMETER :: degsubs=100
1069  REAL*8, PARAMETER :: tmintbl=150.0, tmaxtbl=333.0
1070  INTEGER, PARAMETER :: tablesize=nint(tmaxtbl-tmintbl)*degsubs+1
1071  INTRINSIC nint
1072  INTRINSIC int
1073  esfac = cons_h2omw/cons_airmw
1074  dqsid = 0.0_8
1075  qssid = 0.0_8
1076  DO k=1,lm
1077  tld = tempd(k)
1078  tl = temp(k)
1079  pl = plo(k)
1080  pp = pl*100.0
1081  IF (tl .LE. tmintbl) THEN
1082  ti = tmintbl
1083  tid = 0.0_8
1084  ELSE IF (tl .GE. tmaxtbl - .001) THEN
1085  tid = 0.0_8
1086  ti = tmaxtbl - .001
1087  ELSE
1088  tid = tld
1089  ti = tl
1090  END IF
1091  ttd = degsubs*tid
1092  tt = (ti-tmintbl)*degsubs + 1
1093  it = int(tt)
1094  dqq = estblx(it+1) - estblx(it)
1095  qqd = dqq*ttd
1096  qq = (tt-it)*dqq + estblx(it)
1097  IF (pp .LE. qq) THEN
1098  qsat = max_mixing_ratio
1099  dqsat = 0.0
1100  qsatd = 0.0_8
1101  dqsatd = 0.0_8
1102  ELSE
1103  ddd = -((-((1.0-esfac)*qqd))/(pp-(1.0-esfac)*qq)**2)
1104  dd = 1.0/(pp-(1.0-esfac)*qq)
1105  qsatd = esfac*(qqd*dd+qq*ddd)
1106  qsat = esfac*qq*dd
1107  dqsatd = esfac*degsubs*dqq*pp*(ddd*dd+dd*ddd)
1108  dqsat = esfac*degsubs*dqq*pp*(dd*dd)
1109  END IF
1110  dqsid(k) = dqsatd
1111  dqsi(k) = dqsat
1112  qssid(k) = qsatd
1113  qssi(k) = qsat
1114  END DO
1115 END SUBROUTINE dqsat_ras_d
1116 
1117 ! Differentiation of dqsats_ras in forward (tangent) mode:
1118 ! variations of useful results: dqsi qssi
1119 ! with respect to varying inputs: temp
1120 SUBROUTINE dqsats_ras_d(dqsi, dqsid, qssi, qssid, temp, tempd, plo, &
1121 & estblx, cons_h2omw, cons_airmw)
1122  IMPLICIT NONE
1123 !Inputs
1124  REAL*8 :: temp, plo
1125  REAL*8 :: tempd
1126  REAL*8 :: estblx(:)
1127  REAL*8 :: cons_h2omw, cons_airmw
1128 !Outputs
1129  REAL*8 :: dqsi, qssi
1130  REAL*8 :: dqsid, qssid
1131 !Locals
1132  REAL*8, PARAMETER :: max_mixing_ratio=1.0
1133  REAL*8 :: esfac
1134  REAL*8 :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
1135  REAL*8 :: tld, ttd, tid, dqsatd, qsatd, qqd, ddd
1136  INTEGER :: it
1137  INTEGER, PARAMETER :: degsubs=100
1138  REAL*8, PARAMETER :: tmintbl=150.0, tmaxtbl=333.0
1139  INTEGER, PARAMETER :: tablesize=nint(tmaxtbl-tmintbl)*degsubs+1
1140  INTRINSIC nint
1141  INTRINSIC int
1142  esfac = cons_h2omw/cons_airmw
1143  tld = tempd
1144  tl = temp
1145  pl = plo
1146  pp = pl*100.0
1147  IF (tl .LE. tmintbl) THEN
1148  ti = tmintbl
1149  tid = 0.0_8
1150  ELSE IF (tl .GE. tmaxtbl - .001) THEN
1151  ti = tmaxtbl - .001
1152  tid = 0.0_8
1153  ELSE
1154  tid = tld
1155  ti = tl
1156  END IF
1157  ttd = degsubs*tid
1158  tt = (ti-tmintbl)*degsubs + 1
1159  it = int(tt)
1160  dqq = estblx(it+1) - estblx(it)
1161  qqd = dqq*ttd
1162  qq = (tt-it)*dqq + estblx(it)
1163  IF (pp .LE. qq) THEN
1164  qsat = max_mixing_ratio
1165  dqsat = 0.0
1166  qsatd = 0.0_8
1167  dqsatd = 0.0_8
1168  ELSE
1169  ddd = -((-((1.0-esfac)*qqd))/(pp-(1.0-esfac)*qq)**2)
1170  dd = 1.0/(pp-(1.0-esfac)*qq)
1171  qsatd = esfac*(qqd*dd+qq*ddd)
1172  qsat = esfac*qq*dd
1173  dqsatd = esfac*degsubs*dqq*pp*(ddd*dd+dd*ddd)
1174  dqsat = esfac*degsubs*dqq*pp*(dd*dd)
1175  END IF
1176  dqsid = dqsatd
1177  dqsi = dqsat
1178  qssid = qsatd
1179  qssi = qsat
1180 END SUBROUTINE dqsats_ras_d
1181 
1182 
1183 SUBROUTINE rase0_d(idim, irun, k0, icmin, dt, cons_cp, cons_alhl, &
1184 & cons_grav, cons_rgas, cons_h2omw, cons_airmw, cons_vireps, seedras, &
1185 & sige, kcbl, wgt0, wgt1, frland, ts, tho, thod, qho, qhod, co_auto, ple&
1186 & , rasparams, estblx)
1187  IMPLICIT NONE
1188 !INPUTS
1189  INTEGER, INTENT(IN) :: idim, irun, k0, icmin
1190  REAL*8, DIMENSION(idim, k0 + 1), INTENT(IN) :: ple
1191  REAL*8, DIMENSION(k0 + 1), INTENT(IN) :: sige
1192  REAL*8, INTENT(IN) :: dt, cons_cp, cons_alhl, cons_grav, cons_rgas
1193  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_vireps
1194  INTEGER, DIMENSION(idim), INTENT(IN) :: seedras
1195  INTEGER, DIMENSION(idim), INTENT(IN) :: kcbl
1196  REAL*8, DIMENSION(idim), INTENT(IN) :: ts, frland
1197  REAL*8, DIMENSION(idim), INTENT(IN) :: co_auto
1198  REAL*8, DIMENSION(idim, k0), INTENT(IN) :: wgt0, wgt1
1199  REAL*8, DIMENSION(:), INTENT(IN) :: rasparams
1200  REAL*8, DIMENSION(:), INTENT(IN) :: estblx
1201 !PROGNOSTIC
1202  REAL*8, DIMENSION(idim, k0), INTENT(INOUT) :: tho, qho
1203  REAL*8, DIMENSION(idim, k0), INTENT(INOUT) :: thod, qhod
1204 !LOCALS
1205  INTEGER :: i, ic, l, kk, k
1206 !Parameters
1207  REAL*8, PARAMETER :: onepkap=1.+2./7., daylen=86400.0
1208  REAL*8, PARAMETER :: rhmax=0.9999
1209  REAL*8, PARAMETER :: cbl_qpert=0.0, cbl_tpert=1.0
1210  REAL*8, PARAMETER :: cbl_tpert_mxocn=2.0, cbl_tpert_mxlnd=4.0
1211 !Constants
1212  REAL*8 :: grav, cp, alhl, cpbg, alhi, cpi, gravi, ddt, lbcp
1213 !Rasparams
1214  REAL*8 :: fricfac, cli_crit, rasal1, rasal2
1215  REAL*8 :: friclambda
1216  REAL*8 :: sdqv2, sdqv3, sdqvt1
1217  REAL*8 :: acritfac, pblfrac, autorampb
1218  REAL*8 :: maxdallowed, rhmn, rhmx
1219  REAL*8 :: mxdiam
1220  REAL*8 :: tx2, tx3, akm, acr, alm, tth, qqh, dqx
1221  REAL*8 :: tx2d, akmd, almd
1222  REAL*8 :: wfn, tem, trg, trgexp, evp, wlq, qcc
1223  REAL*8 :: wfnd, temd, trgd
1224  REAL*8 :: cli, te_a, c00_x, cli_crit_x, toki
1225  REAL*8 :: tokid
1226  REAL*8 :: dt_lyr, rate, cvw_x, closs, f2, f3, f4
1227  REAL*8 :: wght0, prcbl, rndu
1228  REAL*8 :: lambda_min, lambda_max
1229  REAL*8 :: tpert, qpert
1230  REAL*8 :: tpertd
1231  REAL*8, DIMENSION(k0) :: poi_sv, qoi_sv
1232  REAL*8, DIMENSION(k0) :: poi_svd, qoi_svd
1233  REAL*8, DIMENSION(k0) :: poi, qoi, dqq, bet, gam, cll
1234  REAL*8, DIMENSION(k0) :: poid, qoid, dqqd, betd, gamd
1235  REAL*8, DIMENSION(k0) :: poi_c, qoi_c
1236  REAL*8, DIMENSION(k0) :: poi_cd, qoi_cd
1237  REAL*8, DIMENSION(k0) :: prh, pri, ght, dpt, dpb, pki
1238  REAL*8, DIMENSION(k0) :: prhd, prid, ghtd, dptd, dpbd, pkid
1239  REAL*8, DIMENSION(k0) :: cln, rns, pol, dm
1240  REAL*8, DIMENSION(k0) :: pold
1241  REAL*8, DIMENSION(k0) :: qst, ssl, rmf, rnn, rn1, rmfc, rmfp
1242  REAL*8, DIMENSION(k0) :: qstd, ssld
1243  REAL*8, DIMENSION(k0) :: gms, eta, gmh, eht, gm1, hcc, rmfd
1244  REAL*8, DIMENSION(k0) :: gmsd, etad, gmhd, ehtd, hccd
1245  REAL*8, DIMENSION(k0) :: hol, hst, qol, zol, hcld, cll0, cllx, clli
1246  REAL*8, DIMENSION(k0) :: hold, hstd, qold, zold
1247  REAL*8, DIMENSION(k0) :: bke, cvw, updfrc
1248  REAL*8, DIMENSION(k0) :: rasal, updfrp, bk2, dll0, dllx
1249  REAL*8, DIMENSION(k0) :: rasald
1250  REAL*8, DIMENSION(k0) :: wght, massf
1251  REAL*8, DIMENSION(k0) :: wghtd
1252  REAL*8, DIMENSION(k0) :: qss, dqs, pf, pk, tempf, zlo
1253  REAL*8, DIMENSION(k0) :: qssd, dqsd, tempfd, zlod
1254  REAL*8, DIMENSION(k0 + 1) :: prj, prs, qht, sht, zet, zle, pke
1255  REAL*8, DIMENSION(k0+1) :: prjd, prsd, qhtd, shtd, zetd, zled
1256  INTRINSIC max
1257  INTRINSIC min
1258  INTRINSIC sqrt
1259  INTRINSIC exp
1260  INTRINSIC sum
1261  REAL*8, DIMENSION(k0+1) :: pwx1
1262  REAL*8 :: pwy1
1263  REAL*8, DIMENSION(k0) :: pwx10
1264  REAL*8 :: pwx11
1265  REAL*8 :: pwr1
1266  REAL*8 :: arg1
1267  REAL*8 :: x1
1268  REAL*8 :: x1d
1269  REAL*8 :: max2
1270  REAL*8 :: max1
1271  REAL*8 :: y1
1272 ! --- 1
1273  fricfac = rasparams(1)
1274 ! --- 4
1275  cli_crit = rasparams(4)
1276 ! --- 5
1277  rasal1 = rasparams(5)
1278 ! --- 6
1279  rasal2 = rasparams(6)
1280 ! --- 11
1281  friclambda = rasparams(11)
1282 ! --- 14
1283  sdqv2 = rasparams(14)
1284 ! --- 15
1285  sdqv3 = rasparams(15)
1286 ! --- 16
1287  sdqvt1 = rasparams(16)
1288 ! --- 17
1289  acritfac = rasparams(17)
1290 ! --- 20
1291  pblfrac = rasparams(20)
1292 ! --- 21
1293  autorampb = rasparams(21)
1294 ! --- 24
1295  rhmn = rasparams(24)
1296 ! --- 24
1297  maxdallowed = rasparams(23)
1298 ! --- 25
1299  rhmx = rasparams(25)
1300  grav = cons_grav
1301  alhl = cons_alhl
1302  cp = cons_cp
1303  cpi = 1.0/cp
1304  alhi = 1.0/alhl
1305  gravi = 1.0/grav
1306  cpbg = cp*gravi
1307  ddt = daylen/dt
1308  lbcp = alhl*cpi
1309  ehtd = 0.0_8
1310  hccd = 0.0_8
1311  dqqd = 0.0_8
1312  dqsd = 0.0_8
1313  ghtd = 0.0_8
1314  hstd = 0.0_8
1315  betd = 0.0_8
1316  qhtd = 0.0_8
1317  qold = 0.0_8
1318  rasald = 0.0_8
1319  etad = 0.0_8
1320  shtd = 0.0_8
1321  gmhd = 0.0_8
1322  qssd = 0.0_8
1323  qstd = 0.0_8
1324  gmsd = 0.0_8
1325  ssld = 0.0_8
1326  zetd = 0.0_8
1327  zold = 0.0_8
1328  gamd = 0.0_8
1329  DO i=1,irun
1330 !CALL FINDBASE
1331  k = kcbl(i)
1332  IF (k .GT. 0) THEN
1333 !Get saturation specific humidity and gradient wrt to T
1334  pwx1 = ple(i, :)/1000.
1335  pwy1 = cons_rgas/cons_cp
1336  pke = pwx1**pwy1
1337  pf = 0.5*(ple(i, 1:k0)+ple(i, 2:k0+1))
1338  pwx10 = pf/1000.
1339  pwy1 = cons_rgas/cons_cp
1340  pk = pwx10**pwy1
1341  tempfd = pk*thod(i, :)
1342  tempf = tho(i, :)*pk
1343  zle = 0.0
1344  zlo = 0.0
1345  zled(k0+1) = 0.0_8
1346  zle(k0+1) = 0.
1347  zlod = 0.0_8
1348  zled = 0.0_8
1349  DO l=k0,1,-1
1350  zled(l) = thod(i, l)*(1.+cons_vireps*qho(i, l)) + tho(i, l)*&
1351 & cons_vireps*qhod(i, l)
1352  zle(l) = tho(i, l)*(1.+cons_vireps*qho(i, l))
1353  zlod(l) = zled(l+1) + cons_cp*(pke(l+1)-pk(l))*zled(l)/cons_grav
1354  zlo(l) = zle(l+1) + cons_cp/cons_grav*(pke(l+1)-pk(l))*zle(l)
1355  zled(l) = zlod(l) + cons_cp*(pk(l)-pke(l))*zled(l)/cons_grav
1356  zle(l) = zlo(l) + cons_cp/cons_grav*(pk(l)-pke(l))*zle(l)
1357  END DO
1358  tpertd = cbl_tpert*(-tempfd(k0)-cons_grav*zlod(k0)/cons_cp)
1359  tpert = cbl_tpert*(ts(i)-(tempf(k0)+cons_grav*zlo(k0)/cons_cp))
1360 !* ( QSSFC - Q(:,:,K0) ) [CBL_QPERT = 0.0]
1361  qpert = cbl_qpert
1362  IF (tpert .LT. 0.0) THEN
1363  tpert = 0.0
1364  tpertd = 0.0_8
1365  ELSE
1366  tpert = tpert
1367  END IF
1368  IF (qpert .LT. 0.0) THEN
1369  qpert = 0.0
1370  ELSE
1371  qpert = qpert
1372  END IF
1373  IF (frland(i) .LT. 0.1) THEN
1374  IF (tpert .GT. cbl_tpert_mxocn) THEN
1375  tpert = cbl_tpert_mxocn
1376  tpertd = 0.0_8
1377  ELSE
1378  tpert = tpert
1379  END IF
1380  ELSE IF (tpert .GT. cbl_tpert_mxlnd) THEN
1381  tpert = cbl_tpert_mxlnd
1382  tpertd = 0.0_8
1383  ELSE
1384  tpert = tpert
1385  END IF
1386  CALL dqsat_ras_d(dqs, dqsd, qss, qssd, tempf, tempfd, pf, k0, &
1387 & estblx, cons_h2omw, cons_airmw)
1388  DO kk=icmin,k+1
1389  prjd(kk) = 0.0_8
1390  prj(kk) = pke(kk)
1391  END DO
1392 ! These initialized here in order not to confuse Valgrind debugger
1393  poi = 0.
1394 ! Do not believe it actually makes any difference.
1395  qoi = 0.
1396  prsd(icmin:k0+1) = 0.0_8
1397  prs(icmin:k0+1) = ple(i, icmin:k0+1)
1398  poid = 0.0_8
1399  poid(icmin:k) = thod(i, icmin:k)
1400  poi(icmin:k) = tho(i, icmin:k)
1401  qoid = 0.0_8
1402  qoid(icmin:k) = qhod(i, icmin:k)
1403  qoi(icmin:k) = qho(i, icmin:k)
1404  qstd(icmin:k) = qssd(icmin:k)
1405  qst(icmin:k) = qss(icmin:k)
1406  dqqd(icmin:k) = dqsd(icmin:k)
1407  dqq(icmin:k) = dqs(icmin:k)
1408 !Mass fraction of each layer below cloud base
1409  massf(:) = wgt0(i, :)
1410 !RESET PRESSURE at bottom edge of CBL
1411  prcbl = prs(k)
1412  DO l=k,k0
1413  prcbl = prcbl + massf(l)*(prs(l+1)-prs(l))
1414  END DO
1415  prsd(k+1) = 0.0_8
1416  prs(k+1) = prcbl
1417  pwx11 = prs(k+1)/1000.
1418  pwy1 = cons_rgas/cons_cp
1419  prjd(k+1) = 0.0_8
1420  prj(k+1) = pwx11**pwy1
1421  DO l=k,icmin,-1
1422  pold(l) = 0.0_8
1423  pol(l) = 0.5*(prs(l)+prs(l+1))
1424  prhd(l) = 0.0_8
1425  prh(l) = (prs(l+1)*prj(l+1)-prs(l)*prj(l))/(onepkap*(prs(l+1)-&
1426 & prs(l)))
1427  pkid(l) = 0.0_8
1428  pki(l) = 1.0/prh(l)
1429  dptd(l) = 0.0_8
1430  dpt(l) = prh(l) - prj(l)
1431  dpbd(l) = 0.0_8
1432  dpb(l) = prj(l+1) - prh(l)
1433  prid(l) = 0.0_8
1434  pri(l) = .01/(prs(l+1)-prs(l))
1435  END DO
1436 !RECALCULATE PROFILE QUAN. IN LOWEST STRAPPED LAYER
1437  IF (k .LE. k0) THEN
1438  poid(k) = 0.0_8
1439  poi(k) = 0.
1440  qoid(k) = 0.0_8
1441  qoi(k) = 0.
1442 !SPECIFY WEIGHTS GIVEN TO EACH LAYER WITHIN SUBCLOUD "SUPERLAYER"
1443  wght = 0.
1444  DO l=k,k0
1445  wghtd(l) = 0.0_8
1446  wght(l) = massf(l)*(ple(i, l+1)-ple(i, l))/(prs(k+1)-prs(k))
1447  END DO
1448  DO l=k,k0
1449  poid(k) = poid(k) + wght(l)*thod(i, l)
1450  poi(k) = poi(k) + wght(l)*tho(i, l)
1451  qoid(k) = qoid(k) + wght(l)*qhod(i, l)
1452  qoi(k) = qoi(k) + wght(l)*qho(i, l)
1453  END DO
1454  CALL dqsats_ras_d(dqq(k), dqqd(k), qst(k), qstd(k), poi(k)*prh(k&
1455 & ), prh(k)*poid(k), pol(k), estblx, cons_h2omw, &
1456 & cons_airmw)
1457  END IF
1458  IF (seedras(i)/1000000. .LT. 1e-6) THEN
1459  rndu = 1e-6
1460  ELSE
1461  rndu = seedras(i)/1000000.
1462  END IF
1463  pwr1 = rndu**(-(1./2.))
1464  mxdiam = maxdallowed*pwr1
1465  DO l=k,icmin,-1
1466 !*
1467  betd(l) = pki(l)*dqqd(l)
1468  bet(l) = dqq(l)*pki(l)
1469 !*
1470  gamd(l) = -(pki(l)*lbcp*dqqd(l)/(1.0+lbcp*dqq(l))**2)
1471  gam(l) = pki(l)/(1.0+lbcp*dqq(l))
1472  IF (l .LT. k) THEN
1473  ghtd(l+1) = dpb(l)*gamd(l) + dpt(l+1)*gamd(l+1)
1474  ght(l+1) = gam(l)*dpb(l) + gam(l+1)*dpt(l+1)
1475  gm1(l+1) = 0.5*lbcp*(dqq(l)/(alhl*(1.0+lbcp*dqq(l)))+dqq(l+1)/&
1476 & (alhl*(1.0+lbcp*dqq(l+1))))
1477  END IF
1478  END DO
1479  rns = 0.
1480  cll = 0.
1481  rmf = 0.
1482  rmfd = 0.
1483  rmfc = 0.
1484  rmfp = 0.
1485  cll0 = 0.
1486  dll0 = 0.
1487  cllx = 0.
1488  dllx = 0.
1489  clli = 0.
1490  poi_svd = poid
1491  poi_sv = poi
1492  qoi_svd = qoid
1493  qoi_sv = qoi
1494  cvw = 0.0
1495  updfrc = 0.0
1496  updfrp = 0.0
1497 ! HOL initialized here in order not to confuse Valgrind debugger
1498  hol = 0.
1499  zetd(k+1) = 0.0_8
1500  zet(k+1) = 0
1501  shtd(k+1) = cp*prj(k+1)*poid(k)
1502  sht(k+1) = cp*poi(k)*prj(k+1)
1503  hold = 0.0_8
1504  DO l=k,icmin,-1
1505  IF (qst(l)*rhmax .GT. qoi(l)) THEN
1506  qold(l) = qoid(l)
1507  qol(l) = qoi(l)
1508  ELSE
1509  qold(l) = rhmax*qstd(l)
1510  qol(l) = qst(l)*rhmax
1511  END IF
1512  IF (0.000 .LT. qol(l)) THEN
1513  qol(l) = qol(l)
1514  ELSE
1515  qold(l) = 0.0_8
1516  qol(l) = 0.000
1517  END IF
1518  ssld(l) = cp*prj(l+1)*poid(l) + grav*zetd(l+1)
1519  ssl(l) = cp*prj(l+1)*poi(l) + grav*zet(l+1)
1520  hold(l) = ssld(l) + alhl*qold(l)
1521  hol(l) = ssl(l) + qol(l)*alhl
1522  hstd(l) = ssld(l) + alhl*qstd(l)
1523  hst(l) = ssl(l) + qst(l)*alhl
1524  temd = (prj(l+1)-prj(l))*cpbg*poid(l)
1525  tem = poi(l)*(prj(l+1)-prj(l))*cpbg
1526  zetd(l) = zetd(l+1) + temd
1527  zet(l) = zet(l+1) + tem
1528  zold(l) = zetd(l+1) + (prj(l+1)-prh(l))*cpbg*poid(l)
1529  zol(l) = zet(l+1) + (prj(l+1)-prh(l))*poi(l)*cpbg
1530  END DO
1531  DO ic=k,icmin+1,-1
1532  alm = 0.
1533  IF (1. .GT. (qoi(k)/qst(k)-rhmn)/(rhmx-rhmn)) THEN
1534  trgd = (qoid(k)*qst(k)-qoi(k)*qstd(k))/qst(k)**2/(rhmx-rhmn)
1535  trg = (qoi(k)/qst(k)-rhmn)/(rhmx-rhmn)
1536  ELSE
1537  trg = 1.
1538  trgd = 0.0_8
1539  END IF
1540  IF (0.0 .LT. (autorampb-sige(ic))/0.2) THEN
1541  y1 = (autorampb-sige(ic))/0.2
1542  ELSE
1543  y1 = 0.0
1544  END IF
1545  IF (1.0 .GT. y1) THEN
1546  f4 = y1
1547  ELSE
1548  f4 = 1.0
1549  END IF
1550  IF (trg .GT. 1.0e-5) THEN
1551 !================>>
1552 !RECOMPUTE SOUNDING UP TO DETRAINMENT LEVEL
1553  poi_cd = poid
1554  poi_c = poi
1555  qoi_cd = qoid
1556  qoi_c = qoi
1557  poi_cd(k) = poi_cd(k) + tpertd
1558  poi_c(k) = poi_c(k) + tpert
1559  qoi_c(k) = qoi_c(k) + qpert
1560  zetd(k+1) = 0.0_8
1561  zet(k+1) = 0.
1562  shtd(k+1) = cp*prj(k+1)*poi_cd(k)
1563  sht(k+1) = cp*poi_c(k)*prj(k+1)
1564  DO l=k,ic,-1
1565  IF (qst(l)*rhmax .GT. qoi_c(l)) THEN
1566  qold(l) = qoi_cd(l)
1567  qol(l) = qoi_c(l)
1568  ELSE
1569  qold(l) = rhmax*qstd(l)
1570  qol(l) = qst(l)*rhmax
1571  END IF
1572  IF (0.000 .LT. qol(l)) THEN
1573  qol(l) = qol(l)
1574  ELSE
1575  qold(l) = 0.0_8
1576  qol(l) = 0.000
1577  END IF
1578  ssld(l) = cp*prj(l+1)*poi_cd(l) + grav*zetd(l+1)
1579  ssl(l) = cp*prj(l+1)*poi_c(l) + grav*zet(l+1)
1580  hold(l) = ssld(l) + alhl*qold(l)
1581  hol(l) = ssl(l) + qol(l)*alhl
1582  hstd(l) = ssld(l) + alhl*qstd(l)
1583  hst(l) = ssl(l) + qst(l)*alhl
1584  temd = (prj(l+1)-prj(l))*cpbg*poi_cd(l)
1585  tem = poi_c(l)*(prj(l+1)-prj(l))*cpbg
1586  zetd(l) = zetd(l+1) + temd
1587  zet(l) = zet(l+1) + tem
1588  zold(l) = zetd(l+1) + (prj(l+1)-prh(l))*cpbg*poi_cd(l)
1589  zol(l) = zet(l+1) + (prj(l+1)-prh(l))*poi_c(l)*cpbg
1590  END DO
1591  DO l=ic+1,k
1592  tem = (prj(l)-prh(l-1))/(prh(l)-prh(l-1))
1593  shtd(l) = ssld(l-1) + tem*(ssld(l)-ssld(l-1))
1594  sht(l) = ssl(l-1) + tem*(ssl(l)-ssl(l-1))
1595  qhtd(l) = .5*(qold(l)+qold(l-1))
1596  qht(l) = .5*(qol(l)+qol(l-1))
1597  END DO
1598 !CALCULATE LAMBDA, ETA, AND WORKFUNCTION
1599  lambda_min = .2/mxdiam
1600  lambda_max = .2/200.
1601  IF (hol(k) .GT. hst(ic)) THEN
1602 !================>>
1603 !LAMBDA CALCULATION: MS-A18
1604  temd = (hstd(ic)-hold(ic))*(zol(ic)-zet(ic+1)) + (hst(ic)-&
1605 & hol(ic))*(zold(ic)-zetd(ic+1))
1606  tem = (hst(ic)-hol(ic))*(zol(ic)-zet(ic+1))
1607  DO l=ic+1,k-1
1608  temd = temd + (hstd(ic)-hold(l))*(zet(l)-zet(l+1)) + (hst(&
1609 & ic)-hol(l))*(zetd(l)-zetd(l+1))
1610  tem = tem + (hst(ic)-hol(l))*(zet(l)-zet(l+1))
1611  END DO
1612  IF (tem .GT. 0.0) THEN
1613 !================>>
1614  almd = ((hold(k)-hstd(ic))*tem-(hol(k)-hst(ic))*temd)/tem&
1615 & **2
1616  alm = (hol(k)-hst(ic))/tem
1617  IF (alm .LE. lambda_max) THEN
1618 !================>>
1619  toki = 1.0
1620  IF (alm .LT. lambda_min) THEN
1621  tokid = 2*alm*almd/lambda_min**2
1622  toki = (alm/lambda_min)**2
1623  ELSE
1624  tokid = 0.0_8
1625  END IF
1626 !ETA CALCULATION: MS-A2
1627  DO l=ic+1,k
1628  etad(l) = almd*(zet(l)-zet(k)) + alm*(zetd(l)-zetd(k))
1629  eta(l) = 1.0 + alm*(zet(l)-zet(k))
1630  END DO
1631  etad(ic) = almd*(zol(ic)-zet(k)) + alm*(zold(ic)-zetd(k)&
1632 & )
1633  eta(ic) = 1.0 + alm*(zol(ic)-zet(k))
1634 !WORKFUNCTION CALCULATION: MS-A22
1635  wfn = 0.0
1636  hccd(k) = hold(k)
1637  hcc(k) = hol(k)
1638  wfnd = 0.0_8
1639  DO l=k-1,ic+1,-1
1640  hccd(l) = hccd(l+1) + (etad(l)-etad(l+1))*hol(l) + (&
1641 & eta(l)-eta(l+1))*hold(l)
1642  hcc(l) = hcc(l+1) + (eta(l)-eta(l+1))*hol(l)
1643  temd = dpb(l)*hccd(l+1) + dpt(l)*hccd(l)
1644  tem = hcc(l+1)*dpb(l) + hcc(l)*dpt(l)
1645  ehtd(l) = dpb(l)*etad(l+1) + dpt(l)*etad(l)
1646  eht(l) = eta(l+1)*dpb(l) + eta(l)*dpt(l)
1647  wfnd = wfnd + (temd-ehtd(l)*hst(l)-eht(l)*hstd(l))*gam&
1648 & (l) + (tem-eht(l)*hst(l))*gamd(l)
1649  wfn = wfn + (tem-eht(l)*hst(l))*gam(l)
1650  END DO
1651  hccd(ic) = hstd(ic)*eta(ic) + hst(ic)*etad(ic)
1652  hcc(ic) = hst(ic)*eta(ic)
1653  wfnd = wfnd + dpb(ic)*((hccd(ic+1)-hstd(ic)*eta(ic+1)-&
1654 & hst(ic)*etad(ic+1))*gam(ic)+(hcc(ic+1)-hst(ic)*eta(ic+&
1655 & 1))*gamd(ic))
1656  wfn = wfn + (hcc(ic+1)-hst(ic)*eta(ic+1))*gam(ic)*dpb(ic&
1657 & )
1658 !VERTICAL VELOCITY/KE CALCULATION (ADDED 12/2001 JTB)
1659  bk2(k) = 0.0
1660  bke(k) = 0.0
1661  hcld(k) = hol(k)
1662  DO l=k-1,ic,-1
1663  hcld(l) = (eta(l+1)*hcld(l+1)+(eta(l)-eta(l+1))*hol(l)&
1664 & )/eta(l)
1665  tem = (hcld(l)-hst(l))*(zet(l)-zet(l+1))/(1.0+lbcp*dqq&
1666 & (l))
1667  bke(l) = bke(l+1) + grav*tem/(cp*prj(l+1)*poi(l))
1668  IF (tem .LT. 0.0) THEN
1669  max1 = 0.0
1670  ELSE
1671  max1 = tem
1672  END IF
1673  bk2(l) = bk2(l+1) + grav*max1/(cp*prj(l+1)*poi(l))
1674  IF (bk2(l) .LT. 0.0) THEN
1675  max2 = 0.0
1676  ELSE
1677  max2 = bk2(l)
1678  END IF
1679  cvw(l) = sqrt(2.0*max2)
1680  END DO
1681 !ALPHA CALCULATION
1682  IF (zet(ic) .LT. 2000.) THEN
1683  rasald(ic) = 0.0_8
1684  rasal(ic) = rasal1
1685  END IF
1686  IF (zet(ic) .GE. 2000.) THEN
1687  rasald(ic) = (rasal2-rasal1)*zetd(ic)/8000.
1688  rasal(ic) = rasal1 + (rasal2-rasal1)*(zet(ic)-2000.)/&
1689 & 8000.
1690  END IF
1691  IF (rasal(ic) .GT. 1.0e5) THEN
1692  rasald(ic) = 0.0_8
1693  rasal(ic) = 1.0e5
1694  ELSE
1695  rasal(ic) = rasal(ic)
1696  END IF
1697  rasald(ic) = -(dt*rasald(ic)/rasal(ic)**2)
1698  rasal(ic) = dt/rasal(ic)
1699  WHERE (cvw(ic:k) .LT. 1.00)
1700  cvw(ic:k) = 1.00
1701  ELSEWHERE
1702  cvw(ic:k) = cvw(ic:k)
1703  END WHERE
1704  CALL acritn(pol(ic), prs(k), acr, acritfac)
1705  IF (wfn .GT. acr) THEN
1706 !================>>
1707  wlq = qol(k)
1708  rnn(k) = 0.
1709  cll0(k) = 0.
1710  DO l=k-1,ic,-1
1711  tem = eta(l) - eta(l+1)
1712  wlq = wlq + tem*qol(l)
1713  IF (l .GT. ic) THEN
1714  tx2 = 0.5*(qst(l)+qst(l-1))*eta(l)
1715  tx3 = 0.5*(hst(l)+hst(l-1))*eta(l)
1716  qcc = tx2 + gm1(l)*(hcc(l)-tx3)
1717  cll0(l) = wlq - qcc
1718  ELSE
1719  cll0(l) = wlq - qst(ic)*eta(ic)
1720  END IF
1721  IF (cll0(l) .LT. 0.00) THEN
1722  cll0(l) = 0.00
1723  ELSE
1724  cll0(l) = cll0(l)
1725  END IF
1726  cli = cll0(l)/eta(l)
1727  te_a = poi(l)*prh(l)
1728  CALL sundq3_ice(te_a, sdqv2, sdqv3, sdqvt1, f2, f3)
1729  c00_x = co_auto(i)*f2*f3*f4
1730  cli_crit_x = cli_crit/(f2*f3)
1731  arg1 = -(cli**2/cli_crit_x**2)
1732  rate = c00_x*(1.0-exp(arg1))
1733  IF (cvw(l) .LT. 1.00) THEN
1734  cvw_x = 1.00
1735  ELSE
1736  cvw_x = cvw(l)
1737  END IF
1738  dt_lyr = (zet(l)-zet(l+1))/cvw_x
1739  closs = cll0(l)*rate*dt_lyr
1740  IF (closs .GT. cll0(l)) THEN
1741  closs = cll0(l)
1742  ELSE
1743  closs = closs
1744  END IF
1745  cll0(l) = cll0(l) - closs
1746  dll0(l) = closs
1747  IF (closs .GT. 0.) THEN
1748  wlq = wlq - closs
1749  rnn(l) = closs
1750  ELSE
1751  rnn(l) = 0.
1752  END IF
1753  END DO
1754  wlq = wlq - qst(ic)*eta(ic)
1755 !CALCULATE GAMMAS AND KERNEL
1756  gmsd(k) = pri(k)*(shtd(k)-ssld(k))
1757  gms(k) = (sht(k)-ssl(k))*pri(k)
1758  gmhd(k) = gmsd(k) + pri(k)*alhl*(qhtd(k)-qold(k))
1759  gmh(k) = gms(k) + (qht(k)-qol(k))*pri(k)*alhl
1760  akmd = dpb(k-1)*(gmhd(k)*gam(k-1)+gmh(k)*gamd(k-1))
1761  akm = gmh(k)*gam(k-1)*dpb(k-1)
1762  tx2d = gmhd(k)
1763  tx2 = gmh(k)
1764  DO l=k-1,ic+1,-1
1765  gmsd(l) = pri(l)*(etad(l)*(sht(l)-ssl(l))+eta(l)*(&
1766 & shtd(l)-ssld(l))+etad(l+1)*(ssl(l)-sht(l+1))+eta(l&
1767 & +1)*(ssld(l)-shtd(l+1)))
1768  gms(l) = (eta(l)*(sht(l)-ssl(l))+eta(l+1)*(ssl(l)-&
1769 & sht(l+1)))*pri(l)
1770  gmhd(l) = gmsd(l) + alhl*pri(l)*(etad(l)*(qht(l)-qol&
1771 & (l))+eta(l)*(qhtd(l)-qold(l))+etad(l+1)*(qol(l)-&
1772 & qht(l+1))+eta(l+1)*(qold(l)-qhtd(l+1)))
1773  gmh(l) = gms(l) + (eta(l)*(qht(l)-qol(l))+eta(l+1)*(&
1774 & qol(l)-qht(l+1)))*alhl*pri(l)
1775  tx2d = tx2d + (etad(l)-etad(l+1))*gmh(l) + (eta(l)-&
1776 & eta(l+1))*gmhd(l)
1777  tx2 = tx2 + (eta(l)-eta(l+1))*gmh(l)
1778  akmd = akmd - pki(l)*(gmsd(l)*eht(l)+gms(l)*ehtd(l))&
1779 & + tx2d*ght(l) + tx2*ghtd(l)
1780  akm = akm - gms(l)*eht(l)*pki(l) + tx2*ght(l)
1781  END DO
1782  gmsd(ic) = pri(ic)*(etad(ic+1)*(ssl(ic)-sht(ic+1))+eta&
1783 & (ic+1)*(ssld(ic)-shtd(ic+1)))
1784  gms(ic) = eta(ic+1)*(ssl(ic)-sht(ic+1))*pri(ic)
1785  akmd = akmd - dpb(ic)*pki(ic)*(gmsd(ic)*eta(ic+1)+gms(&
1786 & ic)*etad(ic+1))
1787  akm = akm - gms(ic)*eta(ic+1)*dpb(ic)*pki(ic)
1788  gmhd(ic) = gmsd(ic) + pri(ic)*(alhl*(etad(ic+1)*(qol(&
1789 & ic)-qht(ic+1))+eta(ic+1)*(qold(ic)-qhtd(ic+1)))+etad&
1790 & (ic)*(hst(ic)-hol(ic))+eta(ic)*(hstd(ic)-hold(ic)))
1791  gmh(ic) = gms(ic) + (eta(ic+1)*(qol(ic)-qht(ic+1))*&
1792 & alhl+eta(ic)*(hst(ic)-hol(ic)))*pri(ic)
1793 !CLOUD BASE MASS FLUX
1794  IF (.NOT.(akm .GE. 0.0 .OR. wlq .LT. 0.0)) THEN
1795 !================>>
1796  wfnd = -((wfnd*akm-(wfn-acr)*akmd)/akm**2)
1797  wfn = -((wfn-acr)/akm)
1798  x1d = (rasald(ic)*wfn+rasal(ic)*wfnd)*trg*toki + &
1799 & rasal(ic)*wfn*(trgd*toki+trg*tokid)
1800  x1 = rasal(ic)*trg*toki*wfn
1801  IF (x1 .GT. (prs(k+1)-prs(k))*(100.*pblfrac)) THEN
1802  wfn = (prs(k+1)-prs(k))*(100.*pblfrac)
1803  wfnd = 0.0_8
1804  ELSE
1805  wfnd = x1d
1806  wfn = x1
1807  END IF
1808 !CUMULATIVE PRECIP AND CLOUD-BASE MASS FLUX FOR OUTPUT
1809  tem = wfn*gravi
1810  cll(ic) = cll(ic) + wlq*tem
1811  rmf(ic) = rmf(ic) + tem
1812  rmfd(ic) = rmfd(ic) + tem*eta(ic)
1813  DO l=ic+1,k
1814  rmfp(l) = tem*eta(l)
1815  rmfc(l) = rmfc(l) + rmfp(l)
1816  dllx(l) = dllx(l) + tem*dll0(l)
1817  IF (cvw(l) .GT. 0.0) THEN
1818  updfrp(l) = rmfp(l)*(ddt/daylen)*1000./(cvw(l)*&
1819 & prs(l))
1820  ELSE
1821  updfrp(l) = 0.0
1822  END IF
1823  clli(l) = cll0(l)/eta(l)
1824  updfrc(l) = updfrc(l) + updfrp(l)
1825  END DO
1826 !THETA AND Q CHANGE DUE TO CLOUD TYPE IC
1827  DO l=ic,k
1828  rns(l) = rns(l) + rnn(l)*tem
1829  gmhd(l) = gmhd(l)*wfn + gmh(l)*wfnd
1830  gmh(l) = gmh(l)*wfn
1831  gmsd(l) = gmsd(l)*wfn + gms(l)*wfnd
1832  gms(l) = gms(l)*wfn
1833  qoid(l) = qoid(l) + alhi*(gmhd(l)-gmsd(l))
1834  qoi(l) = qoi(l) + (gmh(l)-gms(l))*alhi
1835  poid(l) = poid(l) + pki(l)*cpi*gmsd(l)
1836  poi(l) = poi(l) + gms(l)*pki(l)*cpi
1837  qstd(l) = qstd(l) + cpi*(gmsd(l)*bet(l)+gms(l)*&
1838 & betd(l))
1839  qst(l) = qst(l) + gms(l)*bet(l)*cpi
1840  END DO
1841  END IF
1842  END IF
1843  END IF
1844  END IF
1845  END IF
1846  END IF
1847  END DO
1848 !CLOUD LOOP
1849  IF (sum(rmf(icmin:k)) .GT. 0.0) THEN
1850  thod(i, icmin:k-1) = poid(icmin:k-1)
1851  tho(i, icmin:k-1) = poi(icmin:k-1)
1852  qhod(i, icmin:k-1) = qoid(icmin:k-1)
1853  qho(i, icmin:k-1) = qoi(icmin:k-1)
1854 !De-strap tendencies from RAS
1855  wght = wgt1(i, :)
1856 !Scale properly by layer masses
1857  wght0 = 0.
1858  DO l=k,k0
1859  wght0 = wght0 + wght(l)*(ple(i, l+1)-ple(i, l))
1860  END DO
1861  wght0 = (prs(k+1)-prs(k))/wght0
1862  wght = wght0*wght
1863  DO l=k,k0
1864  thod(i, l) = thod(i, l) + wght(l)*(poid(k)-poi_svd(k))
1865  tho(i, l) = tho(i, l) + wght(l)*(poi(k)-poi_sv(k))
1866  qhod(i, l) = qhod(i, l) + wght(l)*(qoid(k)-qoi_svd(k))
1867  qho(i, l) = qho(i, l) + wght(l)*(qoi(k)-qoi_sv(k))
1868  END DO
1869  END IF
1870  END IF
1871  END DO
1872 END SUBROUTINE rase0_d
1873 
1874 SUBROUTINE rase_tracer_d(idim, irun, k0, icmin, dt, cons_cp, cons_alhl, &
1875 & cons_grav, cons_rgas, cons_h2omw, cons_airmw, cons_vireps, seedras, &
1876 & sige, kcbl, wgt0, wgt1, frland, ts, thoin, qhoin, uhoin, vhoin, &
1877 & co_auto, ple, rasparams, estblx, itrcr, xho, xhod, fscav)
1878  IMPLICIT NONE
1879 !INPUTS
1880  INTEGER, INTENT(IN) :: idim, irun, k0, icmin
1881  REAL*8, DIMENSION(idim, k0 + 1), INTENT(IN) :: ple
1882  REAL*8, DIMENSION(k0 + 1), INTENT(IN) :: sige
1883  REAL*8, INTENT(IN) :: dt, cons_cp, cons_alhl, cons_grav, cons_rgas
1884  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_vireps
1885  INTEGER, DIMENSION(idim), INTENT(IN) :: seedras
1886  INTEGER, DIMENSION(idim), INTENT(IN) :: kcbl
1887  REAL*8, DIMENSION(idim), INTENT(IN) :: ts, frland
1888  REAL*8, DIMENSION(idim), INTENT(IN) :: co_auto
1889  REAL*8, DIMENSION(idim, k0), INTENT(IN) :: wgt0, wgt1
1890  REAL*8, DIMENSION(:), INTENT(IN) :: rasparams
1891  REAL*8, DIMENSION(:), INTENT(IN) :: estblx
1892  INTEGER, INTENT(IN) :: itrcr
1893  REAL*8, DIMENSION(itrcr), INTENT(IN) :: fscav
1894  REAL*8, DIMENSION(idim, k0), INTENT(IN) :: thoin, qhoin, uhoin, vhoin
1895 !PROGNOSTIC
1896  REAL*8, DIMENSION(idim, k0, itrcr), INTENT(INOUT) :: xho
1897  REAL*8, DIMENSION(idim, k0, itrcr), INTENT(INOUT) :: xhod
1898 !LOCALS
1899  INTEGER :: i, ic, l, kk, k
1900 !Parameters
1901  REAL*8, PARAMETER :: onepkap=1.+2./7., daylen=86400.0
1902  REAL*8, PARAMETER :: rhmax=0.9999
1903  REAL*8, PARAMETER :: cbl_qpert=0.0, cbl_tpert=1.0
1904  REAL*8, PARAMETER :: cbl_tpert_mxocn=2.0, cbl_tpert_mxlnd=4.0
1905 !Constants
1906  REAL*8 :: grav, cp, alhl, cpbg, alhi, cpi, gravi, ddt, lbcp
1907 !Rasparams
1908  REAL*8 :: fricfac, cli_crit, rasal1, rasal2
1909  REAL*8 :: friclambda
1910  REAL*8 :: sdqv2, sdqv3, sdqvt1
1911  REAL*8 :: acritfac, pblfrac, autorampb
1912  REAL*8 :: maxdallowed, rhmn, rhmx
1913  REAL*8 :: mxdiam
1914  REAL*8 :: tx2, tx3, akm, acr, alm, tth, qqh, dqx
1915  REAL*8 :: wfn, tem, trg, trgexp, evp, wlq, qcc
1916  REAL*8 :: cli, te_a, c00_x, cli_crit_x, toki
1917  REAL*8 :: dt_lyr, rate, cvw_x, closs, f2, f3, f4
1918  REAL*8 :: wght0, prcbl, rndu
1919  REAL*8 :: lambda_min, lambda_max
1920  REAL*8 :: tpert, qpert
1921  REAL*8 :: uht, vht
1922  REAL*8, DIMENSION(k0) :: poi_sv, qoi_sv, uoi_sv, voi_sv
1923  REAL*8, DIMENSION(k0) :: poi, qoi, uoi, voi, dqq, bet, gam, cll
1924  REAL*8, DIMENSION(k0) :: poid, qoid, dqqd, betd, gamd
1925  REAL*8, DIMENSION(k0) :: poi_c, qoi_c
1926  REAL*8, DIMENSION(k0) :: poi_cd, qoi_cd
1927  REAL*8, DIMENSION(k0) :: prh, pri, ght, dpt, dpb, pki
1928  REAL*8, DIMENSION(k0) :: prhd, prid, ghtd, dptd, dpbd, pkid
1929  REAL*8, DIMENSION(k0) :: ucu, vcu
1930  REAL*8, DIMENSION(k0) :: cln, rns, pol
1931  REAL*8, DIMENSION(k0) :: pold
1932  REAL*8, DIMENSION(k0) :: qst, ssl, rmf, rnn, rn1, rmfc, rmfp
1933  REAL*8, DIMENSION(k0) :: qstd, ssld
1934  REAL*8, DIMENSION(k0) :: gms, eta, gmh, eht, gm1, hcc, rmfd
1935  REAL*8, DIMENSION(k0) :: gmsd, etad, gmhd, ehtd, hccd
1936  REAL*8, DIMENSION(k0) :: hol, hst, qol, zol, hcld, cll0, cllx, clli
1937  REAL*8, DIMENSION(k0) :: hold, hstd, qold, zold
1938  REAL*8, DIMENSION(k0) :: bke, cvw, updfrc
1939  REAL*8, DIMENSION(k0) :: rasal, updfrp, bk2, dll0, dllx
1940  REAL*8, DIMENSION(k0) :: rasald
1941  REAL*8, DIMENSION(k0) :: wght, massf
1942  REAL*8, DIMENSION(k0) :: wghtd
1943  REAL*8, DIMENSION(k0) :: qss, dqs, pf, pk, tempf, zlo
1944  REAL*8, DIMENSION(k0) :: zlod
1945  REAL*8, DIMENSION(k0 + 1) :: prj, prs, qht, sht, zet, zle, pke
1946  REAL*8, DIMENSION(k0+1) :: prjd, prsd, qhtd, shtd, zetd, zled
1947  REAL*8, DIMENSION(idim, k0) :: tho, qho, uho, vho
1948 !Tracer scavenging
1949  INTEGER :: itr
1950 !Layer thickness in km
1951  REAL*8 :: delzkm
1952 !Fraction of tracer *not* scavenged
1953  REAL*8 :: fnoscav
1954  REAL*8, DIMENSION(k0, itrcr) :: xoi, xcu, xoi_sv
1955  REAL*8, DIMENSION(k0, itrcr) :: xoid, xcud, xoi_svd
1956  REAL*8, DIMENSION(itrcr) :: xht
1957  REAL*8, DIMENSION(itrcr) :: xhtd
1958  INTRINSIC max
1959  INTRINSIC min
1960  INTRINSIC sqrt
1961  INTRINSIC exp
1962  INTRINSIC sum
1963  REAL*8, DIMENSION(k0+1) :: pwx1
1964  REAL*8 :: pwy1
1965  REAL*8, DIMENSION(k0) :: pwx10
1966  REAL*8 :: pwx11
1967  REAL*8 :: pwr1
1968  REAL*8 :: arg1
1969  REAL*8 :: x5
1970  REAL*8 :: x4
1971  REAL*8 :: x3
1972  REAL*8 :: x2
1973  REAL*8 :: x1
1974  REAL*8 :: max2
1975  REAL*8 :: max1
1976  REAL*8 :: y1
1977 !Pass meteorology to internal arrays so it is not updated
1978  tho = thoin
1979  qho = qhoin
1980  uho = uhoin
1981  vho = vhoin
1982 !Initialize Local Arrays
1983  poi = 0.0
1984  qoi = 0.0
1985  uoi = 0.0
1986  voi = 0.0
1987  dqq = 0.0
1988  bet = 0.0
1989  gam = 0.0
1990  poi_c = 0.0
1991  qoi_c = 0.0
1992  prh = 0.0
1993  pri = 0.0
1994  ght = 0.0
1995  dpt = 0.0
1996  dpb = 0.0
1997  pki = 0.0
1998  ucu = 0.0
1999  vcu = 0.0
2000  cln = 0.0
2001  pol = 0.0
2002  qst = 0.0
2003  ssl = 0.0
2004  rmf = 0.0
2005  rnn = 0.0
2006  rn1 = 0.0
2007  gms = 0.0
2008  eta = 0.0
2009  gmh = 0.0
2010  eht = 0.0
2011  gm1 = 0.0
2012  hcc = 0.0
2013  hol = 0.0
2014  hst = 0.0
2015  qol = 0.0
2016  zol = 0.0
2017  hcld = 0.0
2018  bke = 0.0
2019  cvw = 0.0
2020  updfrc = 0.0
2021  rasal = 0.0
2022  updfrp = 0.0
2023  bk2 = 0.0
2024  wght = 0.0
2025  massf = 0.0
2026  qss = 0.0
2027  dqs = 0.0
2028  pf = 0.0
2029  pk = 0.0
2030  tempf = 0.0
2031  zlo = 0.0
2032  prj = 0.0
2033  prs = 0.0
2034  qht = 0.0
2035  sht = 0.0
2036  zet = 0.0
2037  zle = 0.0
2038  pke = 0.0
2039 ! --- 1
2040  fricfac = rasparams(1)
2041 ! --- 4
2042  cli_crit = rasparams(4)
2043 ! --- 5
2044  rasal1 = rasparams(5)
2045 ! --- 6
2046  rasal2 = rasparams(6)
2047 ! --- 11
2048  friclambda = rasparams(11)
2049 ! --- 14
2050  sdqv2 = rasparams(14)
2051 ! --- 15
2052  sdqv3 = rasparams(15)
2053 ! --- 16
2054  sdqvt1 = rasparams(16)
2055 ! --- 17
2056  acritfac = rasparams(17)
2057 ! --- 20
2058  pblfrac = rasparams(20)
2059 ! --- 21
2060  autorampb = rasparams(21)
2061 ! --- 24
2062  rhmn = rasparams(24)
2063 ! --- 24
2064  maxdallowed = rasparams(23)
2065 ! --- 25
2066  rhmx = rasparams(25)
2067  grav = cons_grav
2068  alhl = cons_alhl
2069  cp = cons_cp
2070  cpi = 1.0/cp
2071  alhi = 1.0/alhl
2072  gravi = 1.0/grav
2073  cpbg = cp*gravi
2074  ddt = daylen/dt
2075  lbcp = alhl*cpi
2076  i = 1
2077 !CALL FINDBASE
2078  k = kcbl(i)
2079  IF (k .GT. 0) THEN
2080 !Get saturation specific humidity and gradient wrt to T
2081  pwx1 = ple(i, :)/1000.
2082  pwy1 = cons_rgas/cons_cp
2083  pke = pwx1**pwy1
2084  pf = 0.5*(ple(i, 1:k0)+ple(i, 2:k0+1))
2085  pwx10 = pf/1000.
2086  pwy1 = cons_rgas/cons_cp
2087  pk = pwx10**pwy1
2088  tempf = tho(i, :)*pk
2089  zle = 0.0
2090  zlo = 0.0
2091  zled(k0+1) = 0.0_8
2092  zle(k0+1) = 0.
2093  DO l=k0,1,-1
2094  zled(l) = 0.0_8
2095  zle(l) = tho(i, l)*(1.+cons_vireps*qho(i, l))
2096  zlod(l) = 0.0_8
2097  zlo(l) = zle(l+1) + cons_cp/cons_grav*(pke(l+1)-pk(l))*zle(l)
2098  zled(l) = 0.0_8
2099  zle(l) = zlo(l) + cons_cp/cons_grav*(pk(l)-pke(l))*zle(l)
2100  END DO
2101  tpert = cbl_tpert*(ts(i)-(tempf(k0)+cons_grav*zlo(k0)/cons_cp))
2102 !* ( QSSFC - Q(:,:,K0) ) [CBL_QPERT = 0.0]
2103  qpert = cbl_qpert
2104  IF (tpert .LT. 0.0) THEN
2105  tpert = 0.0
2106  ELSE
2107  tpert = tpert
2108  END IF
2109  IF (qpert .LT. 0.0) THEN
2110  qpert = 0.0
2111  ELSE
2112  qpert = qpert
2113  END IF
2114  IF (frland(i) .LT. 0.1) THEN
2115  IF (tpert .GT. cbl_tpert_mxocn) THEN
2116  tpert = cbl_tpert_mxocn
2117  ELSE
2118  tpert = tpert
2119  END IF
2120  ELSE IF (tpert .GT. cbl_tpert_mxlnd) THEN
2121  tpert = cbl_tpert_mxlnd
2122  ELSE
2123  tpert = tpert
2124  END IF
2125  CALL dqsat_ras(dqs, qss, tempf, pf, k0, estblx, cons_h2omw, &
2126 & cons_airmw)
2127  DO kk=icmin,k+1
2128  prjd(kk) = 0.0_8
2129  prj(kk) = pke(kk)
2130  END DO
2131  prsd(icmin:k0+1) = 0.0_8
2132  prs(icmin:k0+1) = ple(i, icmin:k0+1)
2133  poid(icmin:k) = 0.0_8
2134  poi(icmin:k) = tho(i, icmin:k)
2135  qoid(icmin:k) = 0.0_8
2136  qoi(icmin:k) = qho(i, icmin:k)
2137  uoi(icmin:k) = uho(i, icmin:k)
2138  voi(icmin:k) = vho(i, icmin:k)
2139  qstd(icmin:k) = 0.0_8
2140  qst(icmin:k) = qss(icmin:k)
2141  dqqd(icmin:k) = 0.0_8
2142  dqq(icmin:k) = dqs(icmin:k)
2143  xoid = 0.0_8
2144 !DO_TRACERS
2145  DO itr=1,itrcr
2146  xoid(icmin:k, itr) = xhod(i, icmin:k, itr)
2147  xoi(icmin:k, itr) = xho(i, icmin:k, itr)
2148  END DO
2149 !Mass fraction of each layer below cloud base
2150  massf(:) = wgt0(i, :)
2151 !RESET PRESSURE at bottom edge of CBL
2152  prcbl = prs(k)
2153  DO l=k,k0
2154  prcbl = prcbl + massf(l)*(prs(l+1)-prs(l))
2155  END DO
2156  prsd(k+1) = 0.0_8
2157  prs(k+1) = prcbl
2158  pwx11 = prs(k+1)/1000.
2159  pwy1 = cons_rgas/cons_cp
2160  prjd(k+1) = 0.0_8
2161  prj(k+1) = pwx11**pwy1
2162  DO l=k,icmin,-1
2163  pold(l) = 0.0_8
2164  pol(l) = 0.5*(prs(l)+prs(l+1))
2165  prhd(l) = 0.0_8
2166  prh(l) = (prs(l+1)*prj(l+1)-prs(l)*prj(l))/(onepkap*(prs(l+1)-prs(&
2167 & l)))
2168  pkid(l) = 0.0_8
2169  pki(l) = 1.0/prh(l)
2170  dptd(l) = 0.0_8
2171  dpt(l) = prh(l) - prj(l)
2172  dpbd(l) = 0.0_8
2173  dpb(l) = prj(l+1) - prh(l)
2174  prid(l) = 0.0_8
2175  pri(l) = .01/(prs(l+1)-prs(l))
2176  END DO
2177 !RECALCULATE PROFILE QUAN. IN LOWEST STRAPPED LAYER
2178  IF (k .LE. k0) THEN
2179  poid(k) = 0.0_8
2180  poi(k) = 0.
2181  qoid(k) = 0.0_8
2182  qoi(k) = 0.
2183  uoi(k) = 0.
2184  voi(k) = 0.
2185 !SPECIFY WEIGHTS GIVEN TO EACH LAYER WITHIN SUBCLOUD "SUPERLAYER"
2186  wght = 0.
2187  DO l=k,k0
2188  wghtd(l) = 0.0_8
2189  wght(l) = massf(l)*(ple(i, l+1)-ple(i, l))/(prs(k+1)-prs(k))
2190  END DO
2191  DO l=k,k0
2192  poid(k) = 0.0_8
2193  poi(k) = poi(k) + wght(l)*tho(i, l)
2194  qoid(k) = 0.0_8
2195  qoi(k) = qoi(k) + wght(l)*qho(i, l)
2196  uoi(k) = uoi(k) + wght(l)*uho(i, l)
2197  voi(k) = voi(k) + wght(l)*vho(i, l)
2198  END DO
2199 !DO_TRACERS
2200  xoid(k, :) = 0.0_8
2201  xoi(k, :) = 0.
2202  DO itr=1,itrcr
2203  DO l=k,k0
2204  xoid(k, itr) = xoid(k, itr) + wght(l)*xhod(i, l, itr)
2205  xoi(k, itr) = xoi(k, itr) + wght(l)*xho(i, l, itr)
2206  END DO
2207  END DO
2208  CALL dqsats_ras(dqq(k), qst(k), poi(k)*prh(k), pol(k), estblx, &
2209 & cons_h2omw, cons_airmw)
2210  END IF
2211  IF (seedras(i)/1000000. .LT. 1e-6) THEN
2212  rndu = 1e-6
2213  ELSE
2214  rndu = seedras(i)/1000000.
2215  END IF
2216  pwr1 = rndu**(-(1./2.))
2217  mxdiam = maxdallowed*pwr1
2218  DO l=k,icmin,-1
2219 !*
2220  betd(l) = 0.0_8
2221  bet(l) = dqq(l)*pki(l)
2222 !*
2223  gamd(l) = 0.0_8
2224  gam(l) = pki(l)/(1.0+lbcp*dqq(l))
2225  IF (l .LT. k) THEN
2226  ghtd(l+1) = 0.0_8
2227  ght(l+1) = gam(l)*dpb(l) + gam(l+1)*dpt(l+1)
2228  gm1(l+1) = 0.5*lbcp*(dqq(l)/(alhl*(1.0+lbcp*dqq(l)))+dqq(l+1)/(&
2229 & alhl*(1.0+lbcp*dqq(l+1))))
2230  END IF
2231  END DO
2232  rns = 0.
2233  cll = 0.
2234  rmf = 0.
2235  rmfd = 0.
2236  rmfc = 0.
2237  rmfp = 0.
2238  cll0 = 0.
2239  dll0 = 0.
2240  cllx = 0.
2241  dllx = 0.
2242  clli = 0.
2243  poi_sv = poi
2244  qoi_sv = qoi
2245  uoi_sv = uoi
2246  voi_sv = voi
2247 !DO_TRACERS
2248  xoi_svd = xoid
2249  xoi_sv = xoi
2250  cvw = 0.0
2251  updfrc = 0.0
2252  updfrp = 0.0
2253 ! HOL initialized here in order not to confuse Valgrind debugger
2254  hol = 0.
2255  zetd(k+1) = 0.0_8
2256  zet(k+1) = 0
2257  shtd(k+1) = 0.0_8
2258  sht(k+1) = cp*poi(k)*prj(k+1)
2259  DO l=k,icmin,-1
2260  IF (qst(l)*rhmax .GT. qoi(l)) THEN
2261  qold(l) = 0.0_8
2262  qol(l) = qoi(l)
2263  ELSE
2264  qold(l) = 0.0_8
2265  qol(l) = qst(l)*rhmax
2266  END IF
2267  IF (0.000 .LT. qol(l)) THEN
2268  qold(l) = 0.0_8
2269  qol(l) = qol(l)
2270  ELSE
2271  qold(l) = 0.0_8
2272  qol(l) = 0.000
2273  END IF
2274  ssld(l) = 0.0_8
2275  ssl(l) = cp*prj(l+1)*poi(l) + grav*zet(l+1)
2276  hold(l) = 0.0_8
2277  hol(l) = ssl(l) + qol(l)*alhl
2278  hstd(l) = 0.0_8
2279  hst(l) = ssl(l) + qst(l)*alhl
2280  tem = poi(l)*(prj(l+1)-prj(l))*cpbg
2281  zetd(l) = 0.0_8
2282  zet(l) = zet(l+1) + tem
2283  zold(l) = 0.0_8
2284  zol(l) = zet(l+1) + (prj(l+1)-prh(l))*poi(l)*cpbg
2285  END DO
2286  xcud = 0.0_8
2287  xhtd = 0.0_8
2288  DO ic=k,icmin+1,-1
2289 !DO_TRACERS
2290  xcud(icmin:, :) = 0.0_8
2291  xcu(icmin:, :) = 0.
2292  ucu(icmin:) = 0.
2293  vcu(icmin:) = 0.
2294  alm = 0.
2295  IF (1. .GT. (qoi(k)/qst(k)-rhmn)/(rhmx-rhmn)) THEN
2296  trg = (qoi(k)/qst(k)-rhmn)/(rhmx-rhmn)
2297  ELSE
2298  trg = 1.
2299  END IF
2300  IF (0.0 .LT. (autorampb-sige(ic))/0.2) THEN
2301  y1 = (autorampb-sige(ic))/0.2
2302  ELSE
2303  y1 = 0.0
2304  END IF
2305  IF (1.0 .GT. y1) THEN
2306  f4 = y1
2307  ELSE
2308  f4 = 1.0
2309  END IF
2310  IF (trg .GT. 1.0e-5) THEN
2311 !================>>
2312 !RECOMPUTE SOUNDING UP TO DETRAINMENT LEVEL
2313  poi_c = poi
2314  qoi_c = qoi
2315  poi_cd(k) = 0.0_8
2316  poi_c(k) = poi_c(k) + tpert
2317  qoi_cd(k) = 0.0_8
2318  qoi_c(k) = qoi_c(k) + qpert
2319  zetd(k+1) = 0.0_8
2320  zet(k+1) = 0.
2321  shtd(k+1) = 0.0_8
2322  sht(k+1) = cp*poi_c(k)*prj(k+1)
2323  DO l=k,ic,-1
2324  IF (qst(l)*rhmax .GT. qoi_c(l)) THEN
2325  qold(l) = 0.0_8
2326  qol(l) = qoi_c(l)
2327  ELSE
2328  qold(l) = 0.0_8
2329  qol(l) = qst(l)*rhmax
2330  END IF
2331  IF (0.000 .LT. qol(l)) THEN
2332  qold(l) = 0.0_8
2333  qol(l) = qol(l)
2334  ELSE
2335  qold(l) = 0.0_8
2336  qol(l) = 0.000
2337  END IF
2338  ssld(l) = 0.0_8
2339  ssl(l) = cp*prj(l+1)*poi_c(l) + grav*zet(l+1)
2340  hold(l) = 0.0_8
2341  hol(l) = ssl(l) + qol(l)*alhl
2342  hstd(l) = 0.0_8
2343  hst(l) = ssl(l) + qst(l)*alhl
2344  tem = poi_c(l)*(prj(l+1)-prj(l))*cpbg
2345  zetd(l) = 0.0_8
2346  zet(l) = zet(l+1) + tem
2347  zold(l) = 0.0_8
2348  zol(l) = zet(l+1) + (prj(l+1)-prh(l))*poi_c(l)*cpbg
2349  END DO
2350  DO l=ic+1,k
2351  tem = (prj(l)-prh(l-1))/(prh(l)-prh(l-1))
2352  shtd(l) = 0.0_8
2353  sht(l) = ssl(l-1) + tem*(ssl(l)-ssl(l-1))
2354  qhtd(l) = 0.0_8
2355  qht(l) = .5*(qol(l)+qol(l-1))
2356  END DO
2357 !CALCULATE LAMBDA, ETA, AND WORKFUNCTION
2358  lambda_min = .2/mxdiam
2359  lambda_max = .2/200.
2360  IF (hol(k) .GT. hst(ic)) THEN
2361 !================>>
2362 !LAMBDA CALCULATION: MS-A18
2363  tem = (hst(ic)-hol(ic))*(zol(ic)-zet(ic+1))
2364  DO l=ic+1,k-1
2365  tem = tem + (hst(ic)-hol(l))*(zet(l)-zet(l+1))
2366  END DO
2367  IF (tem .GT. 0.0) THEN
2368 !================>>
2369  alm = (hol(k)-hst(ic))/tem
2370  IF (alm .LE. lambda_max) THEN
2371 !================>>
2372  toki = 1.0
2373  IF (alm .LT. lambda_min) toki = (alm/lambda_min)**2
2374 !ETA CALCULATION: MS-A2
2375  DO l=ic+1,k
2376  etad(l) = 0.0_8
2377  eta(l) = 1.0 + alm*(zet(l)-zet(k))
2378  END DO
2379  etad(ic) = 0.0_8
2380  eta(ic) = 1.0 + alm*(zol(ic)-zet(k))
2381 !WORKFUNCTION CALCULATION: MS-A22
2382  wfn = 0.0
2383  hccd(k) = 0.0_8
2384  hcc(k) = hol(k)
2385  DO l=k-1,ic+1,-1
2386  hccd(l) = 0.0_8
2387  hcc(l) = hcc(l+1) + (eta(l)-eta(l+1))*hol(l)
2388  tem = hcc(l+1)*dpb(l) + hcc(l)*dpt(l)
2389  ehtd(l) = 0.0_8
2390  eht(l) = eta(l+1)*dpb(l) + eta(l)*dpt(l)
2391  wfn = wfn + (tem-eht(l)*hst(l))*gam(l)
2392  END DO
2393  hccd(ic) = 0.0_8
2394  hcc(ic) = hst(ic)*eta(ic)
2395  wfn = wfn + (hcc(ic+1)-hst(ic)*eta(ic+1))*gam(ic)*dpb(ic)
2396 !VERTICAL VELOCITY/KE CALCULATION (ADDED 12/2001 JTB)
2397  bk2(k) = 0.0
2398  bke(k) = 0.0
2399  hcld(k) = hol(k)
2400  DO l=k-1,ic,-1
2401  hcld(l) = (eta(l+1)*hcld(l+1)+(eta(l)-eta(l+1))*hol(l))/&
2402 & eta(l)
2403  tem = (hcld(l)-hst(l))*(zet(l)-zet(l+1))/(1.0+lbcp*dqq(l&
2404 & ))
2405  bke(l) = bke(l+1) + grav*tem/(cp*prj(l+1)*poi(l))
2406  IF (tem .LT. 0.0) THEN
2407  max1 = 0.0
2408  ELSE
2409  max1 = tem
2410  END IF
2411  bk2(l) = bk2(l+1) + grav*max1/(cp*prj(l+1)*poi(l))
2412  IF (bk2(l) .LT. 0.0) THEN
2413  max2 = 0.0
2414  ELSE
2415  max2 = bk2(l)
2416  END IF
2417  cvw(l) = sqrt(2.0*max2)
2418  END DO
2419 !ALPHA CALCULATION
2420  IF (zet(ic) .LT. 2000.) THEN
2421  rasald(ic) = 0.0_8
2422  rasal(ic) = rasal1
2423  END IF
2424  IF (zet(ic) .GE. 2000.) THEN
2425  rasald(ic) = 0.0_8
2426  rasal(ic) = rasal1 + (rasal2-rasal1)*(zet(ic)-2000.)/&
2427 & 8000.
2428  END IF
2429  IF (rasal(ic) .GT. 1.0e5) THEN
2430  rasald(ic) = 0.0_8
2431  rasal(ic) = 1.0e5
2432  ELSE
2433  rasald(ic) = 0.0_8
2434  rasal(ic) = rasal(ic)
2435  END IF
2436  rasald(ic) = 0.0_8
2437  rasal(ic) = dt/rasal(ic)
2438  DO l=ic,k
2439  IF (cvw(l) .LT. 1.00) THEN
2440  cvw(l) = 1.00
2441  ELSE
2442  cvw(l) = cvw(l)
2443  END IF
2444  END DO
2445  CALL acritn(pol(ic), prs(k), acr, acritfac)
2446  IF (wfn .GT. acr) THEN
2447 !================>>
2448 !DO_TRACERS
2449  DO itr=1,itrcr
2450 !Scavenging of the below cloud tracer
2451  delzkm = (zet(ic)-zet(k))/1000.
2452  x4 = exp(-(fscav(itr)*delzkm))
2453  IF (x4 .GT. 1.) THEN
2454  x1 = 1.
2455  ELSE
2456  x1 = x4
2457  END IF
2458  IF (x1 .LT. 0.) THEN
2459  fnoscav = 0.
2460  ELSE
2461  fnoscav = x1
2462  END IF
2463  xhtd(itr) = fnoscav*xoid(k, itr)
2464  xht(itr) = xoi(k, itr)*fnoscav
2465  END DO
2466  wlq = qol(k)
2467  uht = uoi(k)
2468  vht = voi(k)
2469  rnn(k) = 0.
2470  cll0(k) = 0.
2471  DO l=k-1,ic,-1
2472  tem = eta(l) - eta(l+1)
2473  wlq = wlq + tem*qol(l)
2474  uht = uht + tem*uoi(l)
2475  vht = vht + tem*voi(l)
2476 !DO_TRACERS
2477  DO itr=1,itrcr
2478 !Scavenging of the entrained tracer. Updates transported tracer mass.
2479  delzkm = (zet(ic)-zet(l+1))/1000.
2480  x5 = exp(-(fscav(itr)*delzkm))
2481  IF (x5 .GT. 1.) THEN
2482  x2 = 1.
2483  ELSE
2484  x2 = x5
2485  END IF
2486  IF (x2 .LT. 0.) THEN
2487  fnoscav = 0.
2488  ELSE
2489  fnoscav = x2
2490  END IF
2491  xhtd(itr) = xhtd(itr) + tem*fnoscav*xoid(l, itr)
2492  xht(itr) = xht(itr) + tem*xoi(l, itr)*fnoscav
2493  END DO
2494  IF (l .GT. ic) THEN
2495  tx2 = 0.5*(qst(l)+qst(l-1))*eta(l)
2496  tx3 = 0.5*(hst(l)+hst(l-1))*eta(l)
2497  qcc = tx2 + gm1(l)*(hcc(l)-tx3)
2498  cll0(l) = wlq - qcc
2499  ELSE
2500  cll0(l) = wlq - qst(ic)*eta(ic)
2501  END IF
2502  IF (cll0(l) .LT. 0.00) THEN
2503  cll0(l) = 0.00
2504  ELSE
2505  cll0(l) = cll0(l)
2506  END IF
2507  cli = cll0(l)/eta(l)
2508  te_a = poi(l)*prh(l)
2509  CALL sundq3_ice(te_a, sdqv2, sdqv3, sdqvt1, f2, f3)
2510  c00_x = co_auto(i)*f2*f3*f4
2511  cli_crit_x = cli_crit/(f2*f3)
2512  arg1 = -(cli**2/cli_crit_x**2)
2513  rate = c00_x*(1.0-exp(arg1))
2514  IF (cvw(l) .LT. 1.00) THEN
2515  cvw_x = 1.00
2516  ELSE
2517  cvw_x = cvw(l)
2518  END IF
2519  dt_lyr = (zet(l)-zet(l+1))/cvw_x
2520  closs = cll0(l)*rate*dt_lyr
2521  IF (closs .GT. cll0(l)) THEN
2522  closs = cll0(l)
2523  ELSE
2524  closs = closs
2525  END IF
2526  cll0(l) = cll0(l) - closs
2527  dll0(l) = closs
2528  IF (closs .GT. 0.) THEN
2529  wlq = wlq - closs
2530  rnn(l) = closs
2531  ELSE
2532  rnn(l) = 0.
2533  END IF
2534  END DO
2535  wlq = wlq - qst(ic)*eta(ic)
2536 !CALCULATE GAMMAS AND KERNEL
2537  gmsd(k) = 0.0_8
2538  gms(k) = (sht(k)-ssl(k))*pri(k)
2539  gmhd(k) = 0.0_8
2540  gmh(k) = gms(k) + (qht(k)-qol(k))*pri(k)*alhl
2541  akm = gmh(k)*gam(k-1)*dpb(k-1)
2542  tx2 = gmh(k)
2543  DO l=k-1,ic+1,-1
2544  gmsd(l) = 0.0_8
2545  gms(l) = (eta(l)*(sht(l)-ssl(l))+eta(l+1)*(ssl(l)-sht(&
2546 & l+1)))*pri(l)
2547  gmhd(l) = 0.0_8
2548  gmh(l) = gms(l) + (eta(l)*(qht(l)-qol(l))+eta(l+1)*(&
2549 & qol(l)-qht(l+1)))*alhl*pri(l)
2550  tx2 = tx2 + (eta(l)-eta(l+1))*gmh(l)
2551  akm = akm - gms(l)*eht(l)*pki(l) + tx2*ght(l)
2552  END DO
2553  gmsd(ic) = 0.0_8
2554  gms(ic) = eta(ic+1)*(ssl(ic)-sht(ic+1))*pri(ic)
2555  akm = akm - gms(ic)*eta(ic+1)*dpb(ic)*pki(ic)
2556  gmhd(ic) = 0.0_8
2557  gmh(ic) = gms(ic) + (eta(ic+1)*(qol(ic)-qht(ic+1))*alhl+&
2558 & eta(ic)*(hst(ic)-hol(ic)))*pri(ic)
2559 !CLOUD BASE MASS FLUX
2560  IF (.NOT.(akm .GE. 0.0 .OR. wlq .LT. 0.0)) THEN
2561 !================>>
2562  wfn = -((wfn-acr)/akm)
2563  x3 = rasal(ic)*trg*toki*wfn
2564  IF (x3 .GT. (prs(k+1)-prs(k))*(100.*pblfrac)) THEN
2565  wfn = (prs(k+1)-prs(k))*(100.*pblfrac)
2566  ELSE
2567  wfn = x3
2568  END IF
2569 !CUMULATIVE PRECIP AND CLOUD-BASE MASS FLUX FOR OUTPUT
2570  tem = wfn*gravi
2571  cll(ic) = cll(ic) + wlq*tem
2572  rmf(ic) = rmf(ic) + tem
2573  rmfd(ic) = rmfd(ic) + tem*eta(ic)
2574  DO l=ic+1,k
2575  rmfp(l) = tem*eta(l)
2576  rmfc(l) = rmfc(l) + rmfp(l)
2577  dllx(l) = dllx(l) + tem*dll0(l)
2578  IF (cvw(l) .GT. 0.0) THEN
2579  updfrp(l) = rmfp(l)*(ddt/daylen)*1000./(cvw(l)*prs&
2580 & (l))
2581  ELSE
2582  updfrp(l) = 0.0
2583  END IF
2584  clli(l) = cll0(l)/eta(l)
2585  updfrc(l) = updfrc(l) + updfrp(l)
2586  END DO
2587 !THETA AND Q CHANGE DUE TO CLOUD TYPE IC
2588  DO l=ic,k
2589  rns(l) = rns(l) + rnn(l)*tem
2590  gmhd(l) = 0.0_8
2591  gmh(l) = gmh(l)*wfn
2592  gmsd(l) = 0.0_8
2593  gms(l) = gms(l)*wfn
2594  qoid(l) = 0.0_8
2595  qoi(l) = qoi(l) + (gmh(l)-gms(l))*alhi
2596  poid(l) = 0.0_8
2597  poi(l) = poi(l) + gms(l)*pki(l)*cpi
2598  qstd(l) = 0.0_8
2599  qst(l) = qst(l) + gms(l)*bet(l)*cpi
2600  END DO
2601 !*FRICFAC*0.5
2602  wfn = wfn*0.5*1.0
2603 !DO_TRACERS
2604  tem = wfn*pri(k)
2605  DO itr=1,itrcr
2606  xcud(k, itr) = xcud(k, itr) + tem*(xoid(k-1, itr)-&
2607 & xoid(k, itr))
2608  xcu(k, itr) = xcu(k, itr) + tem*(xoi(k-1, itr)-xoi(k&
2609 & , itr))
2610  END DO
2611  DO itr=1,itrcr
2612  DO l=k-1,ic+1,-1
2613  tem = wfn*pri(l)
2614  xcud(l, itr) = xcud(l, itr) + tem*(eta(l)*(xoid(l-&
2615 & 1, itr)-xoid(l, itr))+eta(l+1)*(xoid(l, itr)-&
2616 & xoid(l+1, itr)))
2617  xcu(l, itr) = xcu(l, itr) + tem*((xoi(l-1, itr)-&
2618 & xoi(l, itr))*eta(l)+(xoi(l, itr)-xoi(l+1, itr))*&
2619 & eta(l+1))
2620  END DO
2621  END DO
2622  tem = wfn*pri(ic)
2623  DO itr=1,itrcr
2624  xcud(ic, itr) = xcud(ic, itr) + tem*(2.*(xhtd(itr)-(&
2625 & eta(ic)-eta(ic+1))*xoid(ic, itr))-eta(ic+1)*(xoid(&
2626 & ic, itr)+xoid(ic+1, itr)))
2627  xcu(ic, itr) = xcu(ic, itr) + (2.*(xht(itr)-xoi(ic, &
2628 & itr)*(eta(ic)-eta(ic+1)))-(xoi(ic, itr)+xoi(ic+1, &
2629 & itr))*eta(ic+1))*tem
2630  END DO
2631  DO itr=1,itrcr
2632  DO l=ic,k
2633  xoid(l, itr) = xoid(l, itr) + xcud(l, itr)
2634  xoi(l, itr) = xoi(l, itr) + xcu(l, itr)
2635  END DO
2636  END DO
2637 !End DO_TRACERS
2638 !CUMULUS FRICTION
2639  IF (fricfac .GT. 0.0) THEN
2640 !================>>
2641  wfn = wfn*fricfac*exp(-(alm/friclambda))
2642  tem = wfn*pri(k)
2643  ucu(k) = ucu(k) + tem*(uoi(k-1)-uoi(k))
2644  vcu(k) = vcu(k) + tem*(voi(k-1)-voi(k))
2645  DO l=k-1,ic+1,-1
2646  tem = wfn*pri(l)
2647  ucu(l) = ucu(l) + tem*((uoi(l-1)-uoi(l))*eta(l)+(&
2648 & uoi(l)-uoi(l+1))*eta(l+1))
2649  vcu(l) = vcu(l) + tem*((voi(l-1)-voi(l))*eta(l)+(&
2650 & voi(l)-voi(l+1))*eta(l+1))
2651  END DO
2652  tem = wfn*pri(ic)
2653  ucu(ic) = ucu(ic) + (2.*(uht-uoi(ic)*(eta(ic)-eta(ic&
2654 & +1)))-(uoi(ic)+uoi(ic+1))*eta(ic+1))*tem
2655  vcu(ic) = vcu(ic) + (2.*(vht-voi(ic)*(eta(ic)-eta(ic&
2656 & +1)))-(voi(ic)+voi(ic+1))*eta(ic+1))*tem
2657  DO l=ic,k
2658  uoi(l) = uoi(l) + ucu(l)
2659  voi(l) = voi(l) + vcu(l)
2660  END DO
2661  END IF
2662  END IF
2663  END IF
2664  END IF
2665  END IF
2666  END IF
2667  END IF
2668  END DO
2669 !CLOUD LOOP
2670  IF (sum(rmf(icmin:k)) .GT. 0.0) THEN
2671  tho(i, icmin:k-1) = poi(icmin:k-1)
2672  qho(i, icmin:k-1) = qoi(icmin:k-1)
2673  uho(i, icmin:k-1) = uoi(icmin:k-1)
2674  vho(i, icmin:k-1) = voi(icmin:k-1)
2675 !De-strap tendencies from RAS
2676  wght = wgt1(i, :)
2677 !Scale properly by layer masses
2678  wght0 = 0.
2679  DO l=k,k0
2680  wght0 = wght0 + wght(l)*(ple(i, l+1)-ple(i, l))
2681  END DO
2682  wght0 = (prs(k+1)-prs(k))/wght0
2683  wght = wght0*wght
2684  DO l=k,k0
2685  tho(i, l) = tho(i, l) + wght(l)*(poi(k)-poi_sv(k))
2686  qho(i, l) = qho(i, l) + wght(l)*(qoi(k)-qoi_sv(k))
2687  uho(i, l) = uho(i, l) + wght(l)*(uoi(k)-uoi_sv(k))
2688  vho(i, l) = vho(i, l) + wght(l)*(voi(k)-voi_sv(k))
2689  END DO
2690 !DO_TRACERS
2691  xhod(i, icmin:k-1, :) = xoid(icmin:k-1, :)
2692  xho(i, icmin:k-1, :) = xoi(icmin:k-1, :)
2693  DO itr=1,itrcr
2694  DO l=k,k0
2695  xhod(i, l, itr) = xhod(i, l, itr) + wght(l)*(xoid(k, itr)-&
2696 & xoi_svd(k, itr))
2697  xho(i, l, itr) = xho(i, l, itr) + wght(l)*(xoi(k, itr)-xoi_sv(&
2698 & k, itr))
2699  END DO
2700  END DO
2701  END IF
2702  END IF
2703 END SUBROUTINE rase_tracer_d
2704 
2705 END MODULE convection_tl
convection::acritn
subroutine, public acritn(PL, PLB, ACR, ACRITFAC)
Definition: convection.F90:641
convection::dqsat_ras
subroutine, public dqsat_ras(DQSi, QSSi, TEMP, PLO, lm, ESTBLX, CONS_H2OMW, CONS_AIRMW)
Definition: convection.F90:706
convection_tl
Definition: convection_tl.F90:1
convection_tl::rase0_d
subroutine, public rase0_d(idim, irun, k0, icmin, dt, cons_cp, cons_alhl, cons_grav, cons_rgas, cons_h2omw, cons_airmw, cons_vireps, seedras, sige, kcbl, wgt0, wgt1, frland, ts, tho, thod, qho, qhod, co_auto, ple, rasparams, estblx)
Definition: convection_tl.F90:1187
convection::dqsats_ras
subroutine, public dqsats_ras(DQSi, QSSi, TEMP, PLO, ESTBLX, CONS_H2OMW, CONS_AIRMW)
Definition: convection.F90:774
convection_tl::rase_d
subroutine, public rase_d(idim, irun, k0, icmin, dt, cons_cp, cons_alhl, cons_grav, cons_rgas, cons_h2omw, cons_airmw, cons_vireps, seedras, sige, kcbl, wgt0, wgt1, frland, ts, tho, thod, qho, qhod, uho, uhod, vho, vhod, co_auto, ple, clw, clwd, flxd, flxdd, cnv_prc3, cnv_prc3d, cnv_updfrc, cnv_updfrcd, rasparams, estblx)
Definition: convection_tl.F90:27
convection_tl::dqsats_ras_d
subroutine dqsats_ras_d(dqsi, dqsid, qssi, qssid, temp, tempd, plo, estblx, cons_h2omw, cons_airmw)
Definition: convection_tl.F90:1122
convection_tl::dqsat_ras_d
subroutine dqsat_ras_d(dqsi, dqsid, qssi, qssid, temp, tempd, plo, lm, estblx, cons_h2omw, cons_airmw)
Definition: convection_tl.F90:1051
convection::sundq3_ice
subroutine, public sundq3_ice(TEMP, RATE2, RATE3, TE1, F2, F3)
Definition: convection.F90:671
convection_tl::sundq3_ice_d
subroutine sundq3_ice_d(temp, tempd, rate2, rate3, te1, f2, f2d, f3)
Definition: convection_tl.F90:1014
max
#define max(a, b)
Definition: mosaic_util.h:33
min
#define min(a, b)
Definition: mosaic_util.h:32
convection_tl::rase_tracer_d
subroutine, public rase_tracer_d(idim, irun, k0, icmin, dt, cons_cp, cons_alhl, cons_grav, cons_rgas, cons_h2omw, cons_airmw, cons_vireps, seedras, sige, kcbl, wgt0, wgt1, frland, ts, thoin, qhoin, uhoin, vhoin, co_auto, ple, rasparams, estblx, itrcr, xho, xhod, fscav)
Definition: convection_tl.F90:1878