GCC Code Coverage Report


Directory: ./
File: phys/convect1.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 180 0.0%
Branches: 0 112 0.0%

Line Branch Exec Source
1
2 ! $Header$
3
4 ✗ SUBROUTINE convect1(len, nd, ndp1, noff, minorig, t, q, qs, u, v, p, ph, &
5 ✗ iflag, ft, fq, fu, fv, precip, cbmf, delt, ma)
6 ! .............................START PROLOGUE............................
7
8 ! SCCS IDENTIFICATION: @(#)convect1.f 1.1 04/21/00
9 ! 19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v
10
11 ! CONFIGURATION IDENTIFICATION: None
12
13 ! MODULE NAME: convect1
14
15 ! DESCRIPTION:
16
17 ! convect1 The Emanuel Cumulus Convection Scheme
18
19 ! CONTRACT NUMBER AND TITLE: None
20
21 ! REFERENCES: Programmers K. Emanuel (MIT), Timothy F. Hogan, M. Peng
22 ! (NRL)
23
24 ! CLASSIFICATION: Unclassified
25
26 ! RESTRICTIONS: None
27
28 ! COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
29
30 ! COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
31 ! Fortran 90: -O vector3,scalar3,task1,aggress,overindex -ei -r 2
32
33 ! LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
34
35 ! USAGE: call convect1(len,nd,noff,minorig,
36 ! & t,q,qs,u,v,
37 ! & p,ph,iflag,ft,
38 ! & fq,fu,fv,precip,cbmf,delt)
39
40 ! PARAMETERS:
41 ! Name Type Usage Description
42 ! ---------- ---------- ------- ----------------------------
43
44 ! len Integer Input first (i) dimension
45 ! nd Integer Input vertical (k) dimension
46 ! ndp1 Integer Input nd + 1
47 ! noff Integer Input integer limit for convection
48 ! (nd-noff)
49 ! minorig Integer Input First level of convection
50 ! t Real Input temperature
51 ! q Real Input specific hum
52 ! qs Real Input sat specific hum
53 ! u Real Input u-wind
54 ! v Real Input v-wind
55 ! p Real Input full level pressure
56 ! ph Real Input half level pressure
57 ! iflag Integer Output iflag on latitude strip
58 ! ft Real Output temp tend
59 ! fq Real Output spec hum tend
60 ! fu Real Output u-wind tend
61 ! fv Real Output v-wind tend
62 ! cbmf Real In/Out cumulus mass flux
63 ! delt Real Input time step
64 ! iflag Integer Output integer flag for Emanuel
65 ! conditions
66
67 ! COMMON BLOCKS:
68 ! Block Name Type Usage Notes
69 ! -------- -------- ---- ------ ------------------------
70
71 ! FILES: None
72
73 ! DATA BASES: None
74
75 ! NON-FILE INPUT/OUTPUT: None
76
77 ! ERROR CONDITIONS: None
78
79 ! ADDITIONAL COMMENTS: None
80
81 ! .................MAINTENANCE SECTION................................
82
83 ! MODULES CALLED:
84 ! Name Description
85 ! convect2 Emanuel cumulus convection tendency calculations
86 ! ------- ----------------------
87 ! LOCAL VARIABLES AND
88 ! STRUCTURES:
89 ! Name Type Description
90 ! ------- ------ -----------
91 ! See Comments Below
92
93 ! i Integer loop index
94 ! k Integer loop index
95
96 ! METHOD:
97
98 ! See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation
99 ! of a
100 ! convective scheme for use in climate models.
101
102 ! FILES: None
103
104 ! INCLUDE FILES: None
105
106 ! MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
107
108 ! ..............................END PROLOGUE.............................
109
110
111 USE dimphy
112 IMPLICIT NONE
113
114 INTEGER len
115 INTEGER nd
116 INTEGER ndp1
117 INTEGER noff
118 REAL t(len, nd)
119 REAL q(len, nd)
120 REAL qs(len, nd)
121 REAL u(len, nd)
122 REAL v(len, nd)
123 REAL p(len, nd)
124 REAL ph(len, ndp1)
125 INTEGER iflag(len)
126 REAL ft(len, nd)
127 REAL fq(len, nd)
128 REAL fu(len, nd)
129 REAL fv(len, nd)
130 REAL precip(len)
131 REAL cbmf(len)
132 REAL ma(len, nd)
133 INTEGER minorig
134 REAL delt, cpd, cpv, cl, rv, rd, lv0, g
135 REAL sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp
136 REAL alpha, entp, coeffs, coeffr, omtrain, cu
137
138 ! -------------------------------------------------------------------
139 ! --- ARGUMENTS
140 ! -------------------------------------------------------------------
141 ! --- On input:
142
143 ! t: Array of absolute temperature (K) of dimension ND, with first
144 ! index corresponding to lowest model level. Note that this array
145 ! will be altered by the subroutine if dry convective adjustment
146 ! occurs and if IPBL is not equal to 0.
147
148 ! q: Array of specific humidity (gm/gm) of dimension ND, with first
149 ! index corresponding to lowest model level. Must be defined
150 ! at same grid levels as T. Note that this array will be altered
151 ! if dry convective adjustment occurs and if IPBL is not equal to 0.
152
153 ! qs: Array of saturation specific humidity of dimension ND, with first
154 ! index corresponding to lowest model level. Must be defined
155 ! at same grid levels as T. Note that this array will be altered
156 ! if dry convective adjustment occurs and if IPBL is not equal to 0.
157
158 ! u: Array of zonal wind velocity (m/s) of dimension ND, witth first
159 ! index corresponding with the lowest model level. Defined at
160 ! same levels as T. Note that this array will be altered if
161 ! dry convective adjustment occurs and if IPBL is not equal to 0.
162
163 ! v: Same as u but for meridional velocity.
164
165 ! tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
166 ! where NTRA is the number of different tracers. If no
167 ! convective tracer transport is needed, define a dummy
168 ! input array of dimension (ND,1). Tracers are defined at
169 ! same vertical levels as T. Note that this array will be altered
170 ! if dry convective adjustment occurs and if IPBL is not equal to 0.
171
172 ! p: Array of pressure (mb) of dimension ND, with first
173 ! index corresponding to lowest model level. Must be defined
174 ! at same grid levels as T.
175
176 ! ph: Array of pressure (mb) of dimension ND+1, with first index
177 ! corresponding to lowest level. These pressures are defined at
178 ! levels intermediate between those of P, T, Q and QS. The first
179 ! value of PH should be greater than (i.e. at a lower level than)
180 ! the first value of the array P.
181
182 ! nl: The maximum number of levels to which convection can penetrate, plus
183 ! 1.
184 ! NL MUST be less than or equal to ND-1.
185
186 ! delt: The model time step (sec) between calls to CONVECT
187
188 ! ----------------------------------------------------------------------------
189 ! --- On Output:
190
191 ! iflag: An output integer whose value denotes the following:
192 ! VALUE INTERPRETATION
193 ! ----- --------------
194 ! 0 Moist convection occurs.
195 ! 1 Moist convection occurs, but a CFL condition
196 ! on the subsidence warming is violated. This
197 ! does not cause the scheme to terminate.
198 ! 2 Moist convection, but no precip because ep(inb) lt 0.0001
199 ! 3 No moist convection because new cbmf is 0 and old cbmf is 0.
200 ! 4 No moist convection; atmosphere is not
201 ! unstable
202 ! 6 No moist convection because ihmin le minorig.
203 ! 7 No moist convection because unreasonable
204 ! parcel level temperature or specific humidity.
205 ! 8 No moist convection: lifted condensation
206 ! level is above the 200 mb level.
207 ! 9 No moist convection: cloud base is higher
208 ! then the level NL-1.
209
210 ! ft: Array of temperature tendency (K/s) of dimension ND, defined at
211 ! same
212 ! grid levels as T, Q, QS and P.
213
214 ! fq: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
215 ! defined at same grid levels as T, Q, QS and P.
216
217 ! fu: Array of forcing of zonal velocity (m/s^2) of dimension ND,
218 ! defined at same grid levels as T.
219
220 ! fv: Same as FU, but for forcing of meridional velocity.
221
222 ! ftra: Array of forcing of tracer content, in tracer mixing ratio per
223 ! second, defined at same levels as T. Dimensioned (ND,NTRA).
224
225 ! precip: Scalar convective precipitation rate (mm/day).
226
227 ! wd: A convective downdraft velocity scale. For use in surface
228 ! flux parameterizations. See convect.ps file for details.
229
230 ! tprime: A convective downdraft temperature perturbation scale (K).
231 ! For use in surface flux parameterizations. See convect.ps
232 ! file for details.
233
234 ! qprime: A convective downdraft specific humidity
235 ! perturbation scale (gm/gm).
236 ! For use in surface flux parameterizations. See convect.ps
237 ! file for details.
238
239 ! cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
240 ! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
241 ! ITS NEXT CALL. That is, the value of CBMF must be "remembered"
242 ! by the calling program between calls to CONVECT.
243
244 ! det: Array of detrainment mass flux of dimension ND.
245
246 ! -------------------------------------------------------------------
247
248 ! Local arrays
249
250 INTEGER nl
251 INTEGER nlp
252 INTEGER nlm
253 INTEGER i, k, n
254 REAL delti
255 REAL rowl
256 REAL clmcpv
257 REAL clmcpd
258 REAL cpdmcp
259 REAL cpvmcpd
260 REAL eps
261 REAL epsi
262 REAL epsim1
263 REAL ginv
264 REAL hrd
265 REAL prccon1
266 INTEGER icbmax
267 ✗ REAL lv(klon, klev)
268 ✗ REAL cpn(klon, klev)
269 ✗ REAL cpx(klon, klev)
270 ✗ REAL tv(klon, klev)
271 ✗ REAL gz(klon, klev)
272 ✗ REAL hm(klon, klev)
273 ✗ REAL h(klon, klev)
274 ✗ REAL work(klon)
275 ✗ INTEGER ihmin(klon)
276 ✗ INTEGER nk(klon)
277 ✗ REAL rh(klon)
278 ✗ REAL chi(klon)
279 ✗ REAL plcl(klon)
280 ✗ INTEGER icb(klon)
281 ✗ REAL tnk(klon)
282 ✗ REAL qnk(klon)
283 ✗ REAL gznk(klon)
284 ✗ REAL pnk(klon)
285 ✗ REAL qsnk(klon)
286 ✗ REAL ticb(klon)
287 ✗ REAL gzicb(klon)
288 ✗ REAL tp(klon, klev)
289 ✗ REAL tvp(klon, klev)
290 ✗ REAL clw(klon, klev)
291
292 ✗ REAL ah0(klon), cpp(klon)
293 REAL tg, qg, s, alv, tc, ahg, denom, es, rg
294
295 INTEGER ncum
296 ✗ INTEGER idcum(klon)
297
298 ✗ cpd = 1005.7
299 ✗ cpv = 1870.0
300 ✗ cl = 4190.0
301 ✗ rv = 461.5
302 ✗ rd = 287.04
303 ✗ lv0 = 2.501E6
304 ✗ g = 9.8
305
306 ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
307 ! *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- ***
308 ! *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO ***
309 ! *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY ***
310 ! *** BETWEEN 0 C AND TLCRIT) ***
311 ! *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT ***
312 ! *** FORMULATION ***
313 ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
314 ! *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
315 ! *** OF CLOUD ***
316 ! *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN ***
317 ! *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW ***
318 ! *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
319 ! *** OF RAIN ***
320 ! *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
321 ! *** OF SNOW ***
322 ! *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM ***
323 ! *** TRANSPORT ***
324 ! *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION ***
325 ! *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC ***
326 ! *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF ***
327 ! *** APPROACH TO QUASI-EQUILIBRIUM ***
328 ! *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) ***
329 ! *** (DAMP MUST BE LESS THAN 1) ***
330
331 ✗ sigs = 0.12
332 ✗ sigd = 0.05
333 ✗ elcrit = 0.0011
334 ✗ tlcrit = -55.0
335 ✗ omtsnow = 5.5
336 ✗ dtmax = 0.9
337 ✗ damp = 0.1
338 ✗ alpha = 0.2
339 ✗ entp = 1.5
340 ✗ coeffs = 0.8
341 ✗ coeffr = 1.0
342 ✗ omtrain = 50.0
343
344 ✗ cu = 0.70
345 damp = 0.1
346
347
348 ! Define nl, nlp, nlm, and delti
349
350 ✗ nl = nd - noff
351 ✗ nlp = nl + 1
352 ✗ nlm = nl - 1
353 delti = 1.0/delt
354
355 ! -------------------------------------------------------------------
356 ! --- SET CONSTANTS
357 ! -------------------------------------------------------------------
358
359 rowl = 1000.0
360 clmcpv = cl - cpv
361 clmcpd = cl - cpd
362 cpdmcp = cpd - cpv
363 cpvmcpd = cpv - cpd
364 eps = rd/rv
365 epsi = 1.0/eps
366 epsim1 = epsi - 1.0
367 ginv = 1.0/g
368 hrd = 0.5*rd
369 prccon1 = 86400.0*1000.0/(rowl*g)
370
371 ! dtmax is the maximum negative temperature perturbation.
372
373 ! =====================================================================
374 ! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS
375 ! =====================================================================
376
377 ✗ DO k = 1, nd
378 ✗ DO i = 1, len
379 ✗ ft(i, k) = 0.0
380 ✗ fq(i, k) = 0.0
381 ✗ fu(i, k) = 0.0
382 ✗ fv(i, k) = 0.0
383 ✗ tvp(i, k) = 0.0
384 ✗ tp(i, k) = 0.0
385 ✗ clw(i, k) = 0.0
386 ✗ gz(i, k) = 0.
387 END DO
388 END DO
389 ✗ DO i = 1, len
390 ✗ precip(i) = 0.0
391 ✗ iflag(i) = 0
392 END DO
393
394 ! =====================================================================
395 ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
396 ! =====================================================================
397 ✗ DO k = 1, nl + 1
398 ✗ DO i = 1, len
399 ✗ lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
400 ✗ cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
401 ✗ cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
402 ✗ tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
403 END DO
404 END DO
405
406 ! gz = phi at the full levels (same as p).
407
408 ✗ DO i = 1, len
409 ✗ gz(i, 1) = 0.0
410 END DO
411 ✗ DO k = 2, nlp
412 ✗ DO i = 1, len
413 gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
414 ✗ k)
415 END DO
416 END DO
417
418 ! h = phi + cpT (dry static energy).
419 ! hm = phi + cp(T-Tbase)+Lq
420
421 ✗ DO k = 1, nlp
422 ✗ DO i = 1, len
423 ✗ h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
424 ✗ hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
425 END DO
426 END DO
427
428 ! -------------------------------------------------------------------
429 ! --- Find level of minimum moist static energy
430 ! --- If level of minimum moist static energy coincides with
431 ! --- or is lower than minimum allowable parcel origin level,
432 ! --- set iflag to 6.
433 ! -------------------------------------------------------------------
434 ✗ DO i = 1, len
435 ✗ work(i) = 1.0E12
436 ✗ ihmin(i) = nl
437 END DO
438 ✗ DO k = 2, nlp
439 ✗ DO i = 1, len
440 ✗ IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
441 ✗ work(i) = hm(i, k)
442 ✗ ihmin(i) = k
443 END IF
444 END DO
445 END DO
446 ✗ DO i = 1, len
447 ✗ ihmin(i) = min(ihmin(i), nlm)
448 ✗ IF (ihmin(i)<=minorig) THEN
449 ✗ iflag(i) = 6
450 END IF
451 END DO
452
453 ! -------------------------------------------------------------------
454 ! --- Find that model level below the level of minimum moist static
455 ! --- energy that has the maximum value of moist static energy
456 ! -------------------------------------------------------------------
457
458 ✗ DO i = 1, len
459 ✗ work(i) = hm(i, minorig)
460 ✗ nk(i) = minorig
461 END DO
462 ✗ DO k = minorig + 1, nl
463 ✗ DO i = 1, len
464 ✗ IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
465 ✗ work(i) = hm(i, k)
466 ✗ nk(i) = k
467 END IF
468 END DO
469 END DO
470 ! -------------------------------------------------------------------
471 ! --- Check whether parcel level temperature and specific humidity
472 ! --- are reasonable
473 ! -------------------------------------------------------------------
474 ✗ DO i = 1, len
475 IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
476 ✗ 400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
477 END DO
478 ! -------------------------------------------------------------------
479 ! --- Calculate lifted condensation level of air at parcel origin level
480 ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
481 ! -------------------------------------------------------------------
482 ✗ DO i = 1, len
483 ✗ tnk(i) = t(i, nk(i))
484 ✗ qnk(i) = q(i, nk(i))
485 ✗ gznk(i) = gz(i, nk(i))
486 ✗ pnk(i) = p(i, nk(i))
487 ✗ qsnk(i) = qs(i, nk(i))
488
489 ✗ rh(i) = qnk(i)/qsnk(i)
490 ✗ rh(i) = min(1.0, rh(i))
491 ✗ chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
492 ✗ plcl(i) = pnk(i)*(rh(i)**chi(i))
493 ✗ IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
494 ✗ ) = 8
495 END DO
496 ! -------------------------------------------------------------------
497 ! --- Calculate first level above lcl (=icb)
498 ! -------------------------------------------------------------------
499 ✗ DO i = 1, len
500 ✗ icb(i) = nlm
501 END DO
502
503 ✗ DO k = minorig, nl
504 ✗ DO i = 1, len
505 ✗ IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
506 END DO
507 END DO
508
509 ✗ DO i = 1, len
510 ✗ IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
511 END DO
512
513 ! Compute icbmax.
514
515 icbmax = 2
516 ✗ DO i = 1, len
517 ✗ icbmax = max(icbmax, icb(i))
518 END DO
519
520 ! -------------------------------------------------------------------
521 ! --- Calculates the lifted parcel virtual temperature at nk,
522 ! --- the actual temperature, and the adiabatic
523 ! --- liquid water content. The procedure is to solve the equation.
524 ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
525 ! -------------------------------------------------------------------
526
527 ✗ DO i = 1, len
528 ✗ tnk(i) = t(i, nk(i))
529 ✗ qnk(i) = q(i, nk(i))
530 ✗ gznk(i) = gz(i, nk(i))
531 ✗ ticb(i) = t(i, icb(i))
532 ✗ gzicb(i) = gz(i, icb(i))
533 END DO
534
535 ! *** Calculate certain parcel quantities, including static energy ***
536
537 ✗ DO i = 1, len
538 ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
539 ✗ 273.15)) + gznk(i)
540 ✗ cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
541 END DO
542
543 ! *** Calculate lifted parcel quantities below cloud base ***
544
545 ✗ DO k = minorig, icbmax - 1
546 ✗ DO i = 1, len
547 ✗ tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
548 ✗ tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
549 END DO
550 END DO
551
552 ! *** Find lifted parcel quantities above cloud base ***
553
554 ✗ DO i = 1, len
555 ✗ tg = ticb(i)
556 ✗ qg = qs(i, icb(i))
557 ✗ alv = lv0 - clmcpv*(ticb(i)-273.15)
558
559 ! First iteration.
560
561 ✗ s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
562 ✗ s = 1./s
563 ✗ ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
564 ✗ tg = tg + s*(ah0(i)-ahg)
565 ✗ tg = max(tg, 35.0)
566 ✗ tc = tg - 273.15
567 ✗ denom = 243.5 + tc
568 ✗ IF (tc>=0.0) THEN
569 ✗ es = 6.112*exp(17.67*tc/denom)
570 ELSE
571 ✗ es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
572 END IF
573 ✗ qg = eps*es/(p(i,icb(i))-es*(1.-eps))
574
575 ! Second iteration.
576
577 ✗ s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
578 ✗ s = 1./s
579 ✗ ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
580 ✗ tg = tg + s*(ah0(i)-ahg)
581 ✗ tg = max(tg, 35.0)
582 ✗ tc = tg - 273.15
583 ✗ denom = 243.5 + tc
584 ✗ IF (tc>=0.0) THEN
585 ✗ es = 6.112*exp(17.67*tc/denom)
586 ELSE
587 ✗ es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
588 END IF
589 ✗ qg = eps*es/(p(i,icb(i))-es*(1.-eps))
590
591 alv = lv0 - clmcpv*(ticb(i)-273.15)
592 ✗ tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
593 ✗ clw(i, icb(i)) = qnk(i) - qg
594 ✗ clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
595 ✗ rg = qg/(1.-qnk(i))
596 ✗ tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
597 END DO
598
599 ✗ DO k = minorig, icbmax
600 ✗ DO i = 1, len
601 ✗ tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
602 END DO
603 END DO
604
605 ! -------------------------------------------------------------------
606 ! --- Test for instability.
607 ! --- If there was no convection at last time step and parcel
608 ! --- is stable at icb, then set iflag to 4.
609 ! -------------------------------------------------------------------
610
611 ✗ DO i = 1, len
612 ✗ IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
613 ✗ icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
614 END DO
615
616 ! =====================================================================
617 ! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
618 ! =====================================================================
619
620 ✗ ncum = 0
621 ✗ DO i = 1, len
622 ✗ IF (iflag(i)==0) THEN
623 ✗ ncum = ncum + 1
624 ✗ idcum(ncum) = i
625 END IF
626 END DO
627
628 ! Call convect2, which compresses the points and computes the heating,
629 ! moistening, velocity mixing, and precipiation.
630
631 ! print*,'cpd avant convect2 ',cpd
632 ✗ IF (ncum>0) THEN
633 CALL convect2(ncum, idcum, len, nd, ndp1, nl, minorig, nk, icb, t, q, qs, &
634 u, v, gz, tv, tp, tvp, clw, h, lv, cpn, p, ph, ft, fq, fu, fv, tnk, &
635 qnk, gznk, plcl, precip, cbmf, iflag, delt, cpd, cpv, cl, rv, rd, lv0, &
636 g, sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp, alpha, entp, &
637 ✗ coeffs, coeffr, omtrain, cu, ma)
638 END IF
639
640 ✗ RETURN
641 END SUBROUTINE convect1
642