GCC Code Coverage Report


Directory: ./
File: phys/cv30_routines.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 1077 0.0%
Branches: 0 942 0.0%

Line Branch Exec Source
1
2 ! $Id: cv30_routines.F90 3577 2019-10-07 13:30:19Z fhourdin $
3
4
5
6 SUBROUTINE cv30_param(nd, delt)
7 IMPLICIT NONE
8
9 ! ------------------------------------------------------------
10 ! Set parameters for convectL for iflag_con = 3
11 ! ------------------------------------------------------------
12
13
14 ! *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
15 ! *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO ***
16 ! *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
17 ! *** EFFICIENCY IS ASSUMED TO BE UNITY ***
18 ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
19 ! *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
20 ! *** OF CLOUD ***
21
22 ! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
23 ! *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
24 ! *** APPROACH TO QUASI-EQUILIBRIUM ***
25 ! *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
26 ! *** (BETA MUST BE LESS THAN OR EQUAL TO 1) ***
27
28 ! *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
29 ! *** APPROACH TO QUASI-EQUILIBRIUM ***
30 ! *** IT MUST BE LESS THAN 0 ***
31
32 include "cv30param.h"
33 include "conema3.h"
34
35 INTEGER nd
36 REAL delt ! timestep (seconds)
37
38 ! noff: integer limit for convection (nd-noff)
39 ! minorig: First level of convection
40
41 ! -- limit levels for convection:
42
43 noff = 1
44 minorig = 1
45 nl = nd - noff
46 nlp = nl + 1
47 nlm = nl - 1
48
49 ! -- "microphysical" parameters:
50
51 sigd = 0.01
52 spfac = 0.15
53 pbcrit = 150.0
54 ptcrit = 500.0
55 ! IM cf. FH epmax = 0.993
56
57 omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
58
59 ! -- misc:
60
61 dtovsh = -0.2 ! dT for overshoot
62 dpbase = -40. ! definition cloud base (400m above LCL)
63 dttrig = 5. ! (loose) condition for triggering
64
65 ! -- rate of approach to quasi-equilibrium:
66
67 dtcrit = -2.0
68 tau = 8000.
69 beta = 1.0 - delt/tau
70 alpha = 1.5E-3*delt/tau
71 ! increase alpha to compensate W decrease:
72 alpha = alpha*1.5
73
74 ! -- interface cloud parameterization:
75
76 delta = 0.01 ! cld
77
78 ! -- interface with boundary-layer (gust factor): (sb)
79
80 betad = 10.0 ! original value (from convect 4.3)
81
82 RETURN
83 END SUBROUTINE cv30_param
84
85 SUBROUTINE cv30_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm, &
86 th)
87 IMPLICIT NONE
88
89 ! =====================================================================
90 ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
91 ! "ori": from convect4.3 (vectorized)
92 ! "convect3": to be exactly consistent with convect3
93 ! =====================================================================
94
95 ! inputs:
96 INTEGER len, nd, ndp1
97 REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
98
99 ! outputs:
100 REAL lv(len, nd), cpn(len, nd), tv(len, nd)
101 REAL gz(len, nd), h(len, nd), hm(len, nd)
102 REAL th(len, nd)
103
104 ! local variables:
105 INTEGER k, i
106 REAL rdcp
107 REAL tvx, tvy ! convect3
108 REAL cpx(len, nd)
109
110 include "cvthermo.h"
111 include "cv30param.h"
112
113
114 ! ori do 110 k=1,nlp
115 DO k = 1, nl ! convect3
116 DO i = 1, len
117 ! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
118 lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
119 cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
120 cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
121 ! ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
122 tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
123 rdcp = (rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i, k)
124 th(i, k) = t(i, k)*(1000.0/p(i,k))**rdcp
125 END DO
126 END DO
127
128 ! gz = phi at the full levels (same as p).
129
130 DO i = 1, len
131 gz(i, 1) = 0.0
132 END DO
133 ! ori do 140 k=2,nlp
134 DO k = 2, nl ! convect3
135 DO i = 1, len
136 tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
137 tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
138 gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy) & !convect3
139 *(p(i,k-1)-p(i,k))/ph(i, k) !convect3
140
141 ! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
142 ! ori & *(p(i,k-1)-p(i,k))/ph(i,k)
143 END DO
144 END DO
145
146 ! h = phi + cpT (dry static energy).
147 ! hm = phi + cp(T-Tbase)+Lq
148
149 ! ori do 170 k=1,nlp
150 DO k = 1, nl ! convect3
151 DO i = 1, len
152 h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
153 hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
154 END DO
155 END DO
156
157 RETURN
158 END SUBROUTINE cv30_prelim
159
160 SUBROUTINE cv30_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, &
161 iflag, tnk, qnk, gznk, plcl)
162 IMPLICIT NONE
163
164 ! ================================================================
165 ! Purpose: CONVECTIVE FEED
166
167 ! Main differences with cv_feed:
168 ! - ph added in input
169 ! - here, nk(i)=minorig
170 ! - icb defined differently (plcl compared with ph instead of p)
171
172 ! Main differences with convect3:
173 ! - we do not compute dplcldt and dplcldr of CLIFT anymore
174 ! - values iflag different (but tests identical)
175 ! - A,B explicitely defined (!...)
176 ! ================================================================
177
178 include "cv30param.h"
179
180 ! inputs:
181 INTEGER len, nd
182 REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
183 REAL hm(len, nd), gz(len, nd)
184 REAL ph(len, nd+1)
185
186 ! outputs:
187 INTEGER iflag(len), nk(len), icb(len), icbmax
188 REAL tnk(len), qnk(len), gznk(len), plcl(len)
189
190 ! local variables:
191 INTEGER i, k
192 INTEGER ihmin(len)
193 REAL work(len)
194 REAL pnk(len), qsnk(len), rh(len), chi(len)
195 REAL a, b ! convect3
196 ! ym
197 plcl = 0.0
198 ! @ !-------------------------------------------------------------------
199 ! @ ! --- Find level of minimum moist static energy
200 ! @ ! --- If level of minimum moist static energy coincides with
201 ! @ ! --- or is lower than minimum allowable parcel origin level,
202 ! @ ! --- set iflag to 6.
203 ! @ !-------------------------------------------------------------------
204 ! @
205 ! @ do 180 i=1,len
206 ! @ work(i)=1.0e12
207 ! @ ihmin(i)=nl
208 ! @ 180 continue
209 ! @ do 200 k=2,nlp
210 ! @ do 190 i=1,len
211 ! @ if((hm(i,k).lt.work(i)).and.
212 ! @ & (hm(i,k).lt.hm(i,k-1)))then
213 ! @ work(i)=hm(i,k)
214 ! @ ihmin(i)=k
215 ! @ endif
216 ! @ 190 continue
217 ! @ 200 continue
218 ! @ do 210 i=1,len
219 ! @ ihmin(i)=min(ihmin(i),nlm)
220 ! @ if(ihmin(i).le.minorig)then
221 ! @ iflag(i)=6
222 ! @ endif
223 ! @ 210 continue
224 ! @ c
225 ! @ !-------------------------------------------------------------------
226 ! @ ! --- Find that model level below the level of minimum moist static
227 ! @ ! --- energy that has the maximum value of moist static energy
228 ! @ !-------------------------------------------------------------------
229 ! @
230 ! @ do 220 i=1,len
231 ! @ work(i)=hm(i,minorig)
232 ! @ nk(i)=minorig
233 ! @ 220 continue
234 ! @ do 240 k=minorig+1,nl
235 ! @ do 230 i=1,len
236 ! @ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
237 ! @ work(i)=hm(i,k)
238 ! @ nk(i)=k
239 ! @ endif
240 ! @ 230 continue
241 ! @ 240 continue
242
243 ! -------------------------------------------------------------------
244 ! --- Origin level of ascending parcels for convect3:
245 ! -------------------------------------------------------------------
246
247 DO i = 1, len
248 nk(i) = minorig
249 END DO
250
251 ! -------------------------------------------------------------------
252 ! --- Check whether parcel level temperature and specific humidity
253 ! --- are reasonable
254 ! -------------------------------------------------------------------
255 DO i = 1, len
256 IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0)) & ! @ & .or.(
257 ! p(i,ihmin(i)).lt.400.0
258 ! ) )
259 .AND. (iflag(i)==0)) iflag(i) = 7
260 END DO
261 ! -------------------------------------------------------------------
262 ! --- Calculate lifted condensation level of air at parcel origin level
263 ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
264 ! -------------------------------------------------------------------
265
266 a = 1669.0 ! convect3
267 b = 122.0 ! convect3
268
269 DO i = 1, len
270
271 IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
272
273 tnk(i) = t(i, nk(i))
274 qnk(i) = q(i, nk(i))
275 gznk(i) = gz(i, nk(i))
276 pnk(i) = p(i, nk(i))
277 qsnk(i) = qs(i, nk(i))
278
279 rh(i) = qnk(i)/qsnk(i)
280 ! ori rh(i)=min(1.0,rh(i)) ! removed for convect3
281 ! ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
282 chi(i) = tnk(i)/(a-b*rh(i)-tnk(i)) ! convect3
283 plcl(i) = pnk(i)*(rh(i)**chi(i))
284 IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag &
285 (i) = 8
286
287 END IF ! iflag=7
288
289 END DO
290
291 ! -------------------------------------------------------------------
292 ! --- Calculate first level above lcl (=icb)
293 ! -------------------------------------------------------------------
294
295 ! @ do 270 i=1,len
296 ! @ icb(i)=nlm
297 ! @ 270 continue
298 ! @c
299 ! @ do 290 k=minorig,nl
300 ! @ do 280 i=1,len
301 ! @ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
302 ! @ & icb(i)=min(icb(i),k)
303 ! @ 280 continue
304 ! @ 290 continue
305 ! @c
306 ! @ do 300 i=1,len
307 ! @ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
308 ! @ 300 continue
309
310 DO i = 1, len
311 icb(i) = nlm
312 END DO
313
314 ! la modification consiste a comparer plcl a ph et non a p:
315 ! icb est defini par : ph(icb)<plcl<ph(icb-1)
316 ! @ do 290 k=minorig,nl
317 DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
318 DO i = 1, len
319 IF (ph(i,k)<plcl(i)) icb(i) = min(icb(i), k)
320 END DO
321 END DO
322
323 DO i = 1, len
324 ! @ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
325 IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
326 END DO
327
328 DO i = 1, len
329 icb(i) = icb(i) - 1 ! icb sup ou egal a 2
330 END DO
331
332 ! Compute icbmax.
333
334 icbmax = 2
335 DO i = 1, len
336 ! ! icbmax=max(icbmax,icb(i))
337 IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
338 END DO
339
340 RETURN
341 END SUBROUTINE cv30_feed
342
343 SUBROUTINE cv30_undilute1(len, nd, t, q, qs, gz, plcl, p, nk, icb, tp, tvp, &
344 clw, icbs)
345 IMPLICIT NONE
346
347 ! ----------------------------------------------------------------
348 ! Equivalent de TLIFT entre NK et ICB+1 inclus
349
350 ! Differences with convect4:
351 ! - specify plcl in input
352 ! - icbs is the first level above LCL (may differ from icb)
353 ! - in the iterations, used x(icbs) instead x(icb)
354 ! - many minor differences in the iterations
355 ! - tvp is computed in only one time
356 ! - icbs: first level above Plcl (IMIN de TLIFT) in output
357 ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
358 ! ----------------------------------------------------------------
359
360 include "cvthermo.h"
361 include "cv30param.h"
362
363 ! inputs:
364 INTEGER len, nd
365 INTEGER nk(len), icb(len)
366 REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
367 REAL p(len, nd)
368 REAL plcl(len) ! convect3
369
370 ! outputs:
371 REAL tp(len, nd), tvp(len, nd), clw(len, nd)
372
373 ! local variables:
374 INTEGER i, k
375 INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
376 REAL tg, qg, alv, s, ahg, tc, denom, es, rg
377 REAL ah0(len), cpp(len)
378 REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
379 REAL qsicb(len) ! convect3
380 REAL cpinv(len) ! convect3
381
382 ! -------------------------------------------------------------------
383 ! --- Calculates the lifted parcel virtual temperature at nk,
384 ! --- the actual temperature, and the adiabatic
385 ! --- liquid water content. The procedure is to solve the equation.
386 ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
387 ! -------------------------------------------------------------------
388
389 DO i = 1, len
390 tnk(i) = t(i, nk(i))
391 qnk(i) = q(i, nk(i))
392 gznk(i) = gz(i, nk(i))
393 ! ori ticb(i)=t(i,icb(i))
394 ! ori gzicb(i)=gz(i,icb(i))
395 END DO
396
397 ! *** Calculate certain parcel quantities, including static energy ***
398
399 DO i = 1, len
400 ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
401 273.15)) + gznk(i)
402 cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
403 cpinv(i) = 1./cpp(i)
404 END DO
405
406 ! *** Calculate lifted parcel quantities below cloud base ***
407
408 DO i = 1, len !convect3
409 icb1(i) = min(max(icb(i), 2), nl)
410 ! if icb is below LCL, start loop at ICB+1:
411 ! (icbs est le premier niveau au-dessus du LCL)
412 icbs(i) = icb1(i) !convect3
413 IF (plcl(i)<p(i,icb1(i))) THEN
414 icbs(i) = min(icbs(i)+1, nl) !convect3
415 END IF
416 END DO !convect3
417
418 DO i = 1, len !convect3
419 ticb(i) = t(i, icbs(i)) !convect3
420 gzicb(i) = gz(i, icbs(i)) !convect3
421 qsicb(i) = qs(i, icbs(i)) !convect3
422 END DO !convect3
423
424
425 ! Re-compute icbsmax (icbsmax2): !convect3
426 ! !convect3
427 icbsmax2 = 2 !convect3
428 DO i = 1, len !convect3
429 icbsmax2 = max(icbsmax2, icbs(i)) !convect3
430 END DO !convect3
431
432 ! initialization outputs:
433
434 DO k = 1, icbsmax2 ! convect3
435 DO i = 1, len ! convect3
436 tp(i, k) = 0.0 ! convect3
437 tvp(i, k) = 0.0 ! convect3
438 clw(i, k) = 0.0 ! convect3
439 END DO ! convect3
440 END DO ! convect3
441
442 ! tp and tvp below cloud base:
443
444 DO k = minorig, icbsmax2 - 1
445 DO i = 1, len
446 tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
447 tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
448 END DO
449 END DO
450
451 ! *** Find lifted parcel quantities above cloud base ***
452
453 DO i = 1, len
454 tg = ticb(i)
455 ! ori qg=qs(i,icb(i))
456 qg = qsicb(i) ! convect3
457 ! debug alv=lv0-clmcpv*(ticb(i)-t0)
458 alv = lv0 - clmcpv*(ticb(i)-273.15)
459
460 ! First iteration.
461
462 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
463 s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
464 +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
465 s = 1./s
466 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
467 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
468 tg = tg + s*(ah0(i)-ahg)
469 ! ori tg=max(tg,35.0)
470 ! debug tc=tg-t0
471 tc = tg - 273.15
472 denom = 243.5 + tc
473 denom = max(denom, 1.0) ! convect3
474 ! ori if(tc.ge.0.0)then
475 es = 6.112*exp(17.67*tc/denom)
476 ! ori else
477 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
478 ! ori endif
479 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
480 qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
481
482 ! Second iteration.
483
484
485 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
486 ! ori s=1./s
487 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
488 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
489 tg = tg + s*(ah0(i)-ahg)
490 ! ori tg=max(tg,35.0)
491 ! debug tc=tg-t0
492 tc = tg - 273.15
493 denom = 243.5 + tc
494 denom = max(denom, 1.0) ! convect3
495 ! ori if(tc.ge.0.0)then
496 es = 6.112*exp(17.67*tc/denom)
497 ! ori else
498 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
499 ! ori end if
500 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
501 qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
502
503 alv = lv0 - clmcpv*(ticb(i)-273.15)
504
505 ! ori c approximation here:
506 ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
507 ! ori & -gz(i,icb(i))-alv*qg)/cpd
508
509 ! convect3: no approximation:
510 tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i))
511
512 ! ori clw(i,icb(i))=qnk(i)-qg
513 ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
514 clw(i, icbs(i)) = qnk(i) - qg
515 clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
516
517 rg = qg/(1.-qnk(i))
518 ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
519 ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
520 tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
521
522 END DO
523
524 ! ori do 380 k=minorig,icbsmax2
525 ! ori do 370 i=1,len
526 ! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
527 ! ori 370 continue
528 ! ori 380 continue
529
530
531 ! -- The following is only for convect3:
532
533 ! * icbs is the first level above the LCL:
534 ! if plcl<p(icb), then icbs=icb+1
535 ! if plcl>p(icb), then icbs=icb
536
537 ! * the routine above computes tvp from minorig to icbs (included).
538
539 ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
540 ! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
541
542 ! * therefore, in the case icbs=icb, we compute tvp at level icb+1
543 ! (tvp at other levels will be computed in cv3_undilute2.F)
544
545
546 DO i = 1, len
547 ticb(i) = t(i, icb(i)+1)
548 gzicb(i) = gz(i, icb(i)+1)
549 qsicb(i) = qs(i, icb(i)+1)
550 END DO
551
552 DO i = 1, len
553 tg = ticb(i)
554 qg = qsicb(i) ! convect3
555 ! debug alv=lv0-clmcpv*(ticb(i)-t0)
556 alv = lv0 - clmcpv*(ticb(i)-273.15)
557
558 ! First iteration.
559
560 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
561 s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
562 +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
563 s = 1./s
564 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
565 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
566 tg = tg + s*(ah0(i)-ahg)
567 ! ori tg=max(tg,35.0)
568 ! debug tc=tg-t0
569 tc = tg - 273.15
570 denom = 243.5 + tc
571 denom = max(denom, 1.0) ! convect3
572 ! ori if(tc.ge.0.0)then
573 es = 6.112*exp(17.67*tc/denom)
574 ! ori else
575 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
576 ! ori endif
577 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
578 qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
579
580 ! Second iteration.
581
582
583 ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
584 ! ori s=1./s
585 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
586 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
587 tg = tg + s*(ah0(i)-ahg)
588 ! ori tg=max(tg,35.0)
589 ! debug tc=tg-t0
590 tc = tg - 273.15
591 denom = 243.5 + tc
592 denom = max(denom, 1.0) ! convect3
593 ! ori if(tc.ge.0.0)then
594 es = 6.112*exp(17.67*tc/denom)
595 ! ori else
596 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
597 ! ori end if
598 ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
599 qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
600
601 alv = lv0 - clmcpv*(ticb(i)-273.15)
602
603 ! ori c approximation here:
604 ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
605 ! ori & -gz(i,icb(i))-alv*qg)/cpd
606
607 ! convect3: no approximation:
608 tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
609
610 ! ori clw(i,icb(i))=qnk(i)-qg
611 ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
612 clw(i, icb(i)+1) = qnk(i) - qg
613 clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
614
615 rg = qg/(1.-qnk(i))
616 ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
617 ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
618 tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
619
620 END DO
621
622 RETURN
623 END SUBROUTINE cv30_undilute1
624
625 SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, &
626 iflag, sig, w0)
627 IMPLICIT NONE
628
629 ! -------------------------------------------------------------------
630 ! --- TRIGGERING
631
632 ! - computes the cloud base
633 ! - triggering (crude in this version)
634 ! - relaxation of sig and w0 when no convection
635
636 ! Caution1: if no convection, we set iflag=4
637 ! (it used to be 0 in convect3)
638
639 ! Caution2: at this stage, tvp (and thus buoy) are know up
640 ! through icb only!
641 ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
642 ! -------------------------------------------------------------------
643
644 include "cv30param.h"
645
646 ! input:
647 INTEGER len, nd
648 INTEGER icb(len)
649 REAL plcl(len), p(len, nd)
650 REAL th(len, nd), tv(len, nd), tvp(len, nd)
651
652 ! output:
653 REAL pbase(len), buoybase(len)
654
655 ! input AND output:
656 INTEGER iflag(len)
657 REAL sig(len, nd), w0(len, nd)
658
659 ! local variables:
660 INTEGER i, k
661 REAL tvpbase, tvbase, tdif, ath, ath1
662
663
664 ! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
665
666 DO i = 1, len
667 pbase(i) = plcl(i) + dpbase
668 tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
669 (p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
670 p(i,icb(i))-p(i,icb(i)+1))
671 tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
672 (p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
673 (i,icb(i))-p(i,icb(i)+1))
674 buoybase(i) = tvpbase - tvbase
675 END DO
676
677
678 ! *** make sure that column is dry adiabatic between the surface ***
679 ! *** and cloud base, and that lifted air is positively buoyant ***
680 ! *** at cloud base ***
681 ! *** if not, return to calling program after resetting ***
682 ! *** sig(i) and w0(i) ***
683
684
685 ! oct3 do 200 i=1,len
686 ! oct3
687 ! oct3 tdif = buoybase(i)
688 ! oct3 ath1 = th(i,1)
689 ! oct3 ath = th(i,icb(i)-1) - dttrig
690 ! oct3
691 ! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
692 ! oct3 do 60 k=1,nl
693 ! oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
694 ! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
695 ! oct3 w0(i,k) = beta*w0(i,k)
696 ! oct3 60 continue
697 ! oct3 iflag(i)=4 ! pour version vectorisee
698 ! oct3c convect3 iflag(i)=0
699 ! oct3cccc return
700 ! oct3 endif
701 ! oct3
702 ! oct3200 continue
703
704 ! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
705
706 DO k = 1, nl
707 DO i = 1, len
708
709 tdif = buoybase(i)
710 ath1 = th(i, 1)
711 ath = th(i, icb(i)-1) - dttrig
712
713 IF (tdif<dtcrit .OR. ath>ath1) THEN
714 sig(i, k) = beta*sig(i, k) - 2.*alpha*tdif*tdif
715 sig(i, k) = amax1(sig(i,k), 0.0)
716 w0(i, k) = beta*w0(i, k)
717 iflag(i) = 4 ! pour version vectorisee
718 ! convect3 iflag(i)=0
719 END IF
720
721 END DO
722 END DO
723
724 ! fin oct3 --
725
726 RETURN
727 END SUBROUTINE cv30_trigger
728
729 SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
730 plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
731 th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
732 iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
733 v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
734 USE print_control_mod, ONLY: lunout
735 IMPLICIT NONE
736
737 include "cv30param.h"
738
739 ! inputs:
740 INTEGER len, ncum, nd, ntra, nloc
741 INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
742 REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
743 REAL pbase1(len), buoybase1(len)
744 REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
745 REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
746 REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
747 REAL tvp1(len, nd), clw1(len, nd)
748 REAL th1(len, nd)
749 REAL sig1(len, nd), w01(len, nd)
750 REAL tra1(len, nd, ntra)
751
752 ! outputs:
753 ! en fait, on a nloc=len pour l'instant (cf cv_driver)
754 INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
755 REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
756 REAL pbase(nloc), buoybase(nloc)
757 REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
758 REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
759 REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
760 REAL tvp(nloc, nd), clw(nloc, nd)
761 REAL th(nloc, nd)
762 REAL sig(nloc, nd), w0(nloc, nd)
763 REAL tra(nloc, nd, ntra)
764
765 ! local variables:
766 INTEGER i, k, nn, j
767
768 CHARACTER (LEN=20) :: modname = 'cv30_compress'
769 CHARACTER (LEN=80) :: abort_message
770
771
772 DO k = 1, nl + 1
773 nn = 0
774 DO i = 1, len
775 IF (iflag1(i)==0) THEN
776 nn = nn + 1
777 sig(nn, k) = sig1(i, k)
778 w0(nn, k) = w01(i, k)
779 t(nn, k) = t1(i, k)
780 q(nn, k) = q1(i, k)
781 qs(nn, k) = qs1(i, k)
782 u(nn, k) = u1(i, k)
783 v(nn, k) = v1(i, k)
784 gz(nn, k) = gz1(i, k)
785 h(nn, k) = h1(i, k)
786 lv(nn, k) = lv1(i, k)
787 cpn(nn, k) = cpn1(i, k)
788 p(nn, k) = p1(i, k)
789 ph(nn, k) = ph1(i, k)
790 tv(nn, k) = tv1(i, k)
791 tp(nn, k) = tp1(i, k)
792 tvp(nn, k) = tvp1(i, k)
793 clw(nn, k) = clw1(i, k)
794 th(nn, k) = th1(i, k)
795 END IF
796 END DO
797 END DO
798
799 ! do 121 j=1,ntra
800 ! do 111 k=1,nd
801 ! nn=0
802 ! do 101 i=1,len
803 ! if(iflag1(i).eq.0)then
804 ! nn=nn+1
805 ! tra(nn,k,j)=tra1(i,k,j)
806 ! endif
807 ! 101 continue
808 ! 111 continue
809 ! 121 continue
810
811 IF (nn/=ncum) THEN
812 WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
813 abort_message = ''
814 CALL abort_physic(modname, abort_message, 1)
815 END IF
816
817 nn = 0
818 DO i = 1, len
819 IF (iflag1(i)==0) THEN
820 nn = nn + 1
821 pbase(nn) = pbase1(i)
822 buoybase(nn) = buoybase1(i)
823 plcl(nn) = plcl1(i)
824 tnk(nn) = tnk1(i)
825 qnk(nn) = qnk1(i)
826 gznk(nn) = gznk1(i)
827 nk(nn) = nk1(i)
828 icb(nn) = icb1(i)
829 icbs(nn) = icbs1(i)
830 iflag(nn) = iflag1(i)
831 END IF
832 END DO
833
834 RETURN
835 END SUBROUTINE cv30_compress
836
837 SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
838 q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
839 ep, sigp, buoy)
840 ! epmax_cape: ajout arguments
841 IMPLICIT NONE
842
843 ! ---------------------------------------------------------------------
844 ! Purpose:
845 ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
846 ! &
847 ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
848 ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
849 ! &
850 ! FIND THE LEVEL OF NEUTRAL BUOYANCY
851
852 ! Main differences convect3/convect4:
853 ! - icbs (input) is the first level above LCL (may differ from icb)
854 ! - many minor differences in the iterations
855 ! - condensed water not removed from tvp in convect3
856 ! - vertical profile of buoyancy computed here (use of buoybase)
857 ! - the determination of inb is different
858 ! - no inb1, only inb in output
859 ! ---------------------------------------------------------------------
860
861 include "cvthermo.h"
862 include "cv30param.h"
863 include "conema3.h"
864
865 ! inputs:
866 INTEGER ncum, nd, nloc
867 INTEGER icb(nloc), icbs(nloc), nk(nloc)
868 REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
869 REAL p(nloc, nd)
870 REAL tnk(nloc), qnk(nloc), gznk(nloc)
871 REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
872 REAL pbase(nloc), buoybase(nloc), plcl(nloc)
873
874 ! outputs:
875 INTEGER inb(nloc)
876 REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
877 REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
878 REAL buoy(nloc, nd)
879
880 ! local variables:
881 INTEGER i, k
882 REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
883 REAL by, defrac, pden
884 REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
885 LOGICAL lcape(nloc)
886
887 ! =====================================================================
888 ! --- SOME INITIALIZATIONS
889 ! =====================================================================
890
891 DO k = 1, nl
892 DO i = 1, ncum
893 ep(i, k) = 0.0
894 sigp(i, k) = spfac
895 END DO
896 END DO
897
898 ! =====================================================================
899 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
900 ! =====================================================================
901
902 ! --- The procedure is to solve the equation.
903 ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
904
905 ! *** Calculate certain parcel quantities, including static energy ***
906
907
908 DO i = 1, ncum
909 ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) & ! debug &
910 ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
911 +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
912 END DO
913
914
915 ! *** Find lifted parcel quantities above cloud base ***
916
917
918 DO k = minorig + 1, nl
919 DO i = 1, ncum
920 ! ori if(k.ge.(icb(i)+1))then
921 IF (k>=(icbs(i)+1)) THEN ! convect3
922 tg = t(i, k)
923 qg = qs(i, k)
924 ! debug alv=lv0-clmcpv*(t(i,k)-t0)
925 alv = lv0 - clmcpv*(t(i,k)-273.15)
926
927 ! First iteration.
928
929 ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
930 s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
931 +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
932 s = 1./s
933 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
934 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
935 tg = tg + s*(ah0(i)-ahg)
936 ! ori tg=max(tg,35.0)
937 ! debug tc=tg-t0
938 tc = tg - 273.15
939 denom = 243.5 + tc
940 denom = max(denom, 1.0) ! convect3
941 ! ori if(tc.ge.0.0)then
942 es = 6.112*exp(17.67*tc/denom)
943 ! ori else
944 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
945 ! ori endif
946 qg = eps*es/(p(i,k)-es*(1.-eps))
947
948 ! Second iteration.
949
950 ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
951 ! ori s=1./s
952 ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
953 ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
954 tg = tg + s*(ah0(i)-ahg)
955 ! ori tg=max(tg,35.0)
956 ! debug tc=tg-t0
957 tc = tg - 273.15
958 denom = 243.5 + tc
959 denom = max(denom, 1.0) ! convect3
960 ! ori if(tc.ge.0.0)then
961 es = 6.112*exp(17.67*tc/denom)
962 ! ori else
963 ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
964 ! ori endif
965 qg = eps*es/(p(i,k)-es*(1.-eps))
966
967 ! debug alv=lv0-clmcpv*(t(i,k)-t0)
968 alv = lv0 - clmcpv*(t(i,k)-273.15)
969 ! print*,'cpd dans convect2 ',cpd
970 ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
971 ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
972
973 ! ori c approximation here:
974 ! ori
975 ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
976
977 ! convect3: no approximation:
978 tp(i, k) = (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
979
980 clw(i, k) = qnk(i) - qg
981 clw(i, k) = max(0.0, clw(i,k))
982 rg = qg/(1.-qnk(i))
983 ! ori tvp(i,k)=tp(i,k)*(1.+rg*epsi)
984 ! convect3: (qg utilise au lieu du vrai mixing ratio rg):
985 tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
986 END IF
987 END DO
988 END DO
989
990 ! =====================================================================
991 ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
992 ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
993 ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
994 ! =====================================================================
995
996 ! ori do 320 k=minorig+1,nl
997 DO k = 1, nl ! convect3
998 DO i = 1, ncum
999 pden = ptcrit - pbcrit
1000 ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1001 ep(i, k) = amax1(ep(i,k), 0.0)
1002 ep(i, k) = amin1(ep(i,k), epmax)
1003 sigp(i, k) = spfac
1004 ! ori if(k.ge.(nk(i)+1))then
1005 ! ori tca=tp(i,k)-t0
1006 ! ori if(tca.ge.0.0)then
1007 ! ori elacrit=elcrit
1008 ! ori else
1009 ! ori elacrit=elcrit*(1.0-tca/tlcrit)
1010 ! ori endif
1011 ! ori elacrit=max(elacrit,0.0)
1012 ! ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
1013 ! ori ep(i,k)=max(ep(i,k),0.0 )
1014 ! ori ep(i,k)=min(ep(i,k),1.0 )
1015 ! ori sigp(i,k)=sigs
1016 ! ori endif
1017 END DO
1018 END DO
1019
1020 ! =====================================================================
1021 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1022 ! --- VIRTUAL TEMPERATURE
1023 ! =====================================================================
1024
1025 ! dans convect3, tvp est calcule en une seule fois, et sans retirer
1026 ! l'eau condensee (~> reversible CAPE)
1027
1028 ! ori do 340 k=minorig+1,nl
1029 ! ori do 330 i=1,ncum
1030 ! ori if(k.ge.(icb(i)+1))then
1031 ! ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
1032 ! oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1033 ! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1034 ! ori endif
1035 ! ori 330 continue
1036 ! ori 340 continue
1037
1038 ! ori do 350 i=1,ncum
1039 ! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1040 ! ori 350 continue
1041
1042 DO i = 1, ncum ! convect3
1043 tp(i, nlp) = tp(i, nl) ! convect3
1044 END DO ! convect3
1045
1046 ! =====================================================================
1047 ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1048 ! =====================================================================
1049
1050 ! -- this is for convect3 only:
1051
1052 ! first estimate of buoyancy:
1053
1054 DO i = 1, ncum
1055 DO k = 1, nl
1056 buoy(i, k) = tvp(i, k) - tv(i, k)
1057 END DO
1058 END DO
1059
1060 ! set buoyancy=buoybase for all levels below base
1061 ! for safety, set buoy(icb)=buoybase
1062
1063 DO i = 1, ncum
1064 DO k = 1, nl
1065 IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1066 buoy(i, k) = buoybase(i)
1067 END IF
1068 END DO
1069 ! IM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
1070 buoy(i, icb(i)) = buoybase(i)
1071 END DO
1072
1073 ! -- end convect3
1074
1075 ! =====================================================================
1076 ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1077 ! --- LEVEL OF NEUTRAL BUOYANCY
1078 ! =====================================================================
1079
1080 ! -- this is for convect3 only:
1081
1082 DO i = 1, ncum
1083 inb(i) = nl - 1
1084 END DO
1085
1086 DO i = 1, ncum
1087 DO k = 1, nl - 1
1088 IF ((k>=icb(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1089 inb(i) = min(inb(i), k)
1090 END IF
1091 END DO
1092 END DO
1093
1094 ! -- end convect3
1095
1096 ! ori do 510 i=1,ncum
1097 ! ori cape(i)=0.0
1098 ! ori capem(i)=0.0
1099 ! ori inb(i)=icb(i)+1
1100 ! ori inb1(i)=inb(i)
1101 ! ori 510 continue
1102
1103 ! Originial Code
1104
1105 ! do 530 k=minorig+1,nl-1
1106 ! do 520 i=1,ncum
1107 ! if(k.ge.(icb(i)+1))then
1108 ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1109 ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1110 ! cape(i)=cape(i)+by
1111 ! if(by.ge.0.0)inb1(i)=k+1
1112 ! if(cape(i).gt.0.0)then
1113 ! inb(i)=k+1
1114 ! capem(i)=cape(i)
1115 ! endif
1116 ! endif
1117 ! 520 continue
1118 ! 530 continue
1119 ! do 540 i=1,ncum
1120 ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1121 ! cape(i)=capem(i)+byp
1122 ! defrac=capem(i)-cape(i)
1123 ! defrac=max(defrac,0.001)
1124 ! frac(i)=-cape(i)/defrac
1125 ! frac(i)=min(frac(i),1.0)
1126 ! frac(i)=max(frac(i),0.0)
1127 ! 540 continue
1128
1129 ! K Emanuel fix
1130
1131 ! call zilch(byp,ncum)
1132 ! do 530 k=minorig+1,nl-1
1133 ! do 520 i=1,ncum
1134 ! if(k.ge.(icb(i)+1))then
1135 ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1136 ! cape(i)=cape(i)+by
1137 ! if(by.ge.0.0)inb1(i)=k+1
1138 ! if(cape(i).gt.0.0)then
1139 ! inb(i)=k+1
1140 ! capem(i)=cape(i)
1141 ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1142 ! endif
1143 ! endif
1144 ! 520 continue
1145 ! 530 continue
1146 ! do 540 i=1,ncum
1147 ! inb(i)=max(inb(i),inb1(i))
1148 ! cape(i)=capem(i)+byp(i)
1149 ! defrac=capem(i)-cape(i)
1150 ! defrac=max(defrac,0.001)
1151 ! frac(i)=-cape(i)/defrac
1152 ! frac(i)=min(frac(i),1.0)
1153 ! frac(i)=max(frac(i),0.0)
1154 ! 540 continue
1155
1156 ! J Teixeira fix
1157
1158 ! ori call zilch(byp,ncum)
1159 ! ori do 515 i=1,ncum
1160 ! ori lcape(i)=.true.
1161 ! ori 515 continue
1162 ! ori do 530 k=minorig+1,nl-1
1163 ! ori do 520 i=1,ncum
1164 ! ori if(cape(i).lt.0.0)lcape(i)=.false.
1165 ! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1166 ! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1167 ! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1168 ! ori cape(i)=cape(i)+by
1169 ! ori if(by.ge.0.0)inb1(i)=k+1
1170 ! ori if(cape(i).gt.0.0)then
1171 ! ori inb(i)=k+1
1172 ! ori capem(i)=cape(i)
1173 ! ori endif
1174 ! ori endif
1175 ! ori 520 continue
1176 ! ori 530 continue
1177 ! ori do 540 i=1,ncum
1178 ! ori cape(i)=capem(i)+byp(i)
1179 ! ori defrac=capem(i)-cape(i)
1180 ! ori defrac=max(defrac,0.001)
1181 ! ori frac(i)=-cape(i)/defrac
1182 ! ori frac(i)=min(frac(i),1.0)
1183 ! ori frac(i)=max(frac(i),0.0)
1184 ! ori 540 continue
1185
1186 ! =====================================================================
1187 ! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1188 ! =====================================================================
1189
1190 ! ym do i=1,ncum*nlp
1191 ! ym hp(i,1)=h(i,1)
1192 ! ym enddo
1193
1194 DO k = 1, nlp
1195 DO i = 1, ncum
1196 hp(i, k) = h(i, k)
1197 END DO
1198 END DO
1199
1200 DO k = minorig + 1, nl
1201 DO i = 1, ncum
1202 IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1203 hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
1204 )
1205 END IF
1206 END DO
1207 END DO
1208
1209 RETURN
1210 END SUBROUTINE cv30_undilute2
1211
1212 SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
1213 sig, w0, cape, m)
1214 IMPLICIT NONE
1215
1216 ! ===================================================================
1217 ! --- CLOSURE OF CONVECT3
1218
1219 ! vectorization: S. Bony
1220 ! ===================================================================
1221
1222 include "cvthermo.h"
1223 include "cv30param.h"
1224
1225 ! input:
1226 INTEGER ncum, nd, nloc
1227 INTEGER icb(nloc), inb(nloc)
1228 REAL pbase(nloc)
1229 REAL p(nloc, nd), ph(nloc, nd+1)
1230 REAL tv(nloc, nd), buoy(nloc, nd)
1231
1232 ! input/output:
1233 REAL sig(nloc, nd), w0(nloc, nd)
1234
1235 ! output:
1236 REAL cape(nloc)
1237 REAL m(nloc, nd)
1238
1239 ! local variables:
1240 INTEGER i, j, k, icbmax
1241 REAL deltap, fac, w, amu
1242 REAL dtmin(nloc, nd), sigold(nloc, nd)
1243
1244 ! -------------------------------------------------------
1245 ! -- Initialization
1246 ! -------------------------------------------------------
1247
1248 DO k = 1, nl
1249 DO i = 1, ncum
1250 m(i, k) = 0.0
1251 END DO
1252 END DO
1253
1254 ! -------------------------------------------------------
1255 ! -- Reset sig(i) and w0(i) for i>inb and i<icb
1256 ! -------------------------------------------------------
1257
1258 ! update sig and w0 above LNB:
1259
1260 DO k = 1, nl - 1
1261 DO i = 1, ncum
1262 IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1263 sig(i, k) = beta*sig(i, k) + 2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb( &
1264 i)))
1265 sig(i, k) = amax1(sig(i,k), 0.0)
1266 w0(i, k) = beta*w0(i, k)
1267 END IF
1268 END DO
1269 END DO
1270
1271 ! compute icbmax:
1272
1273 icbmax = 2
1274 DO i = 1, ncum
1275 icbmax = max(icbmax, icb(i))
1276 END DO
1277
1278 ! update sig and w0 below cloud base:
1279
1280 DO k = 1, icbmax
1281 DO i = 1, ncum
1282 IF (k<=icb(i)) THEN
1283 sig(i, k) = beta*sig(i, k) - 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
1284 sig(i, k) = amax1(sig(i,k), 0.0)
1285 w0(i, k) = beta*w0(i, k)
1286 END IF
1287 END DO
1288 END DO
1289
1290 ! ! if(inb.lt.(nl-1))then
1291 ! ! do 85 i=inb+1,nl-1
1292 ! ! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
1293 ! ! 1 abs(buoy(inb))
1294 ! ! sig(i)=amax1(sig(i),0.0)
1295 ! ! w0(i)=beta*w0(i)
1296 ! ! 85 continue
1297 ! ! end if
1298
1299 ! ! do 87 i=1,icb
1300 ! ! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
1301 ! ! sig(i)=amax1(sig(i),0.0)
1302 ! ! w0(i)=beta*w0(i)
1303 ! ! 87 continue
1304
1305 ! -------------------------------------------------------------
1306 ! -- Reset fractional areas of updrafts and w0 at initial time
1307 ! -- and after 10 time steps of no convection
1308 ! -------------------------------------------------------------
1309
1310 DO k = 1, nl - 1
1311 DO i = 1, ncum
1312 IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1313 sig(i, k) = 0.0
1314 w0(i, k) = 0.0
1315 END IF
1316 END DO
1317 END DO
1318
1319 ! -------------------------------------------------------------
1320 ! -- Calculate convective available potential energy (cape),
1321 ! -- vertical velocity (w), fractional area covered by
1322 ! -- undilute updraft (sig), and updraft mass flux (m)
1323 ! -------------------------------------------------------------
1324
1325 DO i = 1, ncum
1326 cape(i) = 0.0
1327 END DO
1328
1329 ! compute dtmin (minimum buoyancy between ICB and given level k):
1330
1331 DO i = 1, ncum
1332 DO k = 1, nl
1333 dtmin(i, k) = 100.0
1334 END DO
1335 END DO
1336
1337 DO i = 1, ncum
1338 DO k = 1, nl
1339 DO j = minorig, nl
1340 IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
1341 1))) THEN
1342 dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1343 END IF
1344 END DO
1345 END DO
1346 END DO
1347
1348 ! the interval on which cape is computed starts at pbase :
1349 DO k = 1, nl
1350 DO i = 1, ncum
1351
1352 IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1353
1354 deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1355 cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1)
1356 cape(i) = amax1(0.0, cape(i))
1357 sigold(i, k) = sig(i, k)
1358
1359 ! dtmin(i,k)=100.0
1360 ! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1361 ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1362 ! 97 continue
1363
1364 sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
1365 sig(i, k) = amax1(sig(i,k), 0.0)
1366 sig(i, k) = amin1(sig(i,k), 0.01)
1367 fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1368 w = (1.-beta)*fac*sqrt(cape(i)) + beta*w0(i, k)
1369 amu = 0.5*(sig(i,k)+sigold(i,k))*w
1370 m(i, k) = amu*0.007*p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1371 w0(i, k) = w
1372 END IF
1373
1374 END DO
1375 END DO
1376
1377 DO i = 1, ncum
1378 w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1379 m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/ &
1380 (ph(i,icb(i)+1)-ph(i,icb(i)+2))
1381 sig(i, icb(i)) = sig(i, icb(i)+1)
1382 sig(i, icb(i)-1) = sig(i, icb(i))
1383 END DO
1384
1385
1386 ! ! cape=0.0
1387 ! ! do 98 i=icb+1,inb
1388 ! ! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1389 ! ! cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
1390 ! ! dcape=rrd*buoy(i-1)*deltap/p(i-1)
1391 ! ! dlnp=deltap/p(i-1)
1392 ! ! cape=amax1(0.0,cape)
1393 ! ! sigold=sig(i)
1394
1395 ! ! dtmin=100.0
1396 ! ! do 97 j=icb,i-1
1397 ! ! dtmin=amin1(dtmin,buoy(j))
1398 ! ! 97 continue
1399
1400 ! ! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
1401 ! ! sig(i)=amax1(sig(i),0.0)
1402 ! ! sig(i)=amin1(sig(i),0.01)
1403 ! ! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1404 ! ! w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
1405 ! ! amu=0.5*(sig(i)+sigold)*w
1406 ! ! m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
1407 ! ! w0(i)=w
1408 ! ! 98 continue
1409 ! ! w0(icb)=0.5*w0(icb+1)
1410 ! ! m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1411 ! ! sig(icb)=sig(icb+1)
1412 ! ! sig(icb-1)=sig(icb)
1413
1414 RETURN
1415 END SUBROUTINE cv30_closure
1416
1417 SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
1418 u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
1419 vent, sij, elij, ments, qents, traent)
1420 IMPLICIT NONE
1421
1422 ! ---------------------------------------------------------------------
1423 ! a faire:
1424 ! - changer rr(il,1) -> qnk(il)
1425 ! - vectorisation de la partie normalisation des flux (do 789...)
1426 ! ---------------------------------------------------------------------
1427
1428 include "cvthermo.h"
1429 include "cv30param.h"
1430
1431 ! inputs:
1432 INTEGER ncum, nd, na, ntra, nloc
1433 INTEGER icb(nloc), inb(nloc), nk(nloc)
1434 REAL sig(nloc, nd)
1435 REAL qnk(nloc)
1436 REAL ph(nloc, nd+1)
1437 REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1438 REAL u(nloc, nd), v(nloc, nd)
1439 REAL tra(nloc, nd, ntra) ! input of convect3
1440 REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
1441 REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
1442 REAL m(nloc, na) ! input of convect3
1443
1444 ! outputs:
1445 REAL ment(nloc, na, na), qent(nloc, na, na)
1446 REAL uent(nloc, na, na), vent(nloc, na, na)
1447 REAL sij(nloc, na, na), elij(nloc, na, na)
1448 REAL traent(nloc, nd, nd, ntra)
1449 REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
1450 REAL sigij(nloc, nd, nd)
1451
1452 ! local variables:
1453 INTEGER i, j, k, il, im, jm
1454 INTEGER num1, num2
1455 INTEGER nent(nloc, na)
1456 REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1457 REAL alt, smid, sjmin, sjmax, delp, delm
1458 REAL asij(nloc), smax(nloc), scrit(nloc)
1459 REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1460 REAL wgh
1461 REAL zm(nloc, na)
1462 LOGICAL lwork(nloc)
1463
1464 ! =====================================================================
1465 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1466 ! =====================================================================
1467
1468 ! ori do 360 i=1,ncum*nlp
1469 DO j = 1, nl
1470 DO i = 1, ncum
1471 nent(i, j) = 0
1472 ! in convect3, m is computed in cv3_closure
1473 ! ori m(i,1)=0.0
1474 END DO
1475 END DO
1476
1477 ! ori do 400 k=1,nlp
1478 ! ori do 390 j=1,nlp
1479 DO j = 1, nl
1480 DO k = 1, nl
1481 DO i = 1, ncum
1482 qent(i, k, j) = rr(i, j)
1483 uent(i, k, j) = u(i, j)
1484 vent(i, k, j) = v(i, j)
1485 elij(i, k, j) = 0.0
1486 ! ym ment(i,k,j)=0.0
1487 ! ym sij(i,k,j)=0.0
1488 END DO
1489 END DO
1490 END DO
1491
1492 ! ym
1493 ment(1:ncum, 1:nd, 1:nd) = 0.0
1494 sij(1:ncum, 1:nd, 1:nd) = 0.0
1495
1496 ! do k=1,ntra
1497 ! do j=1,nd ! instead nlp
1498 ! do i=1,nd ! instead nlp
1499 ! do il=1,ncum
1500 ! traent(il,i,j,k)=tra(il,j,k)
1501 ! enddo
1502 ! enddo
1503 ! enddo
1504 ! enddo
1505 zm(:, :) = 0.
1506
1507 ! =====================================================================
1508 ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1509 ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1510 ! --- FRACTION (sij)
1511 ! =====================================================================
1512
1513 DO i = minorig + 1, nl
1514
1515 DO j = minorig, nl
1516 DO il = 1, ncum
1517 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
1518 1)) .AND. (j<=inb(il))) THEN
1519
1520 rti = rr(il, 1) - ep(il, i)*clw(il, i)
1521 bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
1522 anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
1523 denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
1524 dei = denom
1525 IF (abs(dei)<0.01) dei = 0.01
1526 sij(il, i, j) = anum/dei
1527 sij(il, i, i) = 1.0
1528 altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
1529 altem = altem/bf2
1530 cwat = clw(il, j)*(1.-ep(il,j))
1531 stemp = sij(il, i, j)
1532 IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1533 anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
1534 denom = denom + lv(il, j)*(rr(il,i)-rti)
1535 IF (abs(denom)<0.01) denom = 0.01
1536 sij(il, i, j) = anum/denom
1537 altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - &
1538 rs(il, j)
1539 altem = altem - (bf2-1.)*cwat
1540 END IF
1541 IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1542 qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti
1543 uent(il, i, j) = sij(il, i, j)*u(il, i) + &
1544 (1.-sij(il,i,j))*u(il, nk(il))
1545 vent(il, i, j) = sij(il, i, j)*v(il, i) + &
1546 (1.-sij(il,i,j))*v(il, nk(il))
1547 ! !!! do k=1,ntra
1548 ! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1549 ! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1550 ! !!! end do
1551 elij(il, i, j) = altem
1552 elij(il, i, j) = amax1(0.0, elij(il,i,j))
1553 ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
1554 nent(il, i) = nent(il, i) + 1
1555 END IF
1556 sij(il, i, j) = amax1(0.0, sij(il,i,j))
1557 sij(il, i, j) = amin1(1.0, sij(il,i,j))
1558 END IF ! new
1559 END DO
1560 END DO
1561
1562 ! do k=1,ntra
1563 ! do j=minorig,nl
1564 ! do il=1,ncum
1565 ! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1566 ! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1567 ! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1568 ! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1569 ! endif
1570 ! enddo
1571 ! enddo
1572 ! enddo
1573
1574
1575 ! *** if no air can entrain at level i assume that updraft detrains
1576 ! ***
1577 ! *** at that level and calculate detrained air flux and properties
1578 ! ***
1579
1580
1581 ! @ do 170 i=icb(il),inb(il)
1582
1583 DO il = 1, ncum
1584 IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
1585 ! @ if(nent(il,i).eq.0)then
1586 ment(il, i, i) = m(il, i)
1587 qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i)
1588 uent(il, i, i) = u(il, nk(il))
1589 vent(il, i, i) = v(il, nk(il))
1590 elij(il, i, i) = clw(il, i)
1591 ! MAF sij(il,i,i)=1.0
1592 sij(il, i, i) = 0.0
1593 END IF
1594 END DO
1595 END DO
1596
1597 ! do j=1,ntra
1598 ! do i=minorig+1,nl
1599 ! do il=1,ncum
1600 ! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1601 ! traent(il,i,i,j)=tra(il,nk(il),j)
1602 ! endif
1603 ! enddo
1604 ! enddo
1605 ! enddo
1606
1607 DO j = minorig, nl
1608 DO i = minorig, nl
1609 DO il = 1, ncum
1610 IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
1611 inb(il))) THEN
1612 sigij(il, i, j) = sij(il, i, j)
1613 END IF
1614 END DO
1615 END DO
1616 END DO
1617 ! @ enddo
1618
1619 ! @170 continue
1620
1621 ! =====================================================================
1622 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
1623 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1624 ! =====================================================================
1625
1626 ! ym call zilch(asum,ncum*nd)
1627 ! ym call zilch(bsum,ncum*nd)
1628 ! ym call zilch(csum,ncum*nd)
1629 CALL zilch(asum, nloc*nd)
1630 CALL zilch(csum, nloc*nd)
1631 CALL zilch(csum, nloc*nd)
1632
1633 DO il = 1, ncum
1634 lwork(il) = .FALSE.
1635 END DO
1636
1637 DO i = minorig + 1, nl
1638
1639 num1 = 0
1640 DO il = 1, ncum
1641 IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
1642 END DO
1643 IF (num1<=0) GO TO 789
1644
1645
1646 DO il = 1, ncum
1647 IF (i>=icb(il) .AND. i<=inb(il)) THEN
1648 lwork(il) = (nent(il,i)/=0)
1649 qp = rr(il, 1) - ep(il, i)*clw(il, i)
1650 anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
1651 (cpv-cpd)*t(il, i)*(qp-rr(il,i))
1652 denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
1653 (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
1654 IF (abs(denom)<0.01) denom = 0.01
1655 scrit(il) = anum/denom
1656 alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
1657 IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
1658 smax(il) = 0.0
1659 asij(il) = 0.0
1660 END IF
1661 END DO
1662
1663 DO j = nl, minorig, -1
1664
1665 num2 = 0
1666 DO il = 1, ncum
1667 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1668 il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
1669 END DO
1670 IF (num2<=0) GO TO 175
1671
1672 DO il = 1, ncum
1673 IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1674 il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
1675
1676 IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
1677 wgh = 1.0
1678 IF (j>i) THEN
1679 sjmax = amax1(sij(il,i,j+1), smax(il))
1680 sjmax = amin1(sjmax, scrit(il))
1681 smax(il) = amax1(sij(il,i,j), smax(il))
1682 sjmin = amax1(sij(il,i,j-1), smax(il))
1683 sjmin = amin1(sjmin, scrit(il))
1684 IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
1685 smid = amin1(sij(il,i,j), scrit(il))
1686 ELSE
1687 sjmax = amax1(sij(il,i,j+1), scrit(il))
1688 smid = amax1(sij(il,i,j), scrit(il))
1689 sjmin = 0.0
1690 IF (j>1) sjmin = sij(il, i, j-1)
1691 sjmin = amax1(sjmin, scrit(il))
1692 END IF
1693 delp = abs(sjmax-smid)
1694 delm = abs(sjmin-smid)
1695 asij(il) = asij(il) + wgh*(delp+delm)
1696 ment(il, i, j) = ment(il, i, j)*(delp+delm)*wgh
1697 END IF
1698 END IF
1699 END DO
1700
1701 175 END DO
1702
1703 DO il = 1, ncum
1704 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1705 asij(il) = amax1(1.0E-16, asij(il))
1706 asij(il) = 1.0/asij(il)
1707 asum(il, i) = 0.0
1708 bsum(il, i) = 0.0
1709 csum(il, i) = 0.0
1710 END IF
1711 END DO
1712
1713 DO j = minorig, nl
1714 DO il = 1, ncum
1715 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1716 il)-1) .AND. j<=inb(il)) THEN
1717 ment(il, i, j) = ment(il, i, j)*asij(il)
1718 END IF
1719 END DO
1720 END DO
1721
1722 DO j = minorig, nl
1723 DO il = 1, ncum
1724 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1725 il)-1) .AND. j<=inb(il)) THEN
1726 asum(il, i) = asum(il, i) + ment(il, i, j)
1727 ment(il, i, j) = ment(il, i, j)*sig(il, j)
1728 bsum(il, i) = bsum(il, i) + ment(il, i, j)
1729 END IF
1730 END DO
1731 END DO
1732
1733 DO il = 1, ncum
1734 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1735 bsum(il, i) = amax1(bsum(il,i), 1.0E-16)
1736 bsum(il, i) = 1.0/bsum(il, i)
1737 END IF
1738 END DO
1739
1740 DO j = minorig, nl
1741 DO il = 1, ncum
1742 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1743 il)-1) .AND. j<=inb(il)) THEN
1744 ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i)
1745 END IF
1746 END DO
1747 END DO
1748
1749 DO j = minorig, nl
1750 DO il = 1, ncum
1751 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1752 il)-1) .AND. j<=inb(il)) THEN
1753 csum(il, i) = csum(il, i) + ment(il, i, j)
1754 END IF
1755 END DO
1756 END DO
1757
1758 DO il = 1, ncum
1759 IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1760 csum(il,i)<m(il,i)) THEN
1761 nent(il, i) = 0
1762 ment(il, i, i) = m(il, i)
1763 qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i)
1764 uent(il, i, i) = u(il, nk(il))
1765 vent(il, i, i) = v(il, nk(il))
1766 elij(il, i, i) = clw(il, i)
1767 ! MAF sij(il,i,i)=1.0
1768 sij(il, i, i) = 0.0
1769 END IF
1770 END DO ! il
1771
1772 ! do j=1,ntra
1773 ! do il=1,ncum
1774 ! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1775 ! : .and. csum(il,i).lt.m(il,i) ) then
1776 ! traent(il,i,i,j)=tra(il,nk(il),j)
1777 ! endif
1778 ! enddo
1779 ! enddo
1780 789 END DO
1781
1782 ! MAF: renormalisation de MENT
1783 DO jm = 1, nd
1784 DO im = 1, nd
1785 DO il = 1, ncum
1786 zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
1787 END DO
1788 END DO
1789 END DO
1790
1791 DO jm = 1, nd
1792 DO im = 1, nd
1793 DO il = 1, ncum
1794 IF (zm(il,im)/=0.) THEN
1795 ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
1796 END IF
1797 END DO
1798 END DO
1799 END DO
1800
1801 DO jm = 1, nd
1802 DO im = 1, nd
1803 DO il = 1, ncum
1804 qents(il, im, jm) = qent(il, im, jm)
1805 ments(il, im, jm) = ment(il, im, jm)
1806 END DO
1807 END DO
1808 END DO
1809
1810 RETURN
1811 END SUBROUTINE cv30_mixing
1812
1813
1814 SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
1815 v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
1816 mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
1817 , wdtraina, wdtrainm) ! 26/08/10 RomP-jyg
1818 IMPLICIT NONE
1819
1820
1821 include "cvthermo.h"
1822 include "cv30param.h"
1823 include "cvflag.h"
1824
1825 ! inputs:
1826 INTEGER ncum, nd, na, ntra, nloc
1827 INTEGER icb(nloc), inb(nloc)
1828 REAL delt, plcl(nloc)
1829 REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1830 REAL u(nloc, nd), v(nloc, nd)
1831 REAL tra(nloc, nd, ntra)
1832 REAL p(nloc, nd), ph(nloc, nd+1)
1833 REAL th(nloc, na), gz(nloc, na)
1834 REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
1835 REAL cpn(nloc, na), tv(nloc, na)
1836 REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
1837
1838 ! outputs:
1839 REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
1840 REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
1841 REAL trap(nloc, na, ntra)
1842 REAL b(nloc, na)
1843 ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
1844 ! lascendance adiabatique et des flux melanges Pa et Pm.
1845 ! Distinction des wdtrain
1846 ! Pa = wdtrainA Pm = wdtrainM
1847 REAL wdtraina(nloc, na), wdtrainm(nloc, na)
1848
1849 ! local variables
1850 INTEGER i, j, k, il, num1
1851 REAL tinv, delti
1852 REAL awat, afac, afac1, afac2, bfac
1853 REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
1854 REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1855 REAL ampmax
1856 REAL lvcp(nloc, na)
1857 REAL wdtrain(nloc)
1858 LOGICAL lwork(nloc)
1859
1860
1861 ! ------------------------------------------------------
1862
1863 delti = 1./delt
1864 tinv = 1./3.
1865
1866 mp(:, :) = 0.
1867
1868 DO i = 1, nl
1869 DO il = 1, ncum
1870 mp(il, i) = 0.0
1871 rp(il, i) = rr(il, i)
1872 up(il, i) = u(il, i)
1873 vp(il, i) = v(il, i)
1874 wt(il, i) = 0.001
1875 water(il, i) = 0.0
1876 evap(il, i) = 0.0
1877 b(il, i) = 0.0
1878 lvcp(il, i) = lv(il, i)/cpn(il, i)
1879 END DO
1880 END DO
1881
1882 ! do k=1,ntra
1883 ! do i=1,nd
1884 ! do il=1,ncum
1885 ! trap(il,i,k)=tra(il,i,k)
1886 ! enddo
1887 ! enddo
1888 ! enddo
1889 ! ! RomP >>>
1890 DO i = 1, nd
1891 DO il = 1, ncum
1892 wdtraina(il, i) = 0.0
1893 wdtrainm(il, i) = 0.0
1894 END DO
1895 END DO
1896 ! ! RomP <<<
1897
1898 ! *** check whether ep(inb)=0, if so, skip precipitating ***
1899 ! *** downdraft calculation ***
1900
1901
1902 DO il = 1, ncum
1903 lwork(il) = .TRUE.
1904 IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE.
1905 END DO
1906
1907 CALL zilch(wdtrain, ncum)
1908
1909 DO i = nl + 1, 1, -1
1910
1911 num1 = 0
1912 DO il = 1, ncum
1913 IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
1914 END DO
1915 IF (num1<=0) GO TO 400
1916
1917
1918 ! *** integrate liquid water equation to find condensed water ***
1919 ! *** and condensed water flux ***
1920
1921
1922
1923 ! *** begin downdraft loop ***
1924
1925
1926
1927 ! *** calculate detrained precipitation ***
1928
1929 DO il = 1, ncum
1930 IF (i<=inb(il) .AND. lwork(il)) THEN
1931 IF (cvflag_grav) THEN
1932 wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
1933 wdtraina(il, i) = wdtrain(il)/grav ! Pa 26/08/10 RomP
1934 ELSE
1935 wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
1936 wdtraina(il, i) = wdtrain(il)/10. ! Pa 26/08/10 RomP
1937 END IF
1938 END IF
1939 END DO
1940
1941 IF (i>1) THEN
1942
1943 DO j = 1, i - 1
1944 DO il = 1, ncum
1945 IF (i<=inb(il) .AND. lwork(il)) THEN
1946 awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
1947 awat = amax1(awat, 0.0)
1948 IF (cvflag_grav) THEN
1949 wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
1950 ELSE
1951 wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
1952 END IF
1953 END IF
1954 END DO
1955 END DO
1956 DO il = 1, ncum
1957 IF (cvflag_grav) THEN
1958 wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) ! Pm 26/08/10 RomP
1959 ELSE
1960 wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) ! Pm 26/08/10 RomP
1961 END IF
1962 END DO
1963
1964 END IF
1965
1966
1967 ! *** find rain water and evaporation using provisional ***
1968 ! *** estimates of rp(i)and rp(i-1) ***
1969
1970
1971 DO il = 1, ncum
1972
1973 IF (i<=inb(il) .AND. lwork(il)) THEN
1974
1975 wt(il, i) = 45.0
1976
1977 IF (i<inb(il)) THEN
1978 rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
1979 i))+gz(il,i+1)-gz(il,i))/lv(il, i)
1980 rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
1981 END IF
1982 rp(il, i) = amax1(rp(il,i), 0.0)
1983 rp(il, i) = amin1(rp(il,i), rs(il,i))
1984 rp(il, inb(il)) = rr(il, inb(il))
1985
1986 IF (i==1) THEN
1987 afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
1988 ELSE
1989 rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
1990 i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
1991 rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
1992 rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
1993 rp(il, i-1) = amax1(rp(il,i-1), 0.0)
1994 afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i) &
1995 )
1996 afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
1997 (1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
1998 afac = 0.5*(afac1+afac2)
1999 END IF
2000 IF (i==inb(il)) afac = 0.0
2001 afac = amax1(afac, 0.0)
2002 bfac = 1./(sigd*wt(il,i))
2003
2004 ! jyg1
2005 ! cc sigt=1.0
2006 ! cc if(i.ge.icb)sigt=sigp(i)
2007 ! prise en compte de la variation progressive de sigt dans
2008 ! les couches icb et icb-1:
2009 ! pour plcl<ph(i+1), pr1=0 & pr2=1
2010 ! pour plcl>ph(i), pr1=1 & pr2=0
2011 ! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2012 ! sur le nuage, et pr2 est la proportion sous la base du
2013 ! nuage.
2014 pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2015 pr1 = max(0., min(1.,pr1))
2016 pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2017 pr2 = max(0., min(1.,pr2))
2018 sigt = sigp(il, i)*pr1 + pr2
2019 ! jyg2
2020
2021 b6 = bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac
2022 c6 = water(il, i+1) + bfac*wdtrain(il) - 50.*sigd*bfac*(ph(il,i)-ph( &
2023 il,i+1))*evap(il, i+1)
2024 IF (c6>0.0) THEN
2025 revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
2026 evap(il, i) = sigt*afac*revap
2027 water(il, i) = revap*revap
2028 ELSE
2029 evap(il, i) = -evap(il, i+1) + 0.02*(wdtrain(il)+sigd*wt(il,i)* &
2030 water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1)))
2031 END IF
2032
2033 ! *** calculate precipitating downdraft mass flux under ***
2034 ! *** hydrostatic approximation ***
2035
2036 IF (i/=1) THEN
2037
2038 tevap = amax1(0.0, evap(il,i))
2039 delth = amax1(0.001, (th(il,i)-th(il,i-1)))
2040 IF (cvflag_grav) THEN
2041 mp(il, i) = 100.*ginv*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/ &
2042 delth
2043 ELSE
2044 mp(il, i) = 10.*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth
2045 END IF
2046
2047 ! *** if hydrostatic assumption fails, ***
2048 ! *** solve cubic difference equation for downdraft theta ***
2049 ! *** and mass flux from two simultaneous differential eqns ***
2050
2051 amfac = sigd*sigd*70.0*ph(il, i)*(p(il,i-1)-p(il,i))* &
2052 (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2053 amp2 = abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
2054 IF (amp2>(0.1*amfac)) THEN
2055 xf = 100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1))
2056 tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i)* &
2057 sigd*th(il,i))
2058 af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv
2059 bf = 2.*(tinv*mp(il,i+1))**3 + tinv*mp(il, i+1)*xf*tf + &
2060 50.*(p(il,i-1)-p(il,i))*xf*tevap
2061 fac2 = 1.0
2062 IF (bf<0.0) fac2 = -1.0
2063 bf = abs(bf)
2064 ur = 0.25*bf*bf - af*af*af*tinv*tinv*tinv
2065 IF (ur>=0.0) THEN
2066 sru = sqrt(ur)
2067 fac = 1.0
2068 IF ((0.5*bf-sru)<0.0) fac = -1.0
2069 mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
2070 fac*(abs(0.5*bf-sru))**tinv
2071 ELSE
2072 d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2073 IF (fac2<0.0) d = 3.14159 - d
2074 mp(il, i) = mp(il, i+1)*tinv + 2.*sqrt(af*tinv)*cos(d*tinv)
2075 END IF
2076 mp(il, i) = amax1(0.0, mp(il,i))
2077
2078 IF (cvflag_grav) THEN
2079 ! jyg : il y a vraisemblablement une erreur dans la ligne 2
2080 ! suivante:
2081 ! il faut diviser par (mp(il,i)*sigd*grav) et non par
2082 ! (mp(il,i)+sigd*0.1).
2083 ! Et il faut bien revoir les facteurs 100.
2084 b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, &
2085 i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2086 )*sigd*th(il,i))
2087 ELSE
2088 b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, &
2089 i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2090 )*sigd*th(il,i))
2091 END IF
2092 b(il, i-1) = amax1(b(il,i-1), 0.0)
2093 END IF
2094
2095 ! *** limit magnitude of mp(i) to meet cfl condition
2096 ! ***
2097
2098 ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
2099 amp2 = 2.0*(ph(il,i-1)-ph(il,i))*delti
2100 ampmax = amin1(ampmax, amp2)
2101 mp(il, i) = amin1(mp(il,i), ampmax)
2102
2103 ! *** force mp to decrease linearly to zero
2104 ! ***
2105 ! *** between cloud base and the surface
2106 ! ***
2107
2108 IF (p(il,i)>p(il,icb(il))) THEN
2109 mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ &
2110 (p(il,1)-p(il,icb(il)))
2111 END IF
2112
2113 END IF ! i.eq.1
2114
2115 ! *** find mixing ratio of precipitating downdraft ***
2116
2117
2118 IF (i/=inb(il)) THEN
2119
2120 rp(il, i) = rr(il, i)
2121
2122 IF (mp(il,i)>mp(il,i+1)) THEN
2123
2124 IF (cvflag_grav) THEN
2125 rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2126 rr(il, i)*(mp(il,i)-mp(il,i+1)) + 100.*ginv*0.5*sigd*(ph(il,i &
2127 )-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2128 ELSE
2129 rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2130 rr(il, i)*(mp(il,i)-mp(il,i+1)) + 5.*sigd*(ph(il,i)-ph(il,i+1 &
2131 ))*(evap(il,i+1)+evap(il,i))
2132 END IF
2133 rp(il, i) = rp(il, i)/mp(il, i)
2134 up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+ &
2135 1))
2136 up(il, i) = up(il, i)/mp(il, i)
2137 vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+ &
2138 1))
2139 vp(il, i) = vp(il, i)/mp(il, i)
2140
2141 ! do j=1,ntra
2142 ! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2143 ! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2144 ! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2145 ! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2146 ! end do
2147
2148 ELSE
2149
2150 IF (mp(il,i+1)>1.0E-16) THEN
2151 IF (cvflag_grav) THEN
2152 rp(il, i) = rp(il, i+1) + 100.*ginv*0.5*sigd*(ph(il,i)-ph(il, &
2153 i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2154 ELSE
2155 rp(il, i) = rp(il, i+1) + 5.*sigd*(ph(il,i)-ph(il,i+1))*(evap &
2156 (il,i+1)+evap(il,i))/mp(il, i+1)
2157 END IF
2158 up(il, i) = up(il, i+1)
2159 vp(il, i) = vp(il, i+1)
2160
2161 ! do j=1,ntra
2162 ! trap(il,i,j)=trap(il,i+1,j)
2163 ! end do
2164
2165 END IF
2166 END IF
2167 rp(il, i) = amin1(rp(il,i), rs(il,i))
2168 rp(il, i) = amax1(rp(il,i), 0.0)
2169
2170 END IF
2171 END IF
2172 END DO
2173
2174 400 END DO
2175
2176 RETURN
2177 END SUBROUTINE cv30_unsat
2178
2179 SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
2180 tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
2181 wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
2182 tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
2183 mike, tls, tps, qcondc, wd)
2184 IMPLICIT NONE
2185
2186 include "cvthermo.h"
2187 include "cv30param.h"
2188 include "cvflag.h"
2189 include "conema3.h"
2190
2191 ! inputs:
2192 INTEGER ncum, nd, na, ntra, nloc
2193 INTEGER icb(nloc), inb(nloc)
2194 REAL delt
2195 REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
2196 REAL tra(nloc, nd, ntra), sig(nloc, nd)
2197 REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
2198 REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
2199 REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
2200 REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
2201 REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
2202 REAL water(nloc, na), evap(nloc, na), b(nloc, na)
2203 REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
2204 ! ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2205 REAL vent(nloc, na, na), elij(nloc, na, na)
2206 INTEGER nent(nloc, na)
2207 REAL traent(nloc, na, na, ntra)
2208 REAL tv(nloc, nd), tvp(nloc, nd)
2209
2210 ! input/output:
2211 INTEGER iflag(nloc)
2212
2213 ! outputs:
2214 REAL precip(nloc)
2215 REAL vprecip(nloc, nd+1)
2216 REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
2217 REAL ftra(nloc, nd, ntra)
2218 REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
2219 REAL dnwd0(nloc, nd), mike(nloc, nd)
2220 REAL tls(nloc, nd), tps(nloc, nd)
2221 REAL qcondc(nloc, nd) ! cld
2222 REAL wd(nloc) ! gust
2223
2224 ! local variables:
2225 INTEGER i, k, il, n, j, num1
2226 REAL rat, awat, delti
2227 REAL ax, bx, cx, dx, ex
2228 REAL cpinv, rdcp, dpinv
2229 REAL lvcp(nloc, na), mke(nloc, na)
2230 REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
2231 ! !! real up1(nloc), dn1(nloc)
2232 REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
2233 REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2234 REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
2235 REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
2236
2237
2238 ! -------------------------------------------------------------
2239
2240 ! initialization:
2241
2242 delti = 1.0/delt
2243
2244 DO il = 1, ncum
2245 precip(il) = 0.0
2246 wd(il) = 0.0 ! gust
2247 vprecip(il, nd+1) = 0.
2248 END DO
2249
2250 DO i = 1, nd
2251 DO il = 1, ncum
2252 vprecip(il, i) = 0.0
2253 ft(il, i) = 0.0
2254 fr(il, i) = 0.0
2255 fu(il, i) = 0.0
2256 fv(il, i) = 0.0
2257 qcondc(il, i) = 0.0 ! cld
2258 qcond(il, i) = 0.0 ! cld
2259 nqcond(il, i) = 0.0 ! cld
2260 END DO
2261 END DO
2262
2263 ! do j=1,ntra
2264 ! do i=1,nd
2265 ! do il=1,ncum
2266 ! ftra(il,i,j)=0.0
2267 ! enddo
2268 ! enddo
2269 ! enddo
2270
2271 DO i = 1, nl
2272 DO il = 1, ncum
2273 lvcp(il, i) = lv(il, i)/cpn(il, i)
2274 END DO
2275 END DO
2276
2277
2278
2279 ! *** calculate surface precipitation in mm/day ***
2280
2281 DO il = 1, ncum
2282 IF (ep(il,inb(il))>=0.0001) THEN
2283 IF (cvflag_grav) THEN
2284 precip(il) = wt(il, 1)*sigd*water(il, 1)*86400.*1000./(rowl*grav)
2285 ELSE
2286 precip(il) = wt(il, 1)*sigd*water(il, 1)*8640.
2287 END IF
2288 END IF
2289 END DO
2290
2291 ! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
2292
2293 ! MAF rajout pour lessivage
2294 DO k = 1, nl
2295 DO il = 1, ncum
2296 IF (k<=inb(il)) THEN
2297 IF (cvflag_grav) THEN
2298 vprecip(il, k) = wt(il, k)*sigd*water(il, k)/grav
2299 ELSE
2300 vprecip(il, k) = wt(il, k)*sigd*water(il, k)/10.
2301 END IF
2302 END IF
2303 END DO
2304 END DO
2305
2306
2307 ! *** Calculate downdraft velocity scale ***
2308 ! *** NE PAS UTILISER POUR L'INSTANT ***
2309
2310 ! ! do il=1,ncum
2311 ! ! wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
2312 ! ! : /(sigd*p(il,icb(il)))
2313 ! ! enddo
2314
2315
2316 ! *** calculate tendencies of lowest level potential temperature ***
2317 ! *** and mixing ratio ***
2318
2319 DO il = 1, ncum
2320 work(il) = 1.0/(ph(il,1)-ph(il,2))
2321 am(il) = 0.0
2322 END DO
2323
2324 DO k = 2, nl
2325 DO il = 1, ncum
2326 IF (k<=inb(il)) THEN
2327 am(il) = am(il) + m(il, k)
2328 END IF
2329 END DO
2330 END DO
2331
2332 DO il = 1, ncum
2333
2334 ! convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4
2335 IF (cvflag_grav) THEN
2336 IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
2337 ft(il, 1) = 0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2338 1))/cpn(il,1))
2339 ELSE
2340 IF ((0.1*work(il)*am(il))>=delti) iflag(il) = 1 !consistency vect
2341 ft(il, 1) = 0.1*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2342 1))/cpn(il,1))
2343 END IF
2344
2345 ft(il, 1) = ft(il, 1) - 0.5*lvcp(il, 1)*sigd*(evap(il,1)+evap(il,2))
2346
2347 IF (cvflag_grav) THEN
2348 ft(il, 1) = ft(il, 1) - 0.009*grav*sigd*mp(il, 2)*t(il, 1)*b(il, 1)* &
2349 work(il)
2350 ELSE
2351 ft(il, 1) = ft(il, 1) - 0.09*sigd*mp(il, 2)*t(il, 1)*b(il, 1)*work(il)
2352 END IF
2353
2354 ft(il, 1) = ft(il, 1) + 0.01*sigd*wt(il, 1)*(cl-cpd)*water(il, 2)*(t(il,2 &
2355 )-t(il,1))*work(il)/cpn(il, 1)
2356
2357 IF (cvflag_grav) THEN
2358 ! jyg1 Correction pour mieux conserver l'eau (conformite avec
2359 ! CONVECT4.3)
2360 ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
2361 ! evap)
2362 fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + &
2363 sigd*0.5*(evap(il,1)+evap(il,2))
2364 ! +tard : +sigd*evap(il,1)
2365
2366 fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
2367
2368 fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il, &
2369 1))+am(il)*(u(il,2)-u(il,1)))
2370 fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
2371 1))+am(il)*(v(il,2)-v(il,1)))
2372 ELSE ! cvflag_grav
2373 fr(il, 1) = 0.1*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + &
2374 sigd*0.5*(evap(il,1)+evap(il,2))
2375 fr(il, 1) = fr(il, 1) + 0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
2376 fu(il, 1) = fu(il, 1) + 0.1*work(il)*(mp(il,2)*(up(il,2)-u(il, &
2377 1))+am(il)*(u(il,2)-u(il,1)))
2378 fv(il, 1) = fv(il, 1) + 0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
2379 1))+am(il)*(v(il,2)-v(il,1)))
2380 END IF ! cvflag_grav
2381
2382 END DO ! il
2383
2384 ! do j=1,ntra
2385 ! do il=1,ncum
2386 ! if (cvflag_grav) then
2387 ! ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
2388 ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2389 ! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2390 ! else
2391 ! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2392 ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2393 ! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2394 ! endif
2395 ! enddo
2396 ! enddo
2397
2398 DO j = 2, nl
2399 DO il = 1, ncum
2400 IF (j<=inb(il)) THEN
2401 IF (cvflag_grav) THEN
2402 fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il, &
2403 j,1)-rr(il,1))
2404 fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il, &
2405 j,1)-u(il,1))
2406 fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il, &
2407 j,1)-v(il,1))
2408 ELSE ! cvflag_grav
2409 fr(il, 1) = fr(il, 1) + 0.1*work(il)*ment(il, j, 1)*(qent(il,j,1)- &
2410 rr(il,1))
2411 fu(il, 1) = fu(il, 1) + 0.1*work(il)*ment(il, j, 1)*(uent(il,j,1)-u &
2412 (il,1))
2413 fv(il, 1) = fv(il, 1) + 0.1*work(il)*ment(il, j, 1)*(vent(il,j,1)-v &
2414 (il,1))
2415 END IF ! cvflag_grav
2416 END IF ! j
2417 END DO
2418 END DO
2419
2420 ! do k=1,ntra
2421 ! do j=2,nl
2422 ! do il=1,ncum
2423 ! if (j.le.inb(il)) then
2424
2425 ! if (cvflag_grav) then
2426 ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
2427 ! : *(traent(il,j,1,k)-tra(il,1,k))
2428 ! else
2429 ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
2430 ! : *(traent(il,j,1,k)-tra(il,1,k))
2431 ! endif
2432
2433 ! endif
2434 ! enddo
2435 ! enddo
2436 ! enddo
2437
2438
2439 ! *** calculate tendencies of potential temperature and mixing ratio ***
2440 ! *** at levels above the lowest level ***
2441
2442 ! *** first find the net saturated updraft and downdraft mass fluxes ***
2443 ! *** through each level ***
2444
2445
2446 DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
2447
2448 num1 = 0
2449 DO il = 1, ncum
2450 IF (i<=inb(il)) num1 = num1 + 1
2451 END DO
2452 IF (num1<=0) GO TO 500
2453
2454 CALL zilch(amp1, ncum)
2455 CALL zilch(ad, ncum)
2456
2457 DO k = i + 1, nl + 1
2458 DO il = 1, ncum
2459 IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN
2460 amp1(il) = amp1(il) + m(il, k)
2461 END IF
2462 END DO
2463 END DO
2464
2465 DO k = 1, i
2466 DO j = i + 1, nl + 1
2467 DO il = 1, ncum
2468 IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
2469 amp1(il) = amp1(il) + ment(il, k, j)
2470 END IF
2471 END DO
2472 END DO
2473 END DO
2474
2475 DO k = 1, i - 1
2476 DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
2477 DO il = 1, ncum
2478 IF (i<=inb(il) .AND. j<=inb(il)) THEN
2479 ad(il) = ad(il) + ment(il, j, k)
2480 END IF
2481 END DO
2482 END DO
2483 END DO
2484
2485 DO il = 1, ncum
2486 IF (i<=inb(il)) THEN
2487 dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2488 cpinv = 1.0/cpn(il, i)
2489
2490 ! convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
2491 IF (cvflag_grav) THEN
2492 IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2493 ELSE
2494 IF ((0.1*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2495 END IF
2496
2497 IF (cvflag_grav) THEN
2498 ft(il, i) = 0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
2499 i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
2500 i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( &
2501 il,i)+evap(il,i+1))
2502 rat = cpn(il, i-1)*cpinv
2503 ft(il, i) = ft(il, i) - 0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
2504 -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2505 ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h( &
2506 il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2507 ELSE ! cvflag_grav
2508 ft(il, i) = 0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
2509 i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
2510 i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( &
2511 il,i)+evap(il,i+1))
2512 rat = cpn(il, i-1)*cpinv
2513 ft(il, i) = ft(il, i) - 0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i)-mp(il &
2514 ,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2515 ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i)+ &
2516 t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2517 END IF ! cvflag_grav
2518
2519
2520 ft(il, i) = ft(il, i) + 0.01*sigd*wt(il, i)*(cl-cpd)*water(il, i+1)*( &
2521 t(il,i+1)-t(il,i))*dpinv*cpinv
2522
2523 IF (cvflag_grav) THEN
2524 fr(il, i) = 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
2525 i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2526 fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
2527 i))-ad(il)*(u(il,i)-u(il,i-1)))
2528 fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
2529 i))-ad(il)*(v(il,i)-v(il,i-1)))
2530 ELSE ! cvflag_grav
2531 fr(il, i) = 0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
2532 i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2533 fu(il, i) = fu(il, i) + 0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
2534 i))-ad(il)*(u(il,i)-u(il,i-1)))
2535 fv(il, i) = fv(il, i) + 0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
2536 i))-ad(il)*(v(il,i)-v(il,i-1)))
2537 END IF ! cvflag_grav
2538
2539 END IF ! i
2540 END DO
2541
2542 ! do k=1,ntra
2543 ! do il=1,ncum
2544 ! if (i.le.inb(il)) then
2545 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2546 ! cpinv=1.0/cpn(il,i)
2547 ! if (cvflag_grav) then
2548 ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
2549 ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2550 ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2551 ! else
2552 ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2553 ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2554 ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2555 ! endif
2556 ! endif
2557 ! enddo
2558 ! enddo
2559
2560 DO k = 1, i - 1
2561 DO il = 1, ncum
2562 IF (i<=inb(il)) THEN
2563 dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2564 cpinv = 1.0/cpn(il, i)
2565
2566 awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
2567 awat = amax1(awat, 0.0)
2568
2569 IF (cvflag_grav) THEN