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
2570 fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
2571 ,i)-awat-rr(il,i))
2572 fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2573 ,i)-u(il,i))
2574 fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
2575 ,i)-v(il,i))
2576 ELSE ! cvflag_grav
2577 fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)- &
2578 awat-rr(il,i))
2579 fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2580 ,i)-u(il,i))
2581 fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
2582 il,i))
2583 END IF ! cvflag_grav
2584
2585 ! (saturated updrafts resulting from mixing) ! cld
2586 qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld
2587 nqcond(il, i) = nqcond(il, i) + 1. ! cld
2588 END IF ! i
2589 END DO
2590 END DO
2591
2592 ! do j=1,ntra
2593 ! do k=1,i-1
2594 ! do il=1,ncum
2595 ! if (i.le.inb(il)) then
2596 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2597 ! cpinv=1.0/cpn(il,i)
2598 ! if (cvflag_grav) then
2599 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2600 ! : *(traent(il,k,i,j)-tra(il,i,j))
2601 ! else
2602 ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2603 ! : *(traent(il,k,i,j)-tra(il,i,j))
2604 ! endif
2605 ! endif
2606 ! enddo
2607 ! enddo
2608 ! enddo
2609
2610 DO k = i, nl + 1
2611 DO il = 1, ncum
2612 IF (i<=inb(il) .AND. k<=inb(il)) THEN
2613 dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2614 cpinv = 1.0/cpn(il, i)
2615
2616 IF (cvflag_grav) THEN
2617 fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
2618 ,i)-rr(il,i))
2619 fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2620 ,i)-u(il,i))
2621 fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
2622 ,i)-v(il,i))
2623 ELSE ! cvflag_grav
2624 fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)-rr &
2625 (il,i))
2626 fu(il, i) = fu(il, i) + 0.1*dpinv*ment(il, k, i)*(uent(il,k,i)-u( &
2627 il,i))
2628 fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
2629 il,i))
2630 END IF ! cvflag_grav
2631 END IF ! i and k
2632 END DO
2633 END DO
2634
2635 ! do j=1,ntra
2636 ! do k=i,nl+1
2637 ! do il=1,ncum
2638 ! if (i.le.inb(il) .and. k.le.inb(il)) then
2639 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2640 ! cpinv=1.0/cpn(il,i)
2641 ! if (cvflag_grav) then
2642 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2643 ! : *(traent(il,k,i,j)-tra(il,i,j))
2644 ! else
2645 ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2646 ! : *(traent(il,k,i,j)-tra(il,i,j))
2647 ! endif
2648 ! endif ! i and k
2649 ! enddo
2650 ! enddo
2651 ! enddo
2652
2653 DO il = 1, ncum
2654 IF (i<=inb(il)) THEN
2655 dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2656 cpinv = 1.0/cpn(il, i)
2657
2658 IF (cvflag_grav) THEN
2659 ! sb: on ne fait pas encore la correction permettant de mieux
2660 ! conserver l'eau:
2661 fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + &
2662 0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il, &
2663 i)-rr(il,i-1)))*dpinv
2664
2665 fu(il, i) = fu(il, i) + 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il, &
2666 i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2667 fv(il, i) = fv(il, i) + 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il, &
2668 i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2669 ELSE ! cvflag_grav
2670 fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + &
2671 0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, &
2672 i-1)))*dpinv
2673 fu(il, i) = fu(il, i) + 0.1*(mp(il,i+1)*(up(il,i+1)-u(il, &
2674 i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2675 fv(il, i) = fv(il, i) + 0.1*(mp(il,i+1)*(vp(il,i+1)-v(il, &
2676 i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2677 END IF ! cvflag_grav
2678
2679 END IF ! i
2680 END DO
2681
2682 ! sb: interface with the cloud parameterization: ! cld
2683
2684 DO k = i + 1, nl
2685 DO il = 1, ncum
2686 IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
2687 ! (saturated downdrafts resulting from mixing) ! cld
2688 qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
2689 nqcond(il, i) = nqcond(il, i) + 1. ! cld
2690 END IF ! cld
2691 END DO ! cld
2692 END DO ! cld
2693
2694 ! (particular case: no detraining level is found) ! cld
2695 DO il = 1, ncum ! cld
2696 IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld
2697 qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
2698 nqcond(il, i) = nqcond(il, i) + 1. ! cld
2699 END IF ! cld
2700 END DO ! cld
2701
2702 DO il = 1, ncum ! cld
2703 IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld
2704 qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
2705 END IF ! cld
2706 END DO
2707
2708 ! do j=1,ntra
2709 ! do il=1,ncum
2710 ! if (i.le.inb(il)) then
2711 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2712 ! cpinv=1.0/cpn(il,i)
2713
2714 ! if (cvflag_grav) then
2715 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
2716 ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2717 ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2718 ! else
2719 ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2720 ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2721 ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2722 ! endif
2723 ! endif ! i
2724 ! enddo
2725 ! enddo
2726
2727 500 END DO
2728
2729
2730 ! *** move the detrainment at level inb down to level inb-1 ***
2731 ! *** in such a way as to preserve the vertically ***
2732 ! *** integrated enthalpy and water tendencies ***
2733
2734 DO il = 1, ncum
2735
2736 ax = 0.1*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t(il, &
2737 inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
2738 inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2739 ft(il, inb(il)) = ft(il, inb(il)) - ax
2740 ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il &
2741 ))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, &
2742 inb(il))))
2743
2744 bx = 0.1*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb( &
2745 il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2746 fr(il, inb(il)) = fr(il, inb(il)) - bx
2747 fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ &
2748 1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2749
2750 cx = 0.1*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il &
2751 )))/(ph(il,inb(il))-ph(il,inb(il)+1))
2752 fu(il, inb(il)) = fu(il, inb(il)) - cx
2753 fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ &
2754 1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2755
2756 dx = 0.1*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il &
2757 )))/(ph(il,inb(il))-ph(il,inb(il)+1))
2758 fv(il, inb(il)) = fv(il, inb(il)) - dx
2759 fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ &
2760 1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2761
2762 END DO
2763
2764 ! do j=1,ntra
2765 ! do il=1,ncum
2766 ! ex=0.1*ment(il,inb(il),inb(il))
2767 ! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2768 ! : /(ph(il,inb(il))-ph(il,inb(il)+1))
2769 ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2770 ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2771 ! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2772 ! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2773 ! enddo
2774 ! enddo
2775
2776
2777 ! *** homoginize tendencies below cloud base ***
2778
2779
2780 DO il = 1, ncum
2781 asum(il) = 0.0
2782 bsum(il) = 0.0
2783 csum(il) = 0.0
2784 dsum(il) = 0.0
2785 END DO
2786
2787 DO i = 1, nl
2788 DO il = 1, ncum
2789 IF (i<=(icb(il)-1)) THEN
2790 asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1))
2791 bsum(il) = bsum(il) + fr(il, i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2792 1)))*(ph(il,i)-ph(il,i+1))
2793 csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2794 1)))*(ph(il,i)-ph(il,i+1))
2795 dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
2796 END IF
2797 END DO
2798 END DO
2799
2800 ! !!! do 700 i=1,icb(il)-1
2801 DO i = 1, nl
2802 DO il = 1, ncum
2803 IF (i<=(icb(il)-1)) THEN
2804 ft(il, i) = asum(il)*t(il, i)/(th(il,i)*dsum(il))
2805 fr(il, i) = bsum(il)/csum(il)
2806 END IF
2807 END DO
2808 END DO
2809
2810
2811 ! *** reset counter and return ***
2812
2813 DO il = 1, ncum
2814 sig(il, nd) = 2.0
2815 END DO
2816
2817
2818 DO i = 1, nd
2819 DO il = 1, ncum
2820 upwd(il, i) = 0.0
2821 dnwd(il, i) = 0.0
2822 END DO
2823 END DO
2824
2825 DO i = 1, nl
2826 DO il = 1, ncum
2827 dnwd0(il, i) = -mp(il, i)
2828 END DO
2829 END DO
2830 DO i = nl + 1, nd
2831 DO il = 1, ncum
2832 dnwd0(il, i) = 0.
2833 END DO
2834 END DO
2835
2836
2837 DO i = 1, nl
2838 DO il = 1, ncum
2839 IF (i>=icb(il) .AND. i<=inb(il)) THEN
2840 upwd(il, i) = 0.0
2841 dnwd(il, i) = 0.0
2842 END IF
2843 END DO
2844 END DO
2845
2846 DO i = 1, nl
2847 DO k = 1, nl
2848 DO il = 1, ncum
2849 up1(il, k, i) = 0.0
2850 dn1(il, k, i) = 0.0
2851 END DO
2852 END DO
2853 END DO
2854
2855 DO i = 1, nl
2856 DO k = i, nl
2857 DO n = 1, i - 1
2858 DO il = 1, ncum
2859 IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
2860 up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
2861 dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
2862 END IF
2863 END DO
2864 END DO
2865 END DO
2866 END DO
2867
2868 DO i = 2, nl
2869 DO k = i, nl
2870 DO il = 1, ncum
2871 ! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
2872 ! then
2873 IF (i<=inb(il) .AND. k<=inb(il)) THEN
2874 upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
2875 dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
2876 END IF
2877 END DO
2878 END DO
2879 END DO
2880
2881
2882 ! !!! DO il=1,ncum
2883 ! !!! do i=icb(il),inb(il)
2884 ! !!!
2885 ! !!! upwd(il,i)=0.0
2886 ! !!! dnwd(il,i)=0.0
2887 ! !!! do k=i,inb(il)
2888 ! !!! up1=0.0
2889 ! !!! dn1=0.0
2890 ! !!! do n=1,i-1
2891 ! !!! up1=up1+ment(il,n,k)
2892 ! !!! dn1=dn1-ment(il,k,n)
2893 ! !!! enddo
2894 ! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2895 ! !!! dnwd(il,i)=dnwd(il,i)+dn1
2896 ! !!! enddo
2897 ! !!! enddo
2898 ! !!!
2899 ! !!! ENDDO
2900
2901 ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2902 ! determination de la variation de flux ascendant entre
2903 ! deux niveau non dilue mike
2904 ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2905
2906 DO i = 1, nl
2907 DO il = 1, ncum
2908 mike(il, i) = m(il, i)
2909 END DO
2910 END DO
2911
2912 DO i = nl + 1, nd
2913 DO il = 1, ncum
2914 mike(il, i) = 0.
2915 END DO
2916 END DO
2917
2918 DO i = 1, nd
2919 DO il = 1, ncum
2920 ma(il, i) = 0
2921 END DO
2922 END DO
2923
2924 DO i = 1, nl
2925 DO j = i, nl
2926 DO il = 1, ncum
2927 ma(il, i) = ma(il, i) + m(il, j)
2928 END DO
2929 END DO
2930 END DO
2931
2932 DO i = nl + 1, nd
2933 DO il = 1, ncum
2934 ma(il, i) = 0.
2935 END DO
2936 END DO
2937
2938 DO i = 1, nl
2939 DO il = 1, ncum
2940 IF (i<=(icb(il)-1)) THEN
2941 ma(il, i) = 0
2942 END IF
2943 END DO
2944 END DO
2945
2946 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2947 ! icb represente de niveau ou se trouve la
2948 ! base du nuage , et inb le top du nuage
2949 ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2950
2951 DO i = 1, nd
2952 DO il = 1, ncum
2953 mke(il, i) = upwd(il, i) + dnwd(il, i)
2954 END DO
2955 END DO
2956
2957 DO i = 1, nd
2958 DO il = 1, ncum
2959 rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il, &
2960 i))+rr(il,i)*cpv)
2961 tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp
2962 tps(il, i) = tp(il, i)
2963 END DO
2964 END DO
2965
2966
2967 ! *** diagnose the in-cloud mixing ratio *** ! cld
2968 ! *** of condensed water *** ! cld
2969 ! ! cld
2970
2971 DO i = 1, nd ! cld
2972 DO il = 1, ncum ! cld
2973 mac(il, i) = 0.0 ! cld
2974 wa(il, i) = 0.0 ! cld
2975 siga(il, i) = 0.0 ! cld
2976 sax(il, i) = 0.0 ! cld
2977 END DO ! cld
2978 END DO ! cld
2979
2980 DO i = minorig, nl ! cld
2981 DO k = i + 1, nl + 1 ! cld
2982 DO il = 1, ncum ! cld
2983 IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld
2984 mac(il, i) = mac(il, i) + m(il, k) ! cld
2985 END IF ! cld
2986 END DO ! cld
2987 END DO ! cld
2988 END DO ! cld
2989
2990 DO i = 1, nl ! cld
2991 DO j = 1, i ! cld
2992 DO il = 1, ncum ! cld
2993 IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
2994 .AND. j>=icb(il)) THEN ! cld
2995 sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
2996 *(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
2997 END IF ! cld
2998 END DO ! cld
2999 END DO ! cld
3000 END DO ! cld
3001
3002 DO i = 1, nl ! cld
3003 DO il = 1, ncum ! cld
3004 IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3005 .AND. sax(il,i)>0.0) THEN ! cld
3006 wa(il, i) = sqrt(2.*sax(il,i)) ! cld
3007 END IF ! cld
3008 END DO ! cld
3009 END DO ! cld
3010
3011 DO i = 1, nl ! cld
3012 DO il = 1, ncum ! cld
3013 IF (wa(il,i)>0.0) & ! cld
3014 siga(il, i) = mac(il, i)/wa(il, i) & ! cld
3015 *rrd*tvp(il, i)/p(il, i)/100./delta ! cld
3016 siga(il, i) = min(siga(il,i), 1.0) ! cld
3017 ! IM cf. FH
3018 IF (iflag_clw==0) THEN
3019 qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) & ! cld
3020 +(1.-siga(il,i))*qcond(il, i) ! cld
3021 ELSE IF (iflag_clw==1) THEN
3022 qcondc(il, i) = qcond(il, i) ! cld
3023 END IF
3024
3025 END DO ! cld
3026 END DO ! cld
3027
3028 RETURN
3029 END SUBROUTINE cv30_yield
3030
3031 ! !RomP >>>
3032 SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
3033 d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
3034 IMPLICIT NONE
3035
3036 include "cv30param.h"
3037
3038 ! inputs:
3039 INTEGER ncum, nd, na, nloc, len
3040 REAL ment(nloc, na, na), sij(nloc, na, na)
3041 REAL clw(nloc, nd), elij(nloc, na, na)
3042 REAL ep(nloc, na)
3043 INTEGER icb(nloc), inb(nloc)
3044 REAL vprecip(nloc, nd+1)
3045 ! ouputs:
3046 REAL da(nloc, na), phi(nloc, na, na)
3047 REAL phi2(nloc, na, na)
3048 REAL d1a(nloc, na), dam(nloc, na)
3049 REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
3050 ! variables pour tracer dans precip de l'AA et des mel
3051 ! local variables:
3052 INTEGER i, j, k, nam1
3053 REAL epm(nloc, na, na)
3054
3055 nam1=na-1 ! Introduced because ep is not defined for j=na
3056 ! variables d'Emanuel : du second indice au troisieme
3057 ! ---> tab(i,k,j) -> de l origine k a l arrivee j
3058 ! ment, sij, elij
3059 ! variables personnelles : du troisieme au second indice
3060 ! ---> tab(i,j,k) -> de k a j
3061 ! phi, phi2
3062
3063 ! initialisations
3064 DO j = 1, na
3065 DO i = 1, ncum
3066 da(i, j) = 0.
3067 d1a(i, j) = 0.
3068 dam(i, j) = 0.
3069 eplamm(i, j) = 0.
3070 END DO
3071 END DO
3072 DO k = 1, na
3073 DO j = 1, na
3074 DO i = 1, ncum
3075 epm(i, j, k) = 0.
3076 epmlmmm(i, j, k) = 0.
3077 phi(i, j, k) = 0.
3078 phi2(i, j, k) = 0.
3079 END DO
3080 END DO
3081 END DO
3082
3083 ! fraction deau condensee dans les melanges convertie en precip : epm
3084 ! et eau condens�e pr�cipit�e dans masse d'air satur� : l_m*dM_m/dzdz.dzdz
3085 DO j = 1, nam1
3086 DO k = 1, j - 1
3087 DO i = 1, ncum
3088 IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3089 ! !jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3090 epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3091 ! !
3092 epm(i, j, k) = max(epm(i,j,k), 0.0)
3093 END IF
3094 END DO
3095 END DO
3096 END DO
3097
3098 DO j = 1, nam1
3099 DO k = 1, nam1
3100 DO i = 1, ncum
3101 IF (k>=icb(i) .AND. k<=inb(i)) THEN
3102 eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- &
3103 sij(i,j,k))
3104 END IF
3105 END DO
3106 END DO
3107 END DO
3108
3109 DO j = 1, nam1
3110 DO k = 1, j - 1
3111 DO i = 1, ncum
3112 IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3113 epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
3114 END IF
3115 END DO
3116 END DO
3117 END DO
3118
3119 ! matrices pour calculer la tendance des concentrations dans cvltr.F90
3120 DO j = 1, nam1
3121 DO k = 1, nam1
3122 DO i = 1, ncum
3123 da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
3124 phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
3125 d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sij(i,k,j))
3126 END DO
3127 END DO
3128 END DO
3129
3130 DO j = 1, nam1
3131 DO k = 1, j - 1
3132 DO i = 1, ncum
3133 dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, j, k)*(1.-ep(i,k))*(1.- &
3134 sij(i,k,j))
3135 phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
3136 END DO
3137 END DO
3138 END DO
3139
3140 RETURN
3141 END SUBROUTINE cv30_tracer
3142 ! RomP <<<
3143
3144 SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
3145 vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
3146 dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
3147 epmlmmm, eplamm, wdtraina, wdtrainm,epmax_diag, iflag1, precip1, vprecip1, evap1, &
3148 ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
3149 dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
3150 elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1,epmax_diag1) ! epmax_cape
3151 IMPLICIT NONE
3152
3153 include "cv30param.h"
3154
3155 ! inputs:
3156 INTEGER len, ncum, nd, ntra, nloc
3157 INTEGER idcum(nloc)
3158 INTEGER iflag(nloc)
3159 INTEGER inb(nloc)
3160 REAL precip(nloc)
3161 REAL vprecip(nloc, nd+1), evap(nloc, nd)
3162 REAL ep(nloc, nd)
3163 REAL sig(nloc, nd), w0(nloc, nd)
3164 REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3165 REAL ftra(nloc, nd, ntra)
3166 REAL ma(nloc, nd)
3167 REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3168 REAL qcondc(nloc, nd)
3169 REAL wd(nloc), cape(nloc)
3170 REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
3171 REAL epmax_diag(nloc) ! epmax_cape
3172 ! RomP >>>
3173 REAL phi2(nloc, nd, nd)
3174 REAL d1a(nloc, nd), dam(nloc, nd)
3175 REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
3176 REAL sij(nloc, nd, nd)
3177 REAL elij(nloc, nd, nd), clw(nloc, nd)
3178 REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
3179 ! RomP <<<
3180
3181 ! outputs:
3182 INTEGER iflag1(len)
3183 INTEGER inb1(len)
3184 REAL precip1(len)
3185 REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP
3186 REAL ep1(len, nd) !<<< RomP
3187 REAL sig1(len, nd), w01(len, nd)
3188 REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
3189 REAL ftra1(len, nd, ntra)
3190 REAL ma1(len, nd)
3191 REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
3192 REAL qcondc1(nloc, nd)
3193 REAL wd1(nloc), cape1(nloc)
3194 REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
3195 REAL epmax_diag1(len) ! epmax_cape
3196 ! RomP >>>
3197 REAL phi21(len, nd, nd)
3198 REAL d1a1(len, nd), dam1(len, nd)
3199 REAL wdtraina1(len, nd), wdtrainm1(len, nd)
3200 REAL sij1(len, nd, nd)
3201 REAL elij1(len, nd, nd), clw1(len, nd)
3202 REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
3203 ! RomP <<<
3204
3205 ! local variables:
3206 INTEGER i, k, j
3207
3208 DO i = 1, ncum
3209 precip1(idcum(i)) = precip(i)
3210 iflag1(idcum(i)) = iflag(i)
3211 wd1(idcum(i)) = wd(i)
3212 inb1(idcum(i)) = inb(i)
3213 cape1(idcum(i)) = cape(i)
3214 epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
3215 END DO
3216
3217 DO k = 1, nl
3218 DO i = 1, ncum
3219 vprecip1(idcum(i), k) = vprecip(i, k)
3220 evap1(idcum(i), k) = evap(i, k) !<<< RomP
3221 sig1(idcum(i), k) = sig(i, k)
3222 w01(idcum(i), k) = w0(i, k)
3223 ft1(idcum(i), k) = ft(i, k)
3224 fq1(idcum(i), k) = fq(i, k)
3225 fu1(idcum(i), k) = fu(i, k)
3226 fv1(idcum(i), k) = fv(i, k)
3227 ma1(idcum(i), k) = ma(i, k)
3228 upwd1(idcum(i), k) = upwd(i, k)
3229 dnwd1(idcum(i), k) = dnwd(i, k)
3230 dnwd01(idcum(i), k) = dnwd0(i, k)
3231 qcondc1(idcum(i), k) = qcondc(i, k)
3232 da1(idcum(i), k) = da(i, k)
3233 mp1(idcum(i), k) = mp(i, k)
3234 ! RomP >>>
3235 ep1(idcum(i), k) = ep(i, k)
3236 d1a1(idcum(i), k) = d1a(i, k)
3237 dam1(idcum(i), k) = dam(i, k)
3238 clw1(idcum(i), k) = clw(i, k)
3239 eplamm1(idcum(i), k) = eplamm(i, k)
3240 wdtraina1(idcum(i), k) = wdtraina(i, k)
3241 wdtrainm1(idcum(i), k) = wdtrainm(i, k)
3242 ! RomP <<<
3243 END DO
3244 END DO
3245
3246 DO i = 1, ncum
3247 sig1(idcum(i), nd) = sig(i, nd)
3248 END DO
3249
3250
3251 ! do 2100 j=1,ntra
3252 ! do 2110 k=1,nd ! oct3
3253 ! do 2120 i=1,ncum
3254 ! ftra1(idcum(i),k,j)=ftra(i,k,j)
3255 ! 2120 continue
3256 ! 2110 continue
3257 ! 2100 continue
3258 DO j = 1, nd
3259 DO k = 1, nd
3260 DO i = 1, ncum
3261 sij1(idcum(i), k, j) = sij(i, k, j)
3262 phi1(idcum(i), k, j) = phi(i, k, j)
3263 phi21(idcum(i), k, j) = phi2(i, k, j)
3264 elij1(idcum(i), k, j) = elij(i, k, j)
3265 epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
3266 END DO
3267 END DO
3268 END DO
3269
3270 RETURN
3271 END SUBROUTINE cv30_uncompress
3272
3273 subroutine cv30_epmax_fn_cape(nloc,ncum,nd &
3274 ,cape,ep,hp,icb,inb,clw,nk,t,h,lv &
3275 ,epmax_diag)
3276 implicit none
3277
3278 ! On fait varier epmax en fn de la cape
3279 ! Il faut donc recalculer ep, et hp qui a d�j� �t� calcul� et
3280 ! qui en d�pend
3281 ! Toutes les autres variables fn de ep sont calcul�es plus bas.
3282
3283 !
3284 ! $Header$
3285 !
3286 ! Thermodynamical constants for convectL:
3287
3288 real cpd, cpv, cl, ci, rrv, rrd, lv0, lf0, g, rowl, t0
3289 real clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl, clmci
3290 real eps, epsi, epsim1
3291 real ginv, hrd
3292 real grav
3293
3294 COMMON /cvthermo/ cpd, cpv, cl, ci, rrv, rrd, lv0, lf0, g, rowl &
3295 ,t0, clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl &
3296 ,clmci, eps, epsi, epsim1, ginv, hrd, grav
3297
3298 !$OMP THREADPRIVATE(/cvthermo/)
3299 !
3300 ! $Header$
3301 !
3302 !------------------------------------------------------------
3303 ! Parameters for convectL, iflag_con=30:
3304 ! (includes - microphysical parameters,
3305 ! - parameters that control the rate of approach
3306 ! to quasi-equilibrium)
3307 ! - noff & minorig (previously in input of convect1)
3308 !------------------------------------------------------------
3309
3310 integer noff, minorig, nl, nlp, nlm
3311 real sigd, spfac
3312 !IM cf. FH : pour compatibilite avec conema3 TEMPORAIRE real pbcrit, ptcrit, epmax
3313 real pbcrit, ptcrit
3314 real omtrain
3315 real dtovsh, dpbase, dttrig
3316 real dtcrit, tau, beta, alpha
3317 real delta
3318 real betad
3319
3320 COMMON /cv30param/ noff, minorig, nl, nlp, nlm &
3321 , sigd, spfac &
3322 !IM cf. FH : pour compatibilite avec conema3 TEMPORAIRE : ,pbcrit, ptcrit, epmax
3323 ,pbcrit, ptcrit &
3324 ,omtrain &
3325 ,dtovsh, dpbase, dttrig &
3326 ,dtcrit, tau, beta, alpha, delta, betad
3327
3328 !$OMP THREADPRIVATE(/cv30param/)
3329 !
3330 ! $Header$
3331 !-- Modified by : Filiberti M-A 06/2005
3332 !
3333 real epmax ! 0.993
3334 real coef_epmax_cape ! 0.993
3335 !jyg<
3336 REAL cvl_comp_threshold ! 0.
3337 !>jyg
3338 logical ok_adj_ema ! F
3339 integer iflag_clw ! 0
3340 integer iflag_cvl_sigd
3341 real cvl_sig2feed ! 0.97
3342
3343 !jyg<
3344 !! common/comconema1/epmax,coef_epmax_cape,ok_adj_ema,iflag_clw,sig1feed,sig2feed
3345 !! common/comconema2/iflag_cvl_sigd
3346 common/comconema1/epmax,coef_epmax_cape, cvl_comp_threshold, cvl_sig2feed
3347 common/comconema2/iflag_cvl_sigd, iflag_clw, ok_adj_ema
3348 !>jyg
3349
3350 ! common/comconema/epmax,coef_epmax_cape,ok_adj_ema,iflag_clw
3351 !$OMP THREADPRIVATE(/comconema1/)
3352 !$OMP THREADPRIVATE(/comconema2/)
3353
3354
3355 ! inputs:
3356 integer ncum, nd, nloc
3357 integer icb(nloc), inb(nloc)
3358 real cape(nloc)
3359 real clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd)
3360 integer nk(nloc)
3361 ! inouts:
3362 real ep(nloc,nd)
3363 real hp(nloc,nd)
3364 ! outputs ou local
3365 real epmax_diag(nloc)
3366 ! locals
3367 integer i,k
3368 real hp_bak(nloc,nd)
3369 CHARACTER (LEN=20) :: modname='cv30_epmax_fn_cape'
3370 CHARACTER (LEN=80) :: abort_message
3371
3372 ! on recalcule ep et hp
3373
3374 if (coef_epmax_cape.gt.1e-12) then
3375 do i=1,ncum
3376 epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
3377 do k=1,nl
3378 ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
3379 ep(i,k)=amax1(ep(i,k),0.0)
3380 ep(i,k)=amin1(ep(i,k),epmax_diag(i))
3381 enddo
3382 enddo
3383
3384 ! On recalcule hp:
3385 do k=1,nl
3386 do i=1,ncum
3387 hp_bak(i,k)=hp(i,k)
3388 enddo
3389 enddo
3390 do k=1,nlp
3391 do i=1,ncum
3392 hp(i,k)=h(i,k)
3393 enddo
3394 enddo
3395 do k=minorig+1,nl
3396 do i=1,ncum
3397 if((k.ge.icb(i)).and.(k.le.inb(i)))then
3398 hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
3399 endif
3400 enddo
3401 enddo !do k=minorig+1,n
3402 ! write(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
3403 do i=1,ncum
3404 do k=1,nl
3405 if (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) then
3406 write(*,*) 'i,k=',i,k
3407 write(*,*) 'coef_epmax_cape=',coef_epmax_cape
3408 write(*,*) 'epmax_diag(i)=',epmax_diag(i)
3409 write(*,*) 'ep(i,k)=',ep(i,k)
3410 write(*,*) 'hp(i,k)=',hp(i,k)
3411 write(*,*) 'hp_bak(i,k)=',hp_bak(i,k)
3412 write(*,*) 'h(i,k)=',h(i,k)
3413 write(*,*) 'nk(i)=',nk(i)
3414 write(*,*) 'h(i,nk(i))=',h(i,nk(i))
3415 write(*,*) 'lv(i,k)=',lv(i,k)
3416 write(*,*) 't(i,k)=',t(i,k)
3417 write(*,*) 'clw(i,k)=',clw(i,k)
3418 write(*,*) 'cpd,cpv=',cpd,cpv
3419 CALL abort_physic(modname,abort_message,0)
3420 endif
3421 enddo !do k=1,nl
3422 enddo !do i=1,ncum
3423 endif !if (coef_epmax_cape.gt.1e-12) then
3424
3425 return
3426 end subroutine cv30_epmax_fn_cape
3427
3428
3429