GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/cv30_routines.F90 Lines: 0 1077 0.0 %
Date: 2023-06-30 12:56:34 Branches: 0 942 0.0 %

Line Branch Exec Source
1
2
! $Id: cv30_routines.F90 4593 2023-06-29 13:55:54Z ymeurdesoif $
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
        INCLUDE "cvthermo.h"
3284
        INCLUDE "cv30param.h"
3285
        INCLUDE "conema3.h"
3286
3287
! inputs:
3288
      integer ncum, nd, nloc
3289
      integer icb(nloc), inb(nloc)
3290
      real cape(nloc)
3291
      real clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd)
3292
      integer nk(nloc)
3293
! inouts:
3294
      real ep(nloc,nd)
3295
      real hp(nloc,nd)
3296
! outputs ou local
3297
      real epmax_diag(nloc)
3298
! locals
3299
      integer i,k
3300
      real hp_bak(nloc,nd)
3301
      CHARACTER (LEN=20) :: modname='cv30_epmax_fn_cape'
3302
      CHARACTER (LEN=80) :: abort_message
3303
3304
        ! on recalcule ep et hp
3305
3306
        if (coef_epmax_cape.gt.1e-12) then
3307
        do i=1,ncum
3308
           epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
3309
           do k=1,nl
3310
                ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
3311
                ep(i,k)=amax1(ep(i,k),0.0)
3312
                ep(i,k)=amin1(ep(i,k),epmax_diag(i))
3313
           enddo
3314
        enddo
3315
3316
! On recalcule hp:
3317
      do k=1,nl
3318
        do i=1,ncum
3319
	  hp_bak(i,k)=hp(i,k)
3320
	enddo
3321
      enddo
3322
      do k=1,nlp
3323
        do i=1,ncum
3324
	  hp(i,k)=h(i,k)
3325
	enddo
3326
      enddo
3327
      do k=minorig+1,nl
3328
       do i=1,ncum
3329
        if((k.ge.icb(i)).and.(k.le.inb(i)))then
3330
          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
3331
        endif
3332
       enddo
3333
      enddo !do k=minorig+1,n
3334
!     write(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
3335
      do i=1,ncum
3336
       do k=1,nl
3337
        if (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) then
3338
           write(*,*) 'i,k=',i,k
3339
           write(*,*) 'coef_epmax_cape=',coef_epmax_cape
3340
           write(*,*) 'epmax_diag(i)=',epmax_diag(i)
3341
           write(*,*) 'ep(i,k)=',ep(i,k)
3342
           write(*,*) 'hp(i,k)=',hp(i,k)
3343
           write(*,*) 'hp_bak(i,k)=',hp_bak(i,k)
3344
           write(*,*) 'h(i,k)=',h(i,k)
3345
           write(*,*) 'nk(i)=',nk(i)
3346
           write(*,*) 'h(i,nk(i))=',h(i,nk(i))
3347
           write(*,*) 'lv(i,k)=',lv(i,k)
3348
           write(*,*) 't(i,k)=',t(i,k)
3349
           write(*,*) 'clw(i,k)=',clw(i,k)
3350
           write(*,*) 'cpd,cpv=',cpd,cpv
3351
           CALL abort_physic(modname,abort_message,1)
3352
        endif
3353
       enddo !do k=1,nl
3354
      enddo !do i=1,ncum
3355
      endif !if (coef_epmax_cape.gt.1e-12) then
3356
3357
      return
3358
      end subroutine cv30_epmax_fn_cape
3359
3360