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