source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phylmd/cv30_routines.F90 @ 5360

Last change on this file since 5360 was 2520, checked in by lguez, 9 years ago

Bug fix: icb1 was not made >= 2.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 98.3 KB
Line 
1
2! $Id: cv30_routines.F90 2520 2016-05-23 08:32:55Z abarral $
3
4
5
6SUBROUTINE 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
83END SUBROUTINE cv30_param
84
85SUBROUTINE 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
158END SUBROUTINE cv30_prelim
159
160SUBROUTINE 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
341END SUBROUTINE cv30_feed
342
343SUBROUTINE 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
623END SUBROUTINE cv30_undilute1
624
625SUBROUTINE 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
727END SUBROUTINE cv30_trigger
728
729SUBROUTINE 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
835END SUBROUTINE cv30_compress
836
837SUBROUTINE 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
1210END SUBROUTINE cv30_undilute2
1211
1212SUBROUTINE 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
1415END SUBROUTINE cv30_closure
1416
1417SUBROUTINE 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
1701175 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
1780789 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
1811END SUBROUTINE cv30_mixing
1812
1813
1814SUBROUTINE 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
2174400 END DO
2175
2176  RETURN
2177END SUBROUTINE cv30_unsat
2178
2179SUBROUTINE 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
2727500 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
3029END SUBROUTINE cv30_yield
3030
3031! !RomP >>>
3032SUBROUTINE 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
3053  REAL epm(nloc, na, na)
3054
3055  ! variables d'Emanuel : du second indice au troisieme
3056  ! --->    tab(i,k,j) -> de l origine k a l arrivee j
3057  ! ment, sij, elij
3058  ! variables personnelles : du troisieme au second indice
3059  ! --->    tab(i,j,k) -> de k a j
3060  ! phi, phi2
3061
3062  ! initialisations
3063  DO j = 1, na
3064    DO i = 1, ncum
3065      da(i, j) = 0.
3066      d1a(i, j) = 0.
3067      dam(i, j) = 0.
3068      eplamm(i, j) = 0.
3069    END DO
3070  END DO
3071  DO k = 1, na
3072    DO j = 1, na
3073      DO i = 1, ncum
3074        epm(i, j, k) = 0.
3075        epmlmmm(i, j, k) = 0.
3076        phi(i, j, k) = 0.
3077        phi2(i, j, k) = 0.
3078      END DO
3079    END DO
3080  END DO
3081
3082  ! fraction deau condensee dans les melanges convertie en precip : epm
3083  ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3084  DO j = 1, na
3085    DO k = 1, j - 1
3086      DO i = 1, ncum
3087        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3088          ! !jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3089          epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3090          ! !
3091          epm(i, j, k) = max(epm(i,j,k), 0.0)
3092        END IF
3093      END DO
3094    END DO
3095  END DO
3096
3097  DO j = 1, na
3098    DO k = 1, na
3099      DO i = 1, ncum
3100        IF (k>=icb(i) .AND. k<=inb(i)) THEN
3101          eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- &
3102            sij(i,j,k))
3103        END IF
3104      END DO
3105    END DO
3106  END DO
3107
3108  DO j = 1, na
3109    DO k = 1, j - 1
3110      DO i = 1, ncum
3111        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3112          epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
3113        END IF
3114      END DO
3115    END DO
3116  END DO
3117
3118  ! matrices pour calculer la tendance des concentrations dans cvltr.F90
3119  DO j = 1, na
3120    DO k = 1, na
3121      DO i = 1, ncum
3122        da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
3123        phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
3124        d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sij(i,k,j))
3125      END DO
3126    END DO
3127  END DO
3128
3129  DO j = 1, na
3130    DO k = 1, j - 1
3131      DO i = 1, ncum
3132        dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, j, k)*(1.-ep(i,k))*(1.- &
3133          sij(i,k,j))
3134        phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
3135      END DO
3136    END DO
3137  END DO
3138
3139  RETURN
3140END SUBROUTINE cv30_tracer
3141! RomP <<<
3142
3143SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
3144    vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
3145    dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
3146    epmlmmm, eplamm, wdtraina, wdtrainm,epmax_diag, iflag1, precip1, vprecip1, evap1, &
3147    ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
3148    dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
3149    elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1,epmax_diag1) ! epmax_cape
3150  IMPLICIT NONE
3151
3152  include "cv30param.h"
3153
3154  ! inputs:
3155  INTEGER len, ncum, nd, ntra, nloc
3156  INTEGER idcum(nloc)
3157  INTEGER iflag(nloc)
3158  INTEGER inb(nloc)
3159  REAL precip(nloc)
3160  REAL vprecip(nloc, nd+1), evap(nloc, nd)
3161  REAL ep(nloc, nd)
3162  REAL sig(nloc, nd), w0(nloc, nd)
3163  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3164  REAL ftra(nloc, nd, ntra)
3165  REAL ma(nloc, nd)
3166  REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3167  REAL qcondc(nloc, nd)
3168  REAL wd(nloc), cape(nloc)
3169  REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
3170  REAL epmax_diag(nloc) ! epmax_cape
3171  ! RomP >>>
3172  REAL phi2(nloc, nd, nd)
3173  REAL d1a(nloc, nd), dam(nloc, nd)
3174  REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
3175  REAL sij(nloc, nd, nd)
3176  REAL elij(nloc, nd, nd), clw(nloc, nd)
3177  REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
3178  ! RomP <<<
3179
3180  ! outputs:
3181  INTEGER iflag1(len)
3182  INTEGER inb1(len)
3183  REAL precip1(len)
3184  REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP
3185  REAL ep1(len, nd) !<<< RomP
3186  REAL sig1(len, nd), w01(len, nd)
3187  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
3188  REAL ftra1(len, nd, ntra)
3189  REAL ma1(len, nd)
3190  REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
3191  REAL qcondc1(nloc, nd)
3192  REAL wd1(nloc), cape1(nloc)
3193  REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
3194  REAL epmax_diag1(len) ! epmax_cape
3195  ! RomP >>>
3196  REAL phi21(len, nd, nd)
3197  REAL d1a1(len, nd), dam1(len, nd)
3198  REAL wdtraina1(len, nd), wdtrainm1(len, nd)
3199  REAL sij1(len, nd, nd)
3200  REAL elij1(len, nd, nd), clw1(len, nd)
3201  REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
3202  ! RomP <<<
3203
3204  ! local variables:
3205  INTEGER i, k, j
3206
3207  DO i = 1, ncum
3208    precip1(idcum(i)) = precip(i)
3209    iflag1(idcum(i)) = iflag(i)
3210    wd1(idcum(i)) = wd(i)
3211    inb1(idcum(i)) = inb(i)
3212    cape1(idcum(i)) = cape(i)
3213    epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
3214  END DO
3215
3216  DO k = 1, nl
3217    DO i = 1, ncum
3218      vprecip1(idcum(i), k) = vprecip(i, k)
3219      evap1(idcum(i), k) = evap(i, k) !<<< RomP
3220      sig1(idcum(i), k) = sig(i, k)
3221      w01(idcum(i), k) = w0(i, k)
3222      ft1(idcum(i), k) = ft(i, k)
3223      fq1(idcum(i), k) = fq(i, k)
3224      fu1(idcum(i), k) = fu(i, k)
3225      fv1(idcum(i), k) = fv(i, k)
3226      ma1(idcum(i), k) = ma(i, k)
3227      upwd1(idcum(i), k) = upwd(i, k)
3228      dnwd1(idcum(i), k) = dnwd(i, k)
3229      dnwd01(idcum(i), k) = dnwd0(i, k)
3230      qcondc1(idcum(i), k) = qcondc(i, k)
3231      da1(idcum(i), k) = da(i, k)
3232      mp1(idcum(i), k) = mp(i, k)
3233      ! RomP >>>
3234      ep1(idcum(i), k) = ep(i, k)
3235      d1a1(idcum(i), k) = d1a(i, k)
3236      dam1(idcum(i), k) = dam(i, k)
3237      clw1(idcum(i), k) = clw(i, k)
3238      eplamm1(idcum(i), k) = eplamm(i, k)
3239      wdtraina1(idcum(i), k) = wdtraina(i, k)
3240      wdtrainm1(idcum(i), k) = wdtrainm(i, k)
3241      ! RomP <<<
3242    END DO
3243  END DO
3244
3245  DO i = 1, ncum
3246    sig1(idcum(i), nd) = sig(i, nd)
3247  END DO
3248
3249
3250  ! do 2100 j=1,ntra
3251  ! do 2110 k=1,nd ! oct3
3252  ! do 2120 i=1,ncum
3253  ! ftra1(idcum(i),k,j)=ftra(i,k,j)
3254  ! 2120     continue
3255  ! 2110    continue
3256  ! 2100   continue
3257  DO j = 1, nd
3258    DO k = 1, nd
3259      DO i = 1, ncum
3260        sij1(idcum(i), k, j) = sij(i, k, j)
3261        phi1(idcum(i), k, j) = phi(i, k, j)
3262        phi21(idcum(i), k, j) = phi2(i, k, j)
3263        elij1(idcum(i), k, j) = elij(i, k, j)
3264        epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
3265      END DO
3266    END DO
3267  END DO
3268
3269  RETURN
3270END SUBROUTINE cv30_uncompress
3271
3272        subroutine cv30_epmax_fn_cape(nloc,ncum,nd &
3273                ,cape,ep,hp,icb,inb,clw,nk,t,h,lv &
3274                ,epmax_diag)
3275        implicit none
3276
3277        ! On fait varier epmax en fn de la cape
3278        ! Il faut donc recalculer ep, et hp qui a déjà été calculé et
3279        ! qui en dépend
3280        ! Toutes les autres variables fn de ep sont calculées plus bas.
3281
3282#include "cvthermo.h"
3283#include "cv30param.h"
3284#include "conema3.h"
3285
3286! inputs:
3287      integer ncum, nd, nloc
3288      integer icb(nloc), inb(nloc)
3289      real cape(nloc)
3290      real clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd)
3291      integer nk(nloc)
3292! inouts:
3293      real ep(nloc,nd)
3294      real hp(nloc,nd)
3295! outputs ou local
3296      real epmax_diag(nloc)
3297! locals
3298      integer i,k   
3299      real hp_bak(nloc,nd)
3300
3301        ! on recalcule ep et hp
3302       
3303        if (coef_epmax_cape.gt.1e-12) then
3304        do i=1,ncum
3305           epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
3306           do k=1,nl
3307                ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
3308                ep(i,k)=amax1(ep(i,k),0.0)
3309                ep(i,k)=amin1(ep(i,k),epmax_diag(i))
3310           enddo
3311        enddo
3312
3313! On recalcule hp:
3314      do k=1,nl
3315        do i=1,ncum
3316          hp_bak(i,k)=hp(i,k)
3317        enddo
3318      enddo
3319      do k=1,nlp
3320        do i=1,ncum
3321          hp(i,k)=h(i,k)
3322        enddo
3323      enddo
3324      do k=minorig+1,nl
3325       do i=1,ncum
3326        if((k.ge.icb(i)).and.(k.le.inb(i)))then
3327          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
3328        endif
3329       enddo
3330      enddo !do k=minorig+1,n
3331!     write(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
3332      do i=1,ncum 
3333       do k=1,nl
3334        if (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) then
3335           write(*,*) 'i,k=',i,k
3336           write(*,*) 'coef_epmax_cape=',coef_epmax_cape
3337           write(*,*) 'epmax_diag(i)=',epmax_diag(i)
3338           write(*,*) 'ep(i,k)=',ep(i,k)
3339           write(*,*) 'hp(i,k)=',hp(i,k)
3340           write(*,*) 'hp_bak(i,k)=',hp_bak(i,k)
3341           write(*,*) 'h(i,k)=',h(i,k)
3342           write(*,*) 'nk(i)=',nk(i)
3343           write(*,*) 'h(i,nk(i))=',h(i,nk(i))
3344           write(*,*) 'lv(i,k)=',lv(i,k)
3345           write(*,*) 't(i,k)=',t(i,k)
3346           write(*,*) 'clw(i,k)=',clw(i,k)
3347           write(*,*) 'cpd,cpv=',cpd,cpv
3348           stop
3349        endif
3350       enddo !do k=1,nl
3351      enddo !do i=1,ncum 
3352      endif !if (coef_epmax_cape.gt.1e-12) then
3353
3354      return
3355      end subroutine cv30_epmax_fn_cape
3356
3357
Note: See TracBrowser for help on using the repository browser.