source: LMDZ5/trunk/libf/phylmd/cv30_routines.F90 @ 2481

Last change on this file since 2481 was 2481, checked in by fhourdin, 9 years ago

Introduction d'une dependance epmax=f(Cape) sur proposition de Camille Risi
  • 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 2481 2016-03-29 13:39:55Z fhourdin $
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) = max(icb(i), 2) !convect3
410    icb1(i) = min(icb(i), nl) !convect3
411    ! if icb is below LCL, start loop at ICB+1:
412    ! (icbs est le premier niveau au-dessus du LCL)
413    icbs(i) = icb1(i) !convect3
414    IF (plcl(i)<p(i,icb1(i))) THEN
415      icbs(i) = min(icbs(i)+1, nl) !convect3
416    END IF
417  END DO !convect3
418
419  DO i = 1, len !convect3
420    ticb(i) = t(i, icbs(i)) !convect3
421    gzicb(i) = gz(i, icbs(i)) !convect3
422    qsicb(i) = qs(i, icbs(i)) !convect3
423  END DO !convect3
424
425
426  ! Re-compute icbsmax (icbsmax2):        !convect3
427  ! !convect3
428  icbsmax2 = 2 !convect3
429  DO i = 1, len !convect3
430    icbsmax2 = max(icbsmax2, icbs(i)) !convect3
431  END DO !convect3
432
433  ! initialization outputs:
434
435  DO k = 1, icbsmax2 ! convect3
436    DO i = 1, len ! convect3
437      tp(i, k) = 0.0 ! convect3
438      tvp(i, k) = 0.0 ! convect3
439      clw(i, k) = 0.0 ! convect3
440    END DO ! convect3
441  END DO ! convect3
442
443  ! tp and tvp below cloud base:
444
445  DO k = minorig, icbsmax2 - 1
446    DO i = 1, len
447      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
448      tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
449    END DO
450  END DO
451
452  ! ***  Find lifted parcel quantities above cloud base    ***
453
454  DO i = 1, len
455    tg = ticb(i)
456    ! ori         qg=qs(i,icb(i))
457    qg = qsicb(i) ! convect3
458    ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
459    alv = lv0 - clmcpv*(ticb(i)-273.15)
460
461    ! First iteration.
462
463    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
464    s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
465      +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
466    s = 1./s
467    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
468    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
469    tg = tg + s*(ah0(i)-ahg)
470    ! ori          tg=max(tg,35.0)
471    ! debug          tc=tg-t0
472    tc = tg - 273.15
473    denom = 243.5 + tc
474    denom = max(denom, 1.0) ! convect3
475    ! ori          if(tc.ge.0.0)then
476    es = 6.112*exp(17.67*tc/denom)
477    ! ori          else
478    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
479    ! ori          endif
480    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
481    qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
482
483    ! Second iteration.
484
485
486    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
487    ! ori          s=1./s
488    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
489    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
490    tg = tg + s*(ah0(i)-ahg)
491    ! ori          tg=max(tg,35.0)
492    ! debug          tc=tg-t0
493    tc = tg - 273.15
494    denom = 243.5 + tc
495    denom = max(denom, 1.0) ! convect3
496    ! ori          if(tc.ge.0.0)then
497    es = 6.112*exp(17.67*tc/denom)
498    ! ori          else
499    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
500    ! ori          end if
501    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
502    qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
503
504    alv = lv0 - clmcpv*(ticb(i)-273.15)
505
506    ! ori c approximation here:
507    ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
508    ! ori     &   -gz(i,icb(i))-alv*qg)/cpd
509
510    ! convect3: no approximation:
511    tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i))
512
513    ! ori         clw(i,icb(i))=qnk(i)-qg
514    ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
515    clw(i, icbs(i)) = qnk(i) - qg
516    clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
517
518    rg = qg/(1.-qnk(i))
519    ! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
520    ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
521    tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
522
523  END DO
524
525  ! ori      do 380 k=minorig,icbsmax2
526  ! ori       do 370 i=1,len
527  ! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
528  ! ori 370   continue
529  ! ori 380  continue
530
531
532  ! -- The following is only for convect3:
533
534  ! * icbs is the first level above the LCL:
535  ! if plcl<p(icb), then icbs=icb+1
536  ! if plcl>p(icb), then icbs=icb
537
538  ! * the routine above computes tvp from minorig to icbs (included).
539
540  ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
541  ! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
542
543  ! * therefore, in the case icbs=icb, we compute tvp at level icb+1
544  ! (tvp at other levels will be computed in cv3_undilute2.F)
545
546
547  DO i = 1, len
548    ticb(i) = t(i, icb(i)+1)
549    gzicb(i) = gz(i, icb(i)+1)
550    qsicb(i) = qs(i, icb(i)+1)
551  END DO
552
553  DO i = 1, len
554    tg = ticb(i)
555    qg = qsicb(i) ! convect3
556    ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
557    alv = lv0 - clmcpv*(ticb(i)-273.15)
558
559    ! First iteration.
560
561    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
562    s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
563      +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
564    s = 1./s
565    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
566    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
567    tg = tg + s*(ah0(i)-ahg)
568    ! ori          tg=max(tg,35.0)
569    ! debug          tc=tg-t0
570    tc = tg - 273.15
571    denom = 243.5 + tc
572    denom = max(denom, 1.0) ! convect3
573    ! ori          if(tc.ge.0.0)then
574    es = 6.112*exp(17.67*tc/denom)
575    ! ori          else
576    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
577    ! ori          endif
578    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
579    qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
580
581    ! Second iteration.
582
583
584    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
585    ! ori          s=1./s
586    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
587    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
588    tg = tg + s*(ah0(i)-ahg)
589    ! ori          tg=max(tg,35.0)
590    ! debug          tc=tg-t0
591    tc = tg - 273.15
592    denom = 243.5 + tc
593    denom = max(denom, 1.0) ! convect3
594    ! ori          if(tc.ge.0.0)then
595    es = 6.112*exp(17.67*tc/denom)
596    ! ori          else
597    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
598    ! ori          end if
599    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
600    qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
601
602    alv = lv0 - clmcpv*(ticb(i)-273.15)
603
604    ! ori c approximation here:
605    ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
606    ! ori     &   -gz(i,icb(i))-alv*qg)/cpd
607
608    ! convect3: no approximation:
609    tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
610
611    ! ori         clw(i,icb(i))=qnk(i)-qg
612    ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
613    clw(i, icb(i)+1) = qnk(i) - qg
614    clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
615
616    rg = qg/(1.-qnk(i))
617    ! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
618    ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
619    tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
620
621  END DO
622
623  RETURN
624END SUBROUTINE cv30_undilute1
625
626SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, &
627    iflag, sig, w0)
628  IMPLICIT NONE
629
630  ! -------------------------------------------------------------------
631  ! --- TRIGGERING
632
633  ! - computes the cloud base
634  ! - triggering (crude in this version)
635  ! - relaxation of sig and w0 when no convection
636
637  ! Caution1: if no convection, we set iflag=4
638  ! (it used to be 0 in convect3)
639
640  ! Caution2: at this stage, tvp (and thus buoy) are know up
641  ! through icb only!
642  ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
643  ! -------------------------------------------------------------------
644
645  include "cv30param.h"
646
647  ! input:
648  INTEGER len, nd
649  INTEGER icb(len)
650  REAL plcl(len), p(len, nd)
651  REAL th(len, nd), tv(len, nd), tvp(len, nd)
652
653  ! output:
654  REAL pbase(len), buoybase(len)
655
656  ! input AND output:
657  INTEGER iflag(len)
658  REAL sig(len, nd), w0(len, nd)
659
660  ! local variables:
661  INTEGER i, k
662  REAL tvpbase, tvbase, tdif, ath, ath1
663
664
665  ! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
666
667  DO i = 1, len
668    pbase(i) = plcl(i) + dpbase
669    tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
670      (p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
671      p(i,icb(i))-p(i,icb(i)+1))
672    tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
673      (p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
674      (i,icb(i))-p(i,icb(i)+1))
675    buoybase(i) = tvpbase - tvbase
676  END DO
677
678
679  ! ***   make sure that column is dry adiabatic between the surface  ***
680  ! ***    and cloud base, and that lifted air is positively buoyant  ***
681  ! ***                         at cloud base                         ***
682  ! ***       if not, return to calling program after resetting       ***
683  ! ***                        sig(i) and w0(i)                       ***
684
685
686  ! oct3      do 200 i=1,len
687  ! oct3
688  ! oct3       tdif = buoybase(i)
689  ! oct3       ath1 = th(i,1)
690  ! oct3       ath  = th(i,icb(i)-1) - dttrig
691  ! oct3
692  ! oct3       if (tdif.lt.dtcrit .or. ath.gt.ath1) then
693  ! oct3         do 60 k=1,nl
694  ! oct3            sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
695  ! oct3            sig(i,k) = AMAX1(sig(i,k),0.0)
696  ! oct3            w0(i,k)  = beta*w0(i,k)
697  ! oct3   60    continue
698  ! oct3         iflag(i)=4 ! pour version vectorisee
699  ! oct3c convect3         iflag(i)=0
700  ! oct3cccc         return
701  ! oct3       endif
702  ! oct3
703  ! oct3200   continue
704
705  ! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
706
707  DO k = 1, nl
708    DO i = 1, len
709
710      tdif = buoybase(i)
711      ath1 = th(i, 1)
712      ath = th(i, icb(i)-1) - dttrig
713
714      IF (tdif<dtcrit .OR. ath>ath1) THEN
715        sig(i, k) = beta*sig(i, k) - 2.*alpha*tdif*tdif
716        sig(i, k) = amax1(sig(i,k), 0.0)
717        w0(i, k) = beta*w0(i, k)
718        iflag(i) = 4 ! pour version vectorisee
719        ! convect3         iflag(i)=0
720      END IF
721
722    END DO
723  END DO
724
725  ! fin oct3 --
726
727  RETURN
728END SUBROUTINE cv30_trigger
729
730SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
731    plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
732    th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
733    iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
734    v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
735  USE print_control_mod, ONLY: lunout
736  IMPLICIT NONE
737
738  include "cv30param.h"
739
740  ! inputs:
741  INTEGER len, ncum, nd, ntra, nloc
742  INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
743  REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
744  REAL pbase1(len), buoybase1(len)
745  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
746  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
747  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
748  REAL tvp1(len, nd), clw1(len, nd)
749  REAL th1(len, nd)
750  REAL sig1(len, nd), w01(len, nd)
751  REAL tra1(len, nd, ntra)
752
753  ! outputs:
754  ! en fait, on a nloc=len pour l'instant (cf cv_driver)
755  INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
756  REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
757  REAL pbase(nloc), buoybase(nloc)
758  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
759  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
760  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
761  REAL tvp(nloc, nd), clw(nloc, nd)
762  REAL th(nloc, nd)
763  REAL sig(nloc, nd), w0(nloc, nd)
764  REAL tra(nloc, nd, ntra)
765
766  ! local variables:
767  INTEGER i, k, nn, j
768
769  CHARACTER (LEN=20) :: modname = 'cv30_compress'
770  CHARACTER (LEN=80) :: abort_message
771
772
773  DO k = 1, nl + 1
774    nn = 0
775    DO i = 1, len
776      IF (iflag1(i)==0) THEN
777        nn = nn + 1
778        sig(nn, k) = sig1(i, k)
779        w0(nn, k) = w01(i, k)
780        t(nn, k) = t1(i, k)
781        q(nn, k) = q1(i, k)
782        qs(nn, k) = qs1(i, k)
783        u(nn, k) = u1(i, k)
784        v(nn, k) = v1(i, k)
785        gz(nn, k) = gz1(i, k)
786        h(nn, k) = h1(i, k)
787        lv(nn, k) = lv1(i, k)
788        cpn(nn, k) = cpn1(i, k)
789        p(nn, k) = p1(i, k)
790        ph(nn, k) = ph1(i, k)
791        tv(nn, k) = tv1(i, k)
792        tp(nn, k) = tp1(i, k)
793        tvp(nn, k) = tvp1(i, k)
794        clw(nn, k) = clw1(i, k)
795        th(nn, k) = th1(i, k)
796      END IF
797    END DO
798  END DO
799
800  ! do 121 j=1,ntra
801  ! do 111 k=1,nd
802  ! nn=0
803  ! do 101 i=1,len
804  ! if(iflag1(i).eq.0)then
805  ! nn=nn+1
806  ! tra(nn,k,j)=tra1(i,k,j)
807  ! endif
808  ! 101  continue
809  ! 111  continue
810  ! 121  continue
811
812  IF (nn/=ncum) THEN
813    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
814    abort_message = ''
815    CALL abort_physic(modname, abort_message, 1)
816  END IF
817
818  nn = 0
819  DO i = 1, len
820    IF (iflag1(i)==0) THEN
821      nn = nn + 1
822      pbase(nn) = pbase1(i)
823      buoybase(nn) = buoybase1(i)
824      plcl(nn) = plcl1(i)
825      tnk(nn) = tnk1(i)
826      qnk(nn) = qnk1(i)
827      gznk(nn) = gznk1(i)
828      nk(nn) = nk1(i)
829      icb(nn) = icb1(i)
830      icbs(nn) = icbs1(i)
831      iflag(nn) = iflag1(i)
832    END IF
833  END DO
834
835  RETURN
836END SUBROUTINE cv30_compress
837
838SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, &
839    q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, &
840    ep, sigp, buoy)
841    ! epmax_cape: ajout arguments
842  IMPLICIT NONE
843
844  ! ---------------------------------------------------------------------
845  ! Purpose:
846  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
847  ! &
848  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
849  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
850  ! &
851  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
852
853  ! Main differences convect3/convect4:
854  ! - icbs (input) is the first level above LCL (may differ from icb)
855  ! - many minor differences in the iterations
856  ! - condensed water not removed from tvp in convect3
857  ! - vertical profile of buoyancy computed here (use of buoybase)
858  ! - the determination of inb is different
859  ! - no inb1, only inb in output
860  ! ---------------------------------------------------------------------
861
862  include "cvthermo.h"
863  include "cv30param.h"
864  include "conema3.h"
865
866  ! inputs:
867  INTEGER ncum, nd, nloc
868  INTEGER icb(nloc), icbs(nloc), nk(nloc)
869  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
870  REAL p(nloc, nd)
871  REAL tnk(nloc), qnk(nloc), gznk(nloc)
872  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
873  REAL pbase(nloc), buoybase(nloc), plcl(nloc)
874
875  ! outputs:
876  INTEGER inb(nloc)
877  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
878  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
879  REAL buoy(nloc, nd)
880
881  ! local variables:
882  INTEGER i, k
883  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
884  REAL by, defrac, pden
885  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
886  LOGICAL lcape(nloc)
887
888  ! =====================================================================
889  ! --- SOME INITIALIZATIONS
890  ! =====================================================================
891
892  DO k = 1, nl
893    DO i = 1, ncum
894      ep(i, k) = 0.0
895      sigp(i, k) = spfac
896    END DO
897  END DO
898
899  ! =====================================================================
900  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
901  ! =====================================================================
902
903  ! ---       The procedure is to solve the equation.
904  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
905
906  ! ***  Calculate certain parcel quantities, including static energy   ***
907
908
909  DO i = 1, ncum
910    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) & ! debug     &
911                                                  ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
912      +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
913  END DO
914
915
916  ! ***  Find lifted parcel quantities above cloud base    ***
917
918
919  DO k = minorig + 1, nl
920    DO i = 1, ncum
921      ! ori         if(k.ge.(icb(i)+1))then
922      IF (k>=(icbs(i)+1)) THEN ! convect3
923        tg = t(i, k)
924        qg = qs(i, k)
925        ! debug       alv=lv0-clmcpv*(t(i,k)-t0)
926        alv = lv0 - clmcpv*(t(i,k)-273.15)
927
928        ! First iteration.
929
930        ! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
931        s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
932          +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
933        s = 1./s
934        ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
935        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
936        tg = tg + s*(ah0(i)-ahg)
937        ! ori          tg=max(tg,35.0)
938        ! debug        tc=tg-t0
939        tc = tg - 273.15
940        denom = 243.5 + tc
941        denom = max(denom, 1.0) ! convect3
942        ! ori          if(tc.ge.0.0)then
943        es = 6.112*exp(17.67*tc/denom)
944        ! ori          else
945        ! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
946        ! ori          endif
947        qg = eps*es/(p(i,k)-es*(1.-eps))
948
949        ! Second iteration.
950
951        ! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
952        ! ori          s=1./s
953        ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
954        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
955        tg = tg + s*(ah0(i)-ahg)
956        ! ori          tg=max(tg,35.0)
957        ! debug        tc=tg-t0
958        tc = tg - 273.15
959        denom = 243.5 + tc
960        denom = max(denom, 1.0) ! convect3
961        ! ori          if(tc.ge.0.0)then
962        es = 6.112*exp(17.67*tc/denom)
963        ! ori          else
964        ! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
965        ! ori          endif
966        qg = eps*es/(p(i,k)-es*(1.-eps))
967
968        ! debug        alv=lv0-clmcpv*(t(i,k)-t0)
969        alv = lv0 - clmcpv*(t(i,k)-273.15)
970        ! print*,'cpd dans convect2 ',cpd
971        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
972        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
973
974        ! ori c approximation here:
975        ! ori
976        ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
977
978        ! convect3: no approximation:
979        tp(i, k) = (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
980
981        clw(i, k) = qnk(i) - qg
982        clw(i, k) = max(0.0, clw(i,k))
983        rg = qg/(1.-qnk(i))
984        ! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
985        ! convect3: (qg utilise au lieu du vrai mixing ratio rg):
986        tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
987      END IF
988    END DO
989  END DO
990
991  ! =====================================================================
992  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
993  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
994  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
995  ! =====================================================================
996
997  ! ori      do 320 k=minorig+1,nl
998  DO k = 1, nl ! convect3
999    DO i = 1, ncum
1000      pden = ptcrit - pbcrit
1001      ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1002      ep(i, k) = amax1(ep(i,k), 0.0)
1003      ep(i, k) = amin1(ep(i,k), epmax)
1004      sigp(i, k) = spfac
1005      ! ori          if(k.ge.(nk(i)+1))then
1006      ! ori            tca=tp(i,k)-t0
1007      ! ori            if(tca.ge.0.0)then
1008      ! ori              elacrit=elcrit
1009      ! ori            else
1010      ! ori              elacrit=elcrit*(1.0-tca/tlcrit)
1011      ! ori            endif
1012      ! ori            elacrit=max(elacrit,0.0)
1013      ! ori            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
1014      ! ori            ep(i,k)=max(ep(i,k),0.0 )
1015      ! ori            ep(i,k)=min(ep(i,k),1.0 )
1016      ! ori            sigp(i,k)=sigs
1017      ! ori          endif
1018    END DO
1019  END DO
1020
1021  ! =====================================================================
1022  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1023  ! --- VIRTUAL TEMPERATURE
1024  ! =====================================================================
1025
1026  ! dans convect3, tvp est calcule en une seule fois, et sans retirer
1027  ! l'eau condensee (~> reversible CAPE)
1028
1029  ! ori      do 340 k=minorig+1,nl
1030  ! ori        do 330 i=1,ncum
1031  ! ori        if(k.ge.(icb(i)+1))then
1032  ! ori          tvp(i,k)=tvp(i,k)*(1.0-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  ! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1035  ! ori        endif
1036  ! ori 330    continue
1037  ! ori 340  continue
1038
1039  ! ori      do 350 i=1,ncum
1040  ! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1041  ! ori 350  continue
1042
1043  DO i = 1, ncum ! convect3
1044    tp(i, nlp) = tp(i, nl) ! convect3
1045  END DO ! convect3
1046
1047  ! =====================================================================
1048  ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1049  ! =====================================================================
1050
1051  ! -- this is for convect3 only:
1052
1053  ! first estimate of buoyancy:
1054
1055  DO i = 1, ncum
1056    DO k = 1, nl
1057      buoy(i, k) = tvp(i, k) - tv(i, k)
1058    END DO
1059  END DO
1060
1061  ! set buoyancy=buoybase for all levels below base
1062  ! for safety, set buoy(icb)=buoybase
1063
1064  DO i = 1, ncum
1065    DO k = 1, nl
1066      IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1067        buoy(i, k) = buoybase(i)
1068      END IF
1069    END DO
1070    ! IM cf. CRio/JYG 270807   buoy(icb(i),k)=buoybase(i)
1071    buoy(i, icb(i)) = buoybase(i)
1072  END DO
1073
1074  ! -- end convect3
1075
1076  ! =====================================================================
1077  ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1078  ! --- LEVEL OF NEUTRAL BUOYANCY
1079  ! =====================================================================
1080
1081  ! -- this is for convect3 only:
1082
1083  DO i = 1, ncum
1084    inb(i) = nl - 1
1085  END DO
1086
1087  DO i = 1, ncum
1088    DO k = 1, nl - 1
1089      IF ((k>=icb(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1090        inb(i) = min(inb(i), k)
1091      END IF
1092    END DO
1093  END DO
1094
1095  ! -- end convect3
1096
1097  ! ori      do 510 i=1,ncum
1098  ! ori        cape(i)=0.0
1099  ! ori        capem(i)=0.0
1100  ! ori        inb(i)=icb(i)+1
1101  ! ori        inb1(i)=inb(i)
1102  ! ori 510  continue
1103
1104  ! Originial Code
1105
1106  ! do 530 k=minorig+1,nl-1
1107  ! do 520 i=1,ncum
1108  ! if(k.ge.(icb(i)+1))then
1109  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1110  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1111  ! cape(i)=cape(i)+by
1112  ! if(by.ge.0.0)inb1(i)=k+1
1113  ! if(cape(i).gt.0.0)then
1114  ! inb(i)=k+1
1115  ! capem(i)=cape(i)
1116  ! endif
1117  ! endif
1118  ! 520    continue
1119  ! 530  continue
1120  ! do 540 i=1,ncum
1121  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1122  ! cape(i)=capem(i)+byp
1123  ! defrac=capem(i)-cape(i)
1124  ! defrac=max(defrac,0.001)
1125  ! frac(i)=-cape(i)/defrac
1126  ! frac(i)=min(frac(i),1.0)
1127  ! frac(i)=max(frac(i),0.0)
1128  ! 540   continue
1129
1130  ! K Emanuel fix
1131
1132  ! call zilch(byp,ncum)
1133  ! do 530 k=minorig+1,nl-1
1134  ! do 520 i=1,ncum
1135  ! if(k.ge.(icb(i)+1))then
1136  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1137  ! cape(i)=cape(i)+by
1138  ! if(by.ge.0.0)inb1(i)=k+1
1139  ! if(cape(i).gt.0.0)then
1140  ! inb(i)=k+1
1141  ! capem(i)=cape(i)
1142  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1143  ! endif
1144  ! endif
1145  ! 520    continue
1146  ! 530  continue
1147  ! do 540 i=1,ncum
1148  ! inb(i)=max(inb(i),inb1(i))
1149  ! cape(i)=capem(i)+byp(i)
1150  ! defrac=capem(i)-cape(i)
1151  ! defrac=max(defrac,0.001)
1152  ! frac(i)=-cape(i)/defrac
1153  ! frac(i)=min(frac(i),1.0)
1154  ! frac(i)=max(frac(i),0.0)
1155  ! 540   continue
1156
1157  ! J Teixeira fix
1158
1159  ! ori      call zilch(byp,ncum)
1160  ! ori      do 515 i=1,ncum
1161  ! ori        lcape(i)=.true.
1162  ! ori 515  continue
1163  ! ori      do 530 k=minorig+1,nl-1
1164  ! ori        do 520 i=1,ncum
1165  ! ori          if(cape(i).lt.0.0)lcape(i)=.false.
1166  ! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
1167  ! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1168  ! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1169  ! ori            cape(i)=cape(i)+by
1170  ! ori            if(by.ge.0.0)inb1(i)=k+1
1171  ! ori            if(cape(i).gt.0.0)then
1172  ! ori              inb(i)=k+1
1173  ! ori              capem(i)=cape(i)
1174  ! ori            endif
1175  ! ori          endif
1176  ! ori 520    continue
1177  ! ori 530  continue
1178  ! ori      do 540 i=1,ncum
1179  ! ori          cape(i)=capem(i)+byp(i)
1180  ! ori          defrac=capem(i)-cape(i)
1181  ! ori          defrac=max(defrac,0.001)
1182  ! ori          frac(i)=-cape(i)/defrac
1183  ! ori          frac(i)=min(frac(i),1.0)
1184  ! ori          frac(i)=max(frac(i),0.0)
1185  ! ori 540  continue
1186
1187  ! =====================================================================
1188  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1189  ! =====================================================================
1190
1191  ! ym      do i=1,ncum*nlp
1192  ! ym       hp(i,1)=h(i,1)
1193  ! ym      enddo
1194
1195  DO k = 1, nlp
1196    DO i = 1, ncum
1197      hp(i, k) = h(i, k)
1198    END DO
1199  END DO
1200
1201  DO k = minorig + 1, nl
1202    DO i = 1, ncum
1203      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1204        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
1205          )
1206      END IF
1207    END DO
1208  END DO
1209
1210  RETURN
1211END SUBROUTINE cv30_undilute2
1212
1213SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
1214    sig, w0, cape, m)
1215  IMPLICIT NONE
1216
1217  ! ===================================================================
1218  ! ---  CLOSURE OF CONVECT3
1219
1220  ! vectorization: S. Bony
1221  ! ===================================================================
1222
1223  include "cvthermo.h"
1224  include "cv30param.h"
1225
1226  ! input:
1227  INTEGER ncum, nd, nloc
1228  INTEGER icb(nloc), inb(nloc)
1229  REAL pbase(nloc)
1230  REAL p(nloc, nd), ph(nloc, nd+1)
1231  REAL tv(nloc, nd), buoy(nloc, nd)
1232
1233  ! input/output:
1234  REAL sig(nloc, nd), w0(nloc, nd)
1235
1236  ! output:
1237  REAL cape(nloc)
1238  REAL m(nloc, nd)
1239
1240  ! local variables:
1241  INTEGER i, j, k, icbmax
1242  REAL deltap, fac, w, amu
1243  REAL dtmin(nloc, nd), sigold(nloc, nd)
1244
1245  ! -------------------------------------------------------
1246  ! -- Initialization
1247  ! -------------------------------------------------------
1248
1249  DO k = 1, nl
1250    DO i = 1, ncum
1251      m(i, k) = 0.0
1252    END DO
1253  END DO
1254
1255  ! -------------------------------------------------------
1256  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
1257  ! -------------------------------------------------------
1258
1259  ! update sig and w0 above LNB:
1260
1261  DO k = 1, nl - 1
1262    DO i = 1, ncum
1263      IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1264        sig(i, k) = beta*sig(i, k) + 2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb( &
1265          i)))
1266        sig(i, k) = amax1(sig(i,k), 0.0)
1267        w0(i, k) = beta*w0(i, k)
1268      END IF
1269    END DO
1270  END DO
1271
1272  ! compute icbmax:
1273
1274  icbmax = 2
1275  DO i = 1, ncum
1276    icbmax = max(icbmax, icb(i))
1277  END DO
1278
1279  ! update sig and w0 below cloud base:
1280
1281  DO k = 1, icbmax
1282    DO i = 1, ncum
1283      IF (k<=icb(i)) THEN
1284        sig(i, k) = beta*sig(i, k) - 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
1285        sig(i, k) = amax1(sig(i,k), 0.0)
1286        w0(i, k) = beta*w0(i, k)
1287      END IF
1288    END DO
1289  END DO
1290
1291  ! !      if(inb.lt.(nl-1))then
1292  ! !         do 85 i=inb+1,nl-1
1293  ! !            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
1294  ! !     1              abs(buoy(inb))
1295  ! !            sig(i)=amax1(sig(i),0.0)
1296  ! !            w0(i)=beta*w0(i)
1297  ! !   85    continue
1298  ! !      end if
1299
1300  ! !      do 87 i=1,icb
1301  ! !         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
1302  ! !         sig(i)=amax1(sig(i),0.0)
1303  ! !         w0(i)=beta*w0(i)
1304  ! !   87 continue
1305
1306  ! -------------------------------------------------------------
1307  ! -- Reset fractional areas of updrafts and w0 at initial time
1308  ! -- and after 10 time steps of no convection
1309  ! -------------------------------------------------------------
1310
1311  DO k = 1, nl - 1
1312    DO i = 1, ncum
1313      IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1314        sig(i, k) = 0.0
1315        w0(i, k) = 0.0
1316      END IF
1317    END DO
1318  END DO
1319
1320  ! -------------------------------------------------------------
1321  ! -- Calculate convective available potential energy (cape),
1322  ! -- vertical velocity (w), fractional area covered by
1323  ! -- undilute updraft (sig), and updraft mass flux (m)
1324  ! -------------------------------------------------------------
1325
1326  DO i = 1, ncum
1327    cape(i) = 0.0
1328  END DO
1329
1330  ! compute dtmin (minimum buoyancy between ICB and given level k):
1331
1332  DO i = 1, ncum
1333    DO k = 1, nl
1334      dtmin(i, k) = 100.0
1335    END DO
1336  END DO
1337
1338  DO i = 1, ncum
1339    DO k = 1, nl
1340      DO j = minorig, nl
1341        IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
1342            1))) THEN
1343          dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1344        END IF
1345      END DO
1346    END DO
1347  END DO
1348
1349  ! the interval on which cape is computed starts at pbase :
1350  DO k = 1, nl
1351    DO i = 1, ncum
1352
1353      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1354
1355        deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1356        cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1)
1357        cape(i) = amax1(0.0, cape(i))
1358        sigold(i, k) = sig(i, k)
1359
1360        ! dtmin(i,k)=100.0
1361        ! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1362        ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1363        ! 97     continue
1364
1365        sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
1366        sig(i, k) = amax1(sig(i,k), 0.0)
1367        sig(i, k) = amin1(sig(i,k), 0.01)
1368        fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1369        w = (1.-beta)*fac*sqrt(cape(i)) + beta*w0(i, k)
1370        amu = 0.5*(sig(i,k)+sigold(i,k))*w
1371        m(i, k) = amu*0.007*p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1372        w0(i, k) = w
1373      END IF
1374
1375    END DO
1376  END DO
1377
1378  DO i = 1, ncum
1379    w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1380    m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/ &
1381      (ph(i,icb(i)+1)-ph(i,icb(i)+2))
1382    sig(i, icb(i)) = sig(i, icb(i)+1)
1383    sig(i, icb(i)-1) = sig(i, icb(i))
1384  END DO
1385
1386
1387  ! !      cape=0.0
1388  ! !      do 98 i=icb+1,inb
1389  ! !         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1390  ! !         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
1391  ! !         dcape=rrd*buoy(i-1)*deltap/p(i-1)
1392  ! !         dlnp=deltap/p(i-1)
1393  ! !         cape=amax1(0.0,cape)
1394  ! !         sigold=sig(i)
1395
1396  ! !         dtmin=100.0
1397  ! !         do 97 j=icb,i-1
1398  ! !            dtmin=amin1(dtmin,buoy(j))
1399  ! !   97    continue
1400
1401  ! !         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
1402  ! !         sig(i)=amax1(sig(i),0.0)
1403  ! !         sig(i)=amin1(sig(i),0.01)
1404  ! !         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1405  ! !         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
1406  ! !         amu=0.5*(sig(i)+sigold)*w
1407  ! !         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
1408  ! !         w0(i)=w
1409  ! !   98 continue
1410  ! !      w0(icb)=0.5*w0(icb+1)
1411  ! !      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1412  ! !      sig(icb)=sig(icb+1)
1413  ! !      sig(icb-1)=sig(icb)
1414
1415  RETURN
1416END SUBROUTINE cv30_closure
1417
1418SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
1419    u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
1420    vent, sij, elij, ments, qents, traent)
1421  IMPLICIT NONE
1422
1423  ! ---------------------------------------------------------------------
1424  ! a faire:
1425  ! - changer rr(il,1) -> qnk(il)
1426  ! - vectorisation de la partie normalisation des flux (do 789...)
1427  ! ---------------------------------------------------------------------
1428
1429  include "cvthermo.h"
1430  include "cv30param.h"
1431
1432  ! inputs:
1433  INTEGER ncum, nd, na, ntra, nloc
1434  INTEGER icb(nloc), inb(nloc), nk(nloc)
1435  REAL sig(nloc, nd)
1436  REAL qnk(nloc)
1437  REAL ph(nloc, nd+1)
1438  REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1439  REAL u(nloc, nd), v(nloc, nd)
1440  REAL tra(nloc, nd, ntra) ! input of convect3
1441  REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
1442  REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
1443  REAL m(nloc, na) ! input of convect3
1444
1445  ! outputs:
1446  REAL ment(nloc, na, na), qent(nloc, na, na)
1447  REAL uent(nloc, na, na), vent(nloc, na, na)
1448  REAL sij(nloc, na, na), elij(nloc, na, na)
1449  REAL traent(nloc, nd, nd, ntra)
1450  REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
1451  REAL sigij(nloc, nd, nd)
1452
1453  ! local variables:
1454  INTEGER i, j, k, il, im, jm
1455  INTEGER num1, num2
1456  INTEGER nent(nloc, na)
1457  REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1458  REAL alt, smid, sjmin, sjmax, delp, delm
1459  REAL asij(nloc), smax(nloc), scrit(nloc)
1460  REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1461  REAL wgh
1462  REAL zm(nloc, na)
1463  LOGICAL lwork(nloc)
1464
1465  ! =====================================================================
1466  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1467  ! =====================================================================
1468
1469  ! ori        do 360 i=1,ncum*nlp
1470  DO j = 1, nl
1471    DO i = 1, ncum
1472      nent(i, j) = 0
1473      ! in convect3, m is computed in cv3_closure
1474      ! ori          m(i,1)=0.0
1475    END DO
1476  END DO
1477
1478  ! ori      do 400 k=1,nlp
1479  ! ori       do 390 j=1,nlp
1480  DO j = 1, nl
1481    DO k = 1, nl
1482      DO i = 1, ncum
1483        qent(i, k, j) = rr(i, j)
1484        uent(i, k, j) = u(i, j)
1485        vent(i, k, j) = v(i, j)
1486        elij(i, k, j) = 0.0
1487        ! ym            ment(i,k,j)=0.0
1488        ! ym            sij(i,k,j)=0.0
1489      END DO
1490    END DO
1491  END DO
1492
1493  ! ym
1494  ment(1:ncum, 1:nd, 1:nd) = 0.0
1495  sij(1:ncum, 1:nd, 1:nd) = 0.0
1496
1497  ! do k=1,ntra
1498  ! do j=1,nd  ! instead nlp
1499  ! do i=1,nd ! instead nlp
1500  ! do il=1,ncum
1501  ! traent(il,i,j,k)=tra(il,j,k)
1502  ! enddo
1503  ! enddo
1504  ! enddo
1505  ! enddo
1506  zm(:, :) = 0.
1507
1508  ! =====================================================================
1509  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1510  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1511  ! --- FRACTION (sij)
1512  ! =====================================================================
1513
1514  DO i = minorig + 1, nl
1515
1516    DO j = minorig, nl
1517      DO il = 1, ncum
1518        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
1519            1)) .AND. (j<=inb(il))) THEN
1520
1521          rti = rr(il, 1) - ep(il, i)*clw(il, i)
1522          bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
1523          anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
1524          denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
1525          dei = denom
1526          IF (abs(dei)<0.01) dei = 0.01
1527          sij(il, i, j) = anum/dei
1528          sij(il, i, i) = 1.0
1529          altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
1530          altem = altem/bf2
1531          cwat = clw(il, j)*(1.-ep(il,j))
1532          stemp = sij(il, i, j)
1533          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1534            anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
1535            denom = denom + lv(il, j)*(rr(il,i)-rti)
1536            IF (abs(denom)<0.01) denom = 0.01
1537            sij(il, i, j) = anum/denom
1538            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - &
1539              rs(il, j)
1540            altem = altem - (bf2-1.)*cwat
1541          END IF
1542          IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1543            qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti
1544            uent(il, i, j) = sij(il, i, j)*u(il, i) + &
1545              (1.-sij(il,i,j))*u(il, nk(il))
1546            vent(il, i, j) = sij(il, i, j)*v(il, i) + &
1547              (1.-sij(il,i,j))*v(il, nk(il))
1548            ! !!!      do k=1,ntra
1549            ! !!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1550            ! !!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
1551            ! !!!      end do
1552            elij(il, i, j) = altem
1553            elij(il, i, j) = amax1(0.0, elij(il,i,j))
1554            ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
1555            nent(il, i) = nent(il, i) + 1
1556          END IF
1557          sij(il, i, j) = amax1(0.0, sij(il,i,j))
1558          sij(il, i, j) = amin1(1.0, sij(il,i,j))
1559        END IF ! new
1560      END DO
1561    END DO
1562
1563    ! do k=1,ntra
1564    ! do j=minorig,nl
1565    ! do il=1,ncum
1566    ! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1567    ! :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1568    ! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1569    ! :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
1570    ! endif
1571    ! enddo
1572    ! enddo
1573    ! enddo
1574
1575
1576    ! ***   if no air can entrain at level i assume that updraft detrains
1577    ! ***
1578    ! ***   at that level and calculate detrained air flux and properties
1579    ! ***
1580
1581
1582    ! @      do 170 i=icb(il),inb(il)
1583
1584    DO il = 1, ncum
1585      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
1586        ! @      if(nent(il,i).eq.0)then
1587        ment(il, i, i) = m(il, i)
1588        qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i)
1589        uent(il, i, i) = u(il, nk(il))
1590        vent(il, i, i) = v(il, nk(il))
1591        elij(il, i, i) = clw(il, i)
1592        ! MAF      sij(il,i,i)=1.0
1593        sij(il, i, i) = 0.0
1594      END IF
1595    END DO
1596  END DO
1597
1598  ! do j=1,ntra
1599  ! do i=minorig+1,nl
1600  ! do il=1,ncum
1601  ! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1602  ! traent(il,i,i,j)=tra(il,nk(il),j)
1603  ! endif
1604  ! enddo
1605  ! enddo
1606  ! enddo
1607
1608  DO j = minorig, nl
1609    DO i = minorig, nl
1610      DO il = 1, ncum
1611        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
1612            inb(il))) THEN
1613          sigij(il, i, j) = sij(il, i, j)
1614        END IF
1615      END DO
1616    END DO
1617  END DO
1618  ! @      enddo
1619
1620  ! @170   continue
1621
1622  ! =====================================================================
1623  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
1624  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
1625  ! =====================================================================
1626
1627  ! ym      call zilch(asum,ncum*nd)
1628  ! ym      call zilch(bsum,ncum*nd)
1629  ! ym      call zilch(csum,ncum*nd)
1630  CALL zilch(asum, nloc*nd)
1631  CALL zilch(csum, nloc*nd)
1632  CALL zilch(csum, nloc*nd)
1633
1634  DO il = 1, ncum
1635    lwork(il) = .FALSE.
1636  END DO
1637
1638  DO i = minorig + 1, nl
1639
1640    num1 = 0
1641    DO il = 1, ncum
1642      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
1643    END DO
1644    IF (num1<=0) GO TO 789
1645
1646
1647    DO il = 1, ncum
1648      IF (i>=icb(il) .AND. i<=inb(il)) THEN
1649        lwork(il) = (nent(il,i)/=0)
1650        qp = rr(il, 1) - ep(il, i)*clw(il, i)
1651        anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
1652          (cpv-cpd)*t(il, i)*(qp-rr(il,i))
1653        denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
1654          (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
1655        IF (abs(denom)<0.01) denom = 0.01
1656        scrit(il) = anum/denom
1657        alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
1658        IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
1659        smax(il) = 0.0
1660        asij(il) = 0.0
1661      END IF
1662    END DO
1663
1664    DO j = nl, minorig, -1
1665
1666      num2 = 0
1667      DO il = 1, ncum
1668        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1669          il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
1670      END DO
1671      IF (num2<=0) GO TO 175
1672
1673      DO il = 1, ncum
1674        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1675            il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
1676
1677          IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
1678            wgh = 1.0
1679            IF (j>i) THEN
1680              sjmax = amax1(sij(il,i,j+1), smax(il))
1681              sjmax = amin1(sjmax, scrit(il))
1682              smax(il) = amax1(sij(il,i,j), smax(il))
1683              sjmin = amax1(sij(il,i,j-1), smax(il))
1684              sjmin = amin1(sjmin, scrit(il))
1685              IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
1686              smid = amin1(sij(il,i,j), scrit(il))
1687            ELSE
1688              sjmax = amax1(sij(il,i,j+1), scrit(il))
1689              smid = amax1(sij(il,i,j), scrit(il))
1690              sjmin = 0.0
1691              IF (j>1) sjmin = sij(il, i, j-1)
1692              sjmin = amax1(sjmin, scrit(il))
1693            END IF
1694            delp = abs(sjmax-smid)
1695            delm = abs(sjmin-smid)
1696            asij(il) = asij(il) + wgh*(delp+delm)
1697            ment(il, i, j) = ment(il, i, j)*(delp+delm)*wgh
1698          END IF
1699        END IF
1700      END DO
1701
1702175 END DO
1703
1704    DO il = 1, ncum
1705      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1706        asij(il) = amax1(1.0E-16, asij(il))
1707        asij(il) = 1.0/asij(il)
1708        asum(il, i) = 0.0
1709        bsum(il, i) = 0.0
1710        csum(il, i) = 0.0
1711      END IF
1712    END DO
1713
1714    DO j = minorig, nl
1715      DO il = 1, ncum
1716        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1717            il)-1) .AND. j<=inb(il)) THEN
1718          ment(il, i, j) = ment(il, i, j)*asij(il)
1719        END IF
1720      END DO
1721    END DO
1722
1723    DO j = minorig, nl
1724      DO il = 1, ncum
1725        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1726            il)-1) .AND. j<=inb(il)) THEN
1727          asum(il, i) = asum(il, i) + ment(il, i, j)
1728          ment(il, i, j) = ment(il, i, j)*sig(il, j)
1729          bsum(il, i) = bsum(il, i) + ment(il, i, j)
1730        END IF
1731      END DO
1732    END DO
1733
1734    DO il = 1, ncum
1735      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1736        bsum(il, i) = amax1(bsum(il,i), 1.0E-16)
1737        bsum(il, i) = 1.0/bsum(il, i)
1738      END IF
1739    END DO
1740
1741    DO j = minorig, nl
1742      DO il = 1, ncum
1743        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1744            il)-1) .AND. j<=inb(il)) THEN
1745          ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i)
1746        END IF
1747      END DO
1748    END DO
1749
1750    DO j = minorig, nl
1751      DO il = 1, ncum
1752        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1753            il)-1) .AND. j<=inb(il)) THEN
1754          csum(il, i) = csum(il, i) + ment(il, i, j)
1755        END IF
1756      END DO
1757    END DO
1758
1759    DO il = 1, ncum
1760      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1761          csum(il,i)<m(il,i)) THEN
1762        nent(il, i) = 0
1763        ment(il, i, i) = m(il, i)
1764        qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i)
1765        uent(il, i, i) = u(il, nk(il))
1766        vent(il, i, i) = v(il, nk(il))
1767        elij(il, i, i) = clw(il, i)
1768        ! MAF        sij(il,i,i)=1.0
1769        sij(il, i, i) = 0.0
1770      END IF
1771    END DO ! il
1772
1773    ! do j=1,ntra
1774    ! do il=1,ncum
1775    ! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1776    ! :     .and. csum(il,i).lt.m(il,i) ) then
1777    ! traent(il,i,i,j)=tra(il,nk(il),j)
1778    ! endif
1779    ! enddo
1780    ! enddo
1781789 END DO
1782
1783  ! MAF: renormalisation de MENT
1784  DO jm = 1, nd
1785    DO im = 1, nd
1786      DO il = 1, ncum
1787        zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
1788      END DO
1789    END DO
1790  END DO
1791
1792  DO jm = 1, nd
1793    DO im = 1, nd
1794      DO il = 1, ncum
1795        IF (zm(il,im)/=0.) THEN
1796          ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
1797        END IF
1798      END DO
1799    END DO
1800  END DO
1801
1802  DO jm = 1, nd
1803    DO im = 1, nd
1804      DO il = 1, ncum
1805        qents(il, im, jm) = qent(il, im, jm)
1806        ments(il, im, jm) = ment(il, im, jm)
1807      END DO
1808    END DO
1809  END DO
1810
1811  RETURN
1812END SUBROUTINE cv30_mixing
1813
1814
1815SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
1816    v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
1817    mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
1818    , wdtraina, wdtrainm) ! 26/08/10  RomP-jyg
1819  IMPLICIT NONE
1820
1821
1822  include "cvthermo.h"
1823  include "cv30param.h"
1824  include "cvflag.h"
1825
1826  ! inputs:
1827  INTEGER ncum, nd, na, ntra, nloc
1828  INTEGER icb(nloc), inb(nloc)
1829  REAL delt, plcl(nloc)
1830  REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1831  REAL u(nloc, nd), v(nloc, nd)
1832  REAL tra(nloc, nd, ntra)
1833  REAL p(nloc, nd), ph(nloc, nd+1)
1834  REAL th(nloc, na), gz(nloc, na)
1835  REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
1836  REAL cpn(nloc, na), tv(nloc, na)
1837  REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
1838
1839  ! outputs:
1840  REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
1841  REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
1842  REAL trap(nloc, na, ntra)
1843  REAL b(nloc, na)
1844  ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
1845  ! lascendance adiabatique et des flux melanges Pa et Pm.
1846  ! Distinction des wdtrain
1847  ! Pa = wdtrainA     Pm = wdtrainM
1848  REAL wdtraina(nloc, na), wdtrainm(nloc, na)
1849
1850  ! local variables
1851  INTEGER i, j, k, il, num1
1852  REAL tinv, delti
1853  REAL awat, afac, afac1, afac2, bfac
1854  REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
1855  REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1856  REAL ampmax
1857  REAL lvcp(nloc, na)
1858  REAL wdtrain(nloc)
1859  LOGICAL lwork(nloc)
1860
1861
1862  ! ------------------------------------------------------
1863
1864  delti = 1./delt
1865  tinv = 1./3.
1866
1867  mp(:, :) = 0.
1868
1869  DO i = 1, nl
1870    DO il = 1, ncum
1871      mp(il, i) = 0.0
1872      rp(il, i) = rr(il, i)
1873      up(il, i) = u(il, i)
1874      vp(il, i) = v(il, i)
1875      wt(il, i) = 0.001
1876      water(il, i) = 0.0
1877      evap(il, i) = 0.0
1878      b(il, i) = 0.0
1879      lvcp(il, i) = lv(il, i)/cpn(il, i)
1880    END DO
1881  END DO
1882
1883  ! do k=1,ntra
1884  ! do i=1,nd
1885  ! do il=1,ncum
1886  ! trap(il,i,k)=tra(il,i,k)
1887  ! enddo
1888  ! enddo
1889  ! enddo
1890  ! ! RomP >>>
1891  DO i = 1, nd
1892    DO il = 1, ncum
1893      wdtraina(il, i) = 0.0
1894      wdtrainm(il, i) = 0.0
1895    END DO
1896  END DO
1897  ! ! RomP <<<
1898
1899  ! ***  check whether ep(inb)=0, if so, skip precipitating    ***
1900  ! ***             downdraft calculation                      ***
1901
1902
1903  DO il = 1, ncum
1904    lwork(il) = .TRUE.
1905    IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE.
1906  END DO
1907
1908  CALL zilch(wdtrain, ncum)
1909
1910  DO i = nl + 1, 1, -1
1911
1912    num1 = 0
1913    DO il = 1, ncum
1914      IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
1915    END DO
1916    IF (num1<=0) GO TO 400
1917
1918
1919    ! ***  integrate liquid water equation to find condensed water   ***
1920    ! ***                and condensed water flux                    ***
1921
1922
1923
1924    ! ***                    begin downdraft loop                    ***
1925
1926
1927
1928    ! ***              calculate detrained precipitation             ***
1929
1930    DO il = 1, ncum
1931      IF (i<=inb(il) .AND. lwork(il)) THEN
1932        IF (cvflag_grav) THEN
1933          wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
1934          wdtraina(il, i) = wdtrain(il)/grav !   Pa  26/08/10   RomP
1935        ELSE
1936          wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
1937          wdtraina(il, i) = wdtrain(il)/10. !   Pa  26/08/10   RomP
1938        END IF
1939      END IF
1940    END DO
1941
1942    IF (i>1) THEN
1943
1944      DO j = 1, i - 1
1945        DO il = 1, ncum
1946          IF (i<=inb(il) .AND. lwork(il)) THEN
1947            awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
1948            awat = amax1(awat, 0.0)
1949            IF (cvflag_grav) THEN
1950              wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
1951            ELSE
1952              wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
1953            END IF
1954          END IF
1955        END DO
1956      END DO
1957      DO il = 1, ncum
1958        IF (cvflag_grav) THEN
1959          wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) !   Pm  26/08/10   RomP
1960        ELSE
1961          wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) !   Pm  26/08/10   RomP
1962        END IF
1963      END DO
1964
1965    END IF
1966
1967
1968    ! ***    find rain water and evaporation using provisional   ***
1969    ! ***              estimates of rp(i)and rp(i-1)             ***
1970
1971
1972    DO il = 1, ncum
1973
1974      IF (i<=inb(il) .AND. lwork(il)) THEN
1975
1976        wt(il, i) = 45.0
1977
1978        IF (i<inb(il)) THEN
1979          rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
1980            i))+gz(il,i+1)-gz(il,i))/lv(il, i)
1981          rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
1982        END IF
1983        rp(il, i) = amax1(rp(il,i), 0.0)
1984        rp(il, i) = amin1(rp(il,i), rs(il,i))
1985        rp(il, inb(il)) = rr(il, inb(il))
1986
1987        IF (i==1) THEN
1988          afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
1989        ELSE
1990          rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
1991            i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
1992          rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
1993          rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
1994          rp(il, i-1) = amax1(rp(il,i-1), 0.0)
1995          afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i) &
1996            )
1997          afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
1998            (1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
1999          afac = 0.5*(afac1+afac2)
2000        END IF
2001        IF (i==inb(il)) afac = 0.0
2002        afac = amax1(afac, 0.0)
2003        bfac = 1./(sigd*wt(il,i))
2004
2005        ! jyg1
2006        ! cc        sigt=1.0
2007        ! cc        if(i.ge.icb)sigt=sigp(i)
2008        ! prise en compte de la variation progressive de sigt dans
2009        ! les couches icb et icb-1:
2010        ! pour plcl<ph(i+1), pr1=0 & pr2=1
2011        ! pour plcl>ph(i),   pr1=1 & pr2=0
2012        ! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2013        ! sur le nuage, et pr2 est la proportion sous la base du
2014        ! nuage.
2015        pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2016        pr1 = max(0., min(1.,pr1))
2017        pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2018        pr2 = max(0., min(1.,pr2))
2019        sigt = sigp(il, i)*pr1 + pr2
2020        ! jyg2
2021
2022        b6 = bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac
2023        c6 = water(il, i+1) + bfac*wdtrain(il) - 50.*sigd*bfac*(ph(il,i)-ph( &
2024          il,i+1))*evap(il, i+1)
2025        IF (c6>0.0) THEN
2026          revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
2027          evap(il, i) = sigt*afac*revap
2028          water(il, i) = revap*revap
2029        ELSE
2030          evap(il, i) = -evap(il, i+1) + 0.02*(wdtrain(il)+sigd*wt(il,i)* &
2031            water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1)))
2032        END IF
2033
2034        ! ***  calculate precipitating downdraft mass flux under     ***
2035        ! ***              hydrostatic approximation                 ***
2036
2037        IF (i/=1) THEN
2038
2039          tevap = amax1(0.0, evap(il,i))
2040          delth = amax1(0.001, (th(il,i)-th(il,i-1)))
2041          IF (cvflag_grav) THEN
2042            mp(il, i) = 100.*ginv*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/ &
2043              delth
2044          ELSE
2045            mp(il, i) = 10.*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth
2046          END IF
2047
2048          ! ***           if hydrostatic assumption fails,             ***
2049          ! ***   solve cubic difference equation for downdraft theta  ***
2050          ! ***  and mass flux from two simultaneous differential eqns ***
2051
2052          amfac = sigd*sigd*70.0*ph(il, i)*(p(il,i-1)-p(il,i))* &
2053            (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2054          amp2 = abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
2055          IF (amp2>(0.1*amfac)) THEN
2056            xf = 100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1))
2057            tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i)* &
2058              sigd*th(il,i))
2059            af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv
2060            bf = 2.*(tinv*mp(il,i+1))**3 + tinv*mp(il, i+1)*xf*tf + &
2061              50.*(p(il,i-1)-p(il,i))*xf*tevap
2062            fac2 = 1.0
2063            IF (bf<0.0) fac2 = -1.0
2064            bf = abs(bf)
2065            ur = 0.25*bf*bf - af*af*af*tinv*tinv*tinv
2066            IF (ur>=0.0) THEN
2067              sru = sqrt(ur)
2068              fac = 1.0
2069              IF ((0.5*bf-sru)<0.0) fac = -1.0
2070              mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
2071                fac*(abs(0.5*bf-sru))**tinv
2072            ELSE
2073              d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2074              IF (fac2<0.0) d = 3.14159 - d
2075              mp(il, i) = mp(il, i+1)*tinv + 2.*sqrt(af*tinv)*cos(d*tinv)
2076            END IF
2077            mp(il, i) = amax1(0.0, mp(il,i))
2078
2079            IF (cvflag_grav) THEN
2080              ! jyg : il y a vraisemblablement une erreur dans la ligne 2
2081              ! suivante:
2082              ! il faut diviser par (mp(il,i)*sigd*grav) et non par
2083              ! (mp(il,i)+sigd*0.1).
2084              ! Et il faut bien revoir les facteurs 100.
2085              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, &
2086                i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2087                )*sigd*th(il,i))
2088            ELSE
2089              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, &
2090                i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2091                )*sigd*th(il,i))
2092            END IF
2093            b(il, i-1) = amax1(b(il,i-1), 0.0)
2094          END IF
2095
2096          ! ***         limit magnitude of mp(i) to meet cfl condition
2097          ! ***
2098
2099          ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
2100          amp2 = 2.0*(ph(il,i-1)-ph(il,i))*delti
2101          ampmax = amin1(ampmax, amp2)
2102          mp(il, i) = amin1(mp(il,i), ampmax)
2103
2104          ! ***      force mp to decrease linearly to zero
2105          ! ***
2106          ! ***       between cloud base and the surface
2107          ! ***
2108
2109          IF (p(il,i)>p(il,icb(il))) THEN
2110            mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ &
2111              (p(il,1)-p(il,icb(il)))
2112          END IF
2113
2114        END IF ! i.eq.1
2115
2116        ! ***       find mixing ratio of precipitating downdraft     ***
2117
2118
2119        IF (i/=inb(il)) THEN
2120
2121          rp(il, i) = rr(il, i)
2122
2123          IF (mp(il,i)>mp(il,i+1)) THEN
2124
2125            IF (cvflag_grav) THEN
2126              rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2127                rr(il, i)*(mp(il,i)-mp(il,i+1)) + 100.*ginv*0.5*sigd*(ph(il,i &
2128                )-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2129            ELSE
2130              rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2131                rr(il, i)*(mp(il,i)-mp(il,i+1)) + 5.*sigd*(ph(il,i)-ph(il,i+1 &
2132                ))*(evap(il,i+1)+evap(il,i))
2133            END IF
2134            rp(il, i) = rp(il, i)/mp(il, i)
2135            up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+ &
2136              1))
2137            up(il, i) = up(il, i)/mp(il, i)
2138            vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+ &
2139              1))
2140            vp(il, i) = vp(il, i)/mp(il, i)
2141
2142            ! do j=1,ntra
2143            ! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2144            ! testmaf     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2145            ! :            +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2146            ! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2147            ! end do
2148
2149          ELSE
2150
2151            IF (mp(il,i+1)>1.0E-16) THEN
2152              IF (cvflag_grav) THEN
2153                rp(il, i) = rp(il, i+1) + 100.*ginv*0.5*sigd*(ph(il,i)-ph(il, &
2154                  i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2155              ELSE
2156                rp(il, i) = rp(il, i+1) + 5.*sigd*(ph(il,i)-ph(il,i+1))*(evap &
2157                  (il,i+1)+evap(il,i))/mp(il, i+1)
2158              END IF
2159              up(il, i) = up(il, i+1)
2160              vp(il, i) = vp(il, i+1)
2161
2162              ! do j=1,ntra
2163              ! trap(il,i,j)=trap(il,i+1,j)
2164              ! end do
2165
2166            END IF
2167          END IF
2168          rp(il, i) = amin1(rp(il,i), rs(il,i))
2169          rp(il, i) = amax1(rp(il,i), 0.0)
2170
2171        END IF
2172      END IF
2173    END DO
2174
2175400 END DO
2176
2177  RETURN
2178END SUBROUTINE cv30_unsat
2179
2180SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
2181    tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
2182    wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
2183    tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
2184    mike, tls, tps, qcondc, wd)
2185  IMPLICIT NONE
2186
2187  include "cvthermo.h"
2188  include "cv30param.h"
2189  include "cvflag.h"
2190  include "conema3.h"
2191
2192  ! inputs:
2193  INTEGER ncum, nd, na, ntra, nloc
2194  INTEGER icb(nloc), inb(nloc)
2195  REAL delt
2196  REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
2197  REAL tra(nloc, nd, ntra), sig(nloc, nd)
2198  REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
2199  REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
2200  REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
2201  REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
2202  REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
2203  REAL water(nloc, na), evap(nloc, na), b(nloc, na)
2204  REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
2205  ! ym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2206  REAL vent(nloc, na, na), elij(nloc, na, na)
2207  INTEGER nent(nloc, na)
2208  REAL traent(nloc, na, na, ntra)
2209  REAL tv(nloc, nd), tvp(nloc, nd)
2210
2211  ! input/output:
2212  INTEGER iflag(nloc)
2213
2214  ! outputs:
2215  REAL precip(nloc)
2216  REAL vprecip(nloc, nd+1)
2217  REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
2218  REAL ftra(nloc, nd, ntra)
2219  REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
2220  REAL dnwd0(nloc, nd), mike(nloc, nd)
2221  REAL tls(nloc, nd), tps(nloc, nd)
2222  REAL qcondc(nloc, nd) ! cld
2223  REAL wd(nloc) ! gust
2224
2225  ! local variables:
2226  INTEGER i, k, il, n, j, num1
2227  REAL rat, awat, delti
2228  REAL ax, bx, cx, dx, ex
2229  REAL cpinv, rdcp, dpinv
2230  REAL lvcp(nloc, na), mke(nloc, na)
2231  REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
2232  ! !!      real up1(nloc), dn1(nloc)
2233  REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
2234  REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2235  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
2236  REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
2237
2238
2239  ! -------------------------------------------------------------
2240
2241  ! initialization:
2242
2243  delti = 1.0/delt
2244
2245  DO il = 1, ncum
2246    precip(il) = 0.0
2247    wd(il) = 0.0 ! gust
2248    vprecip(il, nd+1) = 0.
2249  END DO
2250
2251  DO i = 1, nd
2252    DO il = 1, ncum
2253      vprecip(il, i) = 0.0
2254      ft(il, i) = 0.0
2255      fr(il, i) = 0.0
2256      fu(il, i) = 0.0
2257      fv(il, i) = 0.0
2258      qcondc(il, i) = 0.0 ! cld
2259      qcond(il, i) = 0.0 ! cld
2260      nqcond(il, i) = 0.0 ! cld
2261    END DO
2262  END DO
2263
2264  ! do j=1,ntra
2265  ! do i=1,nd
2266  ! do il=1,ncum
2267  ! ftra(il,i,j)=0.0
2268  ! enddo
2269  ! enddo
2270  ! enddo
2271
2272  DO i = 1, nl
2273    DO il = 1, ncum
2274      lvcp(il, i) = lv(il, i)/cpn(il, i)
2275    END DO
2276  END DO
2277
2278
2279
2280  ! ***  calculate surface precipitation in mm/day     ***
2281
2282  DO il = 1, ncum
2283    IF (ep(il,inb(il))>=0.0001) THEN
2284      IF (cvflag_grav) THEN
2285        precip(il) = wt(il, 1)*sigd*water(il, 1)*86400.*1000./(rowl*grav)
2286      ELSE
2287        precip(il) = wt(il, 1)*sigd*water(il, 1)*8640.
2288      END IF
2289    END IF
2290  END DO
2291
2292  ! ***  CALCULATE VERTICAL PROFILE OF  PRECIPITATIONs IN kg/m2/s  ===
2293
2294  ! MAF rajout pour lessivage
2295  DO k = 1, nl
2296    DO il = 1, ncum
2297      IF (k<=inb(il)) THEN
2298        IF (cvflag_grav) THEN
2299          vprecip(il, k) = wt(il, k)*sigd*water(il, k)/grav
2300        ELSE
2301          vprecip(il, k) = wt(il, k)*sigd*water(il, k)/10.
2302        END IF
2303      END IF
2304    END DO
2305  END DO
2306
2307
2308  ! ***  Calculate downdraft velocity scale    ***
2309  ! ***  NE PAS UTILISER POUR L'INSTANT ***
2310
2311  ! !      do il=1,ncum
2312  ! !        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
2313  ! !     :                                  /(sigd*p(il,icb(il)))
2314  ! !      enddo
2315
2316
2317  ! ***  calculate tendencies of lowest level potential temperature  ***
2318  ! ***                      and mixing ratio                        ***
2319
2320  DO il = 1, ncum
2321    work(il) = 1.0/(ph(il,1)-ph(il,2))
2322    am(il) = 0.0
2323  END DO
2324
2325  DO k = 2, nl
2326    DO il = 1, ncum
2327      IF (k<=inb(il)) THEN
2328        am(il) = am(il) + m(il, k)
2329      END IF
2330    END DO
2331  END DO
2332
2333  DO il = 1, ncum
2334
2335    ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
2336    IF (cvflag_grav) THEN
2337      IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
2338      ft(il, 1) = 0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2339        1))/cpn(il,1))
2340    ELSE
2341      IF ((0.1*work(il)*am(il))>=delti) iflag(il) = 1 !consistency vect
2342      ft(il, 1) = 0.1*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2343        1))/cpn(il,1))
2344    END IF
2345
2346    ft(il, 1) = ft(il, 1) - 0.5*lvcp(il, 1)*sigd*(evap(il,1)+evap(il,2))
2347
2348    IF (cvflag_grav) THEN
2349      ft(il, 1) = ft(il, 1) - 0.009*grav*sigd*mp(il, 2)*t(il, 1)*b(il, 1)* &
2350        work(il)
2351    ELSE
2352      ft(il, 1) = ft(il, 1) - 0.09*sigd*mp(il, 2)*t(il, 1)*b(il, 1)*work(il)
2353    END IF
2354
2355    ft(il, 1) = ft(il, 1) + 0.01*sigd*wt(il, 1)*(cl-cpd)*water(il, 2)*(t(il,2 &
2356      )-t(il,1))*work(il)/cpn(il, 1)
2357
2358    IF (cvflag_grav) THEN
2359      ! jyg1  Correction pour mieux conserver l'eau (conformite avec
2360      ! CONVECT4.3)
2361      ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
2362      ! evap)
2363      fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + &
2364        sigd*0.5*(evap(il,1)+evap(il,2))
2365      ! +tard     :          +sigd*evap(il,1)
2366
2367      fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
2368
2369      fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il, &
2370        1))+am(il)*(u(il,2)-u(il,1)))
2371      fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
2372        1))+am(il)*(v(il,2)-v(il,1)))
2373    ELSE ! cvflag_grav
2374      fr(il, 1) = 0.1*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + &
2375        sigd*0.5*(evap(il,1)+evap(il,2))
2376      fr(il, 1) = fr(il, 1) + 0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
2377      fu(il, 1) = fu(il, 1) + 0.1*work(il)*(mp(il,2)*(up(il,2)-u(il, &
2378        1))+am(il)*(u(il,2)-u(il,1)))
2379      fv(il, 1) = fv(il, 1) + 0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
2380        1))+am(il)*(v(il,2)-v(il,1)))
2381    END IF ! cvflag_grav
2382
2383  END DO ! il
2384
2385  ! do j=1,ntra
2386  ! do il=1,ncum
2387  ! if (cvflag_grav) then
2388  ! ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
2389  ! :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2390  ! :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
2391  ! else
2392  ! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2393  ! :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2394  ! :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
2395  ! endif
2396  ! enddo
2397  ! enddo
2398
2399  DO j = 2, nl
2400    DO il = 1, ncum
2401      IF (j<=inb(il)) THEN
2402        IF (cvflag_grav) THEN
2403          fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il, &
2404            j,1)-rr(il,1))
2405          fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il, &
2406            j,1)-u(il,1))
2407          fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il, &
2408            j,1)-v(il,1))
2409        ELSE ! cvflag_grav
2410          fr(il, 1) = fr(il, 1) + 0.1*work(il)*ment(il, j, 1)*(qent(il,j,1)- &
2411            rr(il,1))
2412          fu(il, 1) = fu(il, 1) + 0.1*work(il)*ment(il, j, 1)*(uent(il,j,1)-u &
2413            (il,1))
2414          fv(il, 1) = fv(il, 1) + 0.1*work(il)*ment(il, j, 1)*(vent(il,j,1)-v &
2415            (il,1))
2416        END IF ! cvflag_grav
2417      END IF ! j
2418    END DO
2419  END DO
2420
2421  ! do k=1,ntra
2422  ! do j=2,nl
2423  ! do il=1,ncum
2424  ! if (j.le.inb(il)) then
2425
2426  ! if (cvflag_grav) then
2427  ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
2428  ! :                *(traent(il,j,1,k)-tra(il,1,k))
2429  ! else
2430  ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
2431  ! :                *(traent(il,j,1,k)-tra(il,1,k))
2432  ! endif
2433
2434  ! endif
2435  ! enddo
2436  ! enddo
2437  ! enddo
2438
2439
2440  ! ***  calculate tendencies of potential temperature and mixing ratio  ***
2441  ! ***               at levels above the lowest level                   ***
2442
2443  ! ***  first find the net saturated updraft and downdraft mass fluxes  ***
2444  ! ***                      through each level                          ***
2445
2446
2447  DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
2448
2449    num1 = 0
2450    DO il = 1, ncum
2451      IF (i<=inb(il)) num1 = num1 + 1
2452    END DO
2453    IF (num1<=0) GO TO 500
2454
2455    CALL zilch(amp1, ncum)
2456    CALL zilch(ad, ncum)
2457
2458    DO k = i + 1, nl + 1
2459      DO il = 1, ncum
2460        IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN
2461          amp1(il) = amp1(il) + m(il, k)
2462        END IF
2463      END DO
2464    END DO
2465
2466    DO k = 1, i
2467      DO j = i + 1, nl + 1
2468        DO il = 1, ncum
2469          IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
2470            amp1(il) = amp1(il) + ment(il, k, j)
2471          END IF
2472        END DO
2473      END DO
2474    END DO
2475
2476    DO k = 1, i - 1
2477      DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
2478        DO il = 1, ncum
2479          IF (i<=inb(il) .AND. j<=inb(il)) THEN
2480            ad(il) = ad(il) + ment(il, j, k)
2481          END IF
2482        END DO
2483      END DO
2484    END DO
2485
2486    DO il = 1, ncum
2487      IF (i<=inb(il)) THEN
2488        dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2489        cpinv = 1.0/cpn(il, i)
2490
2491        ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
2492        IF (cvflag_grav) THEN
2493          IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2494        ELSE
2495          IF ((0.1*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2496        END IF
2497
2498        IF (cvflag_grav) THEN
2499          ft(il, i) = 0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
2500            i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
2501            i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( &
2502            il,i)+evap(il,i+1))
2503          rat = cpn(il, i-1)*cpinv
2504          ft(il, i) = ft(il, i) - 0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
2505            -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2506          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h( &
2507            il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2508        ELSE ! cvflag_grav
2509          ft(il, i) = 0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
2510            i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
2511            i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( &
2512            il,i)+evap(il,i+1))
2513          rat = cpn(il, i-1)*cpinv
2514          ft(il, i) = ft(il, i) - 0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i)-mp(il &
2515            ,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2516          ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i)+ &
2517            t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2518        END IF ! cvflag_grav
2519
2520
2521        ft(il, i) = ft(il, i) + 0.01*sigd*wt(il, i)*(cl-cpd)*water(il, i+1)*( &
2522          t(il,i+1)-t(il,i))*dpinv*cpinv
2523
2524        IF (cvflag_grav) THEN
2525          fr(il, i) = 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
2526            i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2527          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
2528            i))-ad(il)*(u(il,i)-u(il,i-1)))
2529          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
2530            i))-ad(il)*(v(il,i)-v(il,i-1)))
2531        ELSE ! cvflag_grav
2532          fr(il, i) = 0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
2533            i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2534          fu(il, i) = fu(il, i) + 0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
2535            i))-ad(il)*(u(il,i)-u(il,i-1)))
2536          fv(il, i) = fv(il, i) + 0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
2537            i))-ad(il)*(v(il,i)-v(il,i-1)))
2538        END IF ! cvflag_grav
2539
2540      END IF ! i
2541    END DO
2542
2543    ! do k=1,ntra
2544    ! do il=1,ncum
2545    ! if (i.le.inb(il)) then
2546    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2547    ! cpinv=1.0/cpn(il,i)
2548    ! if (cvflag_grav) then
2549    ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
2550    ! :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2551    ! :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2552    ! else
2553    ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2554    ! :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2555    ! :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2556    ! endif
2557    ! endif
2558    ! enddo
2559    ! enddo
2560
2561    DO k = 1, i - 1
2562      DO il = 1, ncum
2563        IF (i<=inb(il)) THEN
2564          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2565          cpinv = 1.0/cpn(il, i)
2566
2567          awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
2568          awat = amax1(awat, 0.0)
2569
2570          IF (cvflag_grav) THEN
2571            fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
2572              ,i)-awat-rr(il,i))
2573            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2574              ,i)-u(il,i))
2575            fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
2576              ,i)-v(il,i))
2577          ELSE ! cvflag_grav
2578            fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)- &
2579              awat-rr(il,i))
2580            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2581              ,i)-u(il,i))
2582            fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
2583              il,i))
2584          END IF ! cvflag_grav
2585
2586          ! (saturated updrafts resulting from mixing)        ! cld
2587          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld
2588          nqcond(il, i) = nqcond(il, i) + 1. ! cld
2589        END IF ! i
2590      END DO
2591    END DO
2592
2593    ! do j=1,ntra
2594    ! do k=1,i-1
2595    ! do il=1,ncum
2596    ! if (i.le.inb(il)) then
2597    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2598    ! cpinv=1.0/cpn(il,i)
2599    ! if (cvflag_grav) then
2600    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2601    ! :        *(traent(il,k,i,j)-tra(il,i,j))
2602    ! else
2603    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2604    ! :        *(traent(il,k,i,j)-tra(il,i,j))
2605    ! endif
2606    ! endif
2607    ! enddo
2608    ! enddo
2609    ! enddo
2610
2611    DO k = i, nl + 1
2612      DO il = 1, ncum
2613        IF (i<=inb(il) .AND. k<=inb(il)) THEN
2614          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2615          cpinv = 1.0/cpn(il, i)
2616
2617          IF (cvflag_grav) THEN
2618            fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
2619              ,i)-rr(il,i))
2620            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2621              ,i)-u(il,i))
2622            fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
2623              ,i)-v(il,i))
2624          ELSE ! cvflag_grav
2625            fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)-rr &
2626              (il,i))
2627            fu(il, i) = fu(il, i) + 0.1*dpinv*ment(il, k, i)*(uent(il,k,i)-u( &
2628              il,i))
2629            fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
2630              il,i))
2631          END IF ! cvflag_grav
2632        END IF ! i and k
2633      END DO
2634    END DO
2635
2636    ! do j=1,ntra
2637    ! do k=i,nl+1
2638    ! do il=1,ncum
2639    ! if (i.le.inb(il) .and. k.le.inb(il)) then
2640    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2641    ! cpinv=1.0/cpn(il,i)
2642    ! if (cvflag_grav) then
2643    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2644    ! :         *(traent(il,k,i,j)-tra(il,i,j))
2645    ! else
2646    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2647    ! :             *(traent(il,k,i,j)-tra(il,i,j))
2648    ! endif
2649    ! endif ! i and k
2650    ! enddo
2651    ! enddo
2652    ! enddo
2653
2654    DO il = 1, ncum
2655      IF (i<=inb(il)) THEN
2656        dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2657        cpinv = 1.0/cpn(il, i)
2658
2659        IF (cvflag_grav) THEN
2660          ! sb: on ne fait pas encore la correction permettant de mieux
2661          ! conserver l'eau:
2662          fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + &
2663            0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il, &
2664            i)-rr(il,i-1)))*dpinv
2665
2666          fu(il, i) = fu(il, i) + 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il, &
2667            i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2668          fv(il, i) = fv(il, i) + 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il, &
2669            i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2670        ELSE ! cvflag_grav
2671          fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + &
2672            0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, &
2673            i-1)))*dpinv
2674          fu(il, i) = fu(il, i) + 0.1*(mp(il,i+1)*(up(il,i+1)-u(il, &
2675            i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2676          fv(il, i) = fv(il, i) + 0.1*(mp(il,i+1)*(vp(il,i+1)-v(il, &
2677            i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2678        END IF ! cvflag_grav
2679
2680      END IF ! i
2681    END DO
2682
2683    ! sb: interface with the cloud parameterization:          ! cld
2684
2685    DO k = i + 1, nl
2686      DO il = 1, ncum
2687        IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
2688          ! (saturated downdrafts resulting from mixing)            ! cld
2689          qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
2690          nqcond(il, i) = nqcond(il, i) + 1. ! cld
2691        END IF ! cld
2692      END DO ! cld
2693    END DO ! cld
2694
2695    ! (particular case: no detraining level is found)         ! cld
2696    DO il = 1, ncum ! cld
2697      IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld
2698        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
2699        nqcond(il, i) = nqcond(il, i) + 1. ! cld
2700      END IF ! cld
2701    END DO ! cld
2702
2703    DO il = 1, ncum ! cld
2704      IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld
2705        qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
2706      END IF ! cld
2707    END DO
2708
2709    ! do j=1,ntra
2710    ! do il=1,ncum
2711    ! if (i.le.inb(il)) then
2712    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2713    ! cpinv=1.0/cpn(il,i)
2714
2715    ! if (cvflag_grav) then
2716    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
2717    ! :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2718    ! :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2719    ! else
2720    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2721    ! :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2722    ! :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2723    ! endif
2724    ! endif ! i
2725    ! enddo
2726    ! enddo
2727
2728500 END DO
2729
2730
2731  ! ***   move the detrainment at level inb down to level inb-1   ***
2732  ! ***        in such a way as to preserve the vertically        ***
2733  ! ***          integrated enthalpy and water tendencies         ***
2734
2735  DO il = 1, ncum
2736
2737    ax = 0.1*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t(il, &
2738      inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
2739      inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2740    ft(il, inb(il)) = ft(il, inb(il)) - ax
2741    ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il &
2742      ))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, &
2743      inb(il))))
2744
2745    bx = 0.1*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb( &
2746      il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2747    fr(il, inb(il)) = fr(il, inb(il)) - bx
2748    fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ &
2749      1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2750
2751    cx = 0.1*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il &
2752      )))/(ph(il,inb(il))-ph(il,inb(il)+1))
2753    fu(il, inb(il)) = fu(il, inb(il)) - cx
2754    fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ &
2755      1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2756
2757    dx = 0.1*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il &
2758      )))/(ph(il,inb(il))-ph(il,inb(il)+1))
2759    fv(il, inb(il)) = fv(il, inb(il)) - dx
2760    fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ &
2761      1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2762
2763  END DO
2764
2765  ! do j=1,ntra
2766  ! do il=1,ncum
2767  ! ex=0.1*ment(il,inb(il),inb(il))
2768  ! :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2769  ! :      /(ph(il,inb(il))-ph(il,inb(il)+1))
2770  ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2771  ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2772  ! :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2773  ! :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
2774  ! enddo
2775  ! enddo
2776
2777
2778  ! ***    homoginize tendencies below cloud base    ***
2779
2780
2781  DO il = 1, ncum
2782    asum(il) = 0.0
2783    bsum(il) = 0.0
2784    csum(il) = 0.0
2785    dsum(il) = 0.0
2786  END DO
2787
2788  DO i = 1, nl
2789    DO il = 1, ncum
2790      IF (i<=(icb(il)-1)) THEN
2791        asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1))
2792        bsum(il) = bsum(il) + fr(il, i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2793          1)))*(ph(il,i)-ph(il,i+1))
2794        csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2795          1)))*(ph(il,i)-ph(il,i+1))
2796        dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
2797      END IF
2798    END DO
2799  END DO
2800
2801  ! !!!      do 700 i=1,icb(il)-1
2802  DO i = 1, nl
2803    DO il = 1, ncum
2804      IF (i<=(icb(il)-1)) THEN
2805        ft(il, i) = asum(il)*t(il, i)/(th(il,i)*dsum(il))
2806        fr(il, i) = bsum(il)/csum(il)
2807      END IF
2808    END DO
2809  END DO
2810
2811
2812  ! ***           reset counter and return           ***
2813
2814  DO il = 1, ncum
2815    sig(il, nd) = 2.0
2816  END DO
2817
2818
2819  DO i = 1, nd
2820    DO il = 1, ncum
2821      upwd(il, i) = 0.0
2822      dnwd(il, i) = 0.0
2823    END DO
2824  END DO
2825
2826  DO i = 1, nl
2827    DO il = 1, ncum
2828      dnwd0(il, i) = -mp(il, i)
2829    END DO
2830  END DO
2831  DO i = nl + 1, nd
2832    DO il = 1, ncum
2833      dnwd0(il, i) = 0.
2834    END DO
2835  END DO
2836
2837
2838  DO i = 1, nl
2839    DO il = 1, ncum
2840      IF (i>=icb(il) .AND. i<=inb(il)) THEN
2841        upwd(il, i) = 0.0
2842        dnwd(il, i) = 0.0
2843      END IF
2844    END DO
2845  END DO
2846
2847  DO i = 1, nl
2848    DO k = 1, nl
2849      DO il = 1, ncum
2850        up1(il, k, i) = 0.0
2851        dn1(il, k, i) = 0.0
2852      END DO
2853    END DO
2854  END DO
2855
2856  DO i = 1, nl
2857    DO k = i, nl
2858      DO n = 1, i - 1
2859        DO il = 1, ncum
2860          IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
2861            up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
2862            dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
2863          END IF
2864        END DO
2865      END DO
2866    END DO
2867  END DO
2868
2869  DO i = 2, nl
2870    DO k = i, nl
2871      DO il = 1, ncum
2872        ! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
2873        ! then
2874        IF (i<=inb(il) .AND. k<=inb(il)) THEN
2875          upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
2876          dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
2877        END IF
2878      END DO
2879    END DO
2880  END DO
2881
2882
2883  ! !!!      DO il=1,ncum
2884  ! !!!      do i=icb(il),inb(il)
2885  ! !!!
2886  ! !!!      upwd(il,i)=0.0
2887  ! !!!      dnwd(il,i)=0.0
2888  ! !!!      do k=i,inb(il)
2889  ! !!!      up1=0.0
2890  ! !!!      dn1=0.0
2891  ! !!!      do n=1,i-1
2892  ! !!!      up1=up1+ment(il,n,k)
2893  ! !!!      dn1=dn1-ment(il,k,n)
2894  ! !!!      enddo
2895  ! !!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
2896  ! !!!      dnwd(il,i)=dnwd(il,i)+dn1
2897  ! !!!      enddo
2898  ! !!!      enddo
2899  ! !!!
2900  ! !!!      ENDDO
2901
2902  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2903  ! determination de la variation de flux ascendant entre
2904  ! deux niveau non dilue mike
2905  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2906
2907  DO i = 1, nl
2908    DO il = 1, ncum
2909      mike(il, i) = m(il, i)
2910    END DO
2911  END DO
2912
2913  DO i = nl + 1, nd
2914    DO il = 1, ncum
2915      mike(il, i) = 0.
2916    END DO
2917  END DO
2918
2919  DO i = 1, nd
2920    DO il = 1, ncum
2921      ma(il, i) = 0
2922    END DO
2923  END DO
2924
2925  DO i = 1, nl
2926    DO j = i, nl
2927      DO il = 1, ncum
2928        ma(il, i) = ma(il, i) + m(il, j)
2929      END DO
2930    END DO
2931  END DO
2932
2933  DO i = nl + 1, nd
2934    DO il = 1, ncum
2935      ma(il, i) = 0.
2936    END DO
2937  END DO
2938
2939  DO i = 1, nl
2940    DO il = 1, ncum
2941      IF (i<=(icb(il)-1)) THEN
2942        ma(il, i) = 0
2943      END IF
2944    END DO
2945  END DO
2946
2947  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2948  ! icb represente de niveau ou se trouve la
2949  ! base du nuage , et inb le top du nuage
2950  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2951
2952  DO i = 1, nd
2953    DO il = 1, ncum
2954      mke(il, i) = upwd(il, i) + dnwd(il, i)
2955    END DO
2956  END DO
2957
2958  DO i = 1, nd
2959    DO il = 1, ncum
2960      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il, &
2961        i))+rr(il,i)*cpv)
2962      tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp
2963      tps(il, i) = tp(il, i)
2964    END DO
2965  END DO
2966
2967
2968  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
2969  ! ***           of condensed water         ***            ! cld
2970  ! ! cld
2971
2972  DO i = 1, nd ! cld
2973    DO il = 1, ncum ! cld
2974      mac(il, i) = 0.0 ! cld
2975      wa(il, i) = 0.0 ! cld
2976      siga(il, i) = 0.0 ! cld
2977      sax(il, i) = 0.0 ! cld
2978    END DO ! cld
2979  END DO ! cld
2980
2981  DO i = minorig, nl ! cld
2982    DO k = i + 1, nl + 1 ! cld
2983      DO il = 1, ncum ! cld
2984        IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld
2985          mac(il, i) = mac(il, i) + m(il, k) ! cld
2986        END IF ! cld
2987      END DO ! cld
2988    END DO ! cld
2989  END DO ! cld
2990
2991  DO i = 1, nl ! cld
2992    DO j = 1, i ! cld
2993      DO il = 1, ncum ! cld
2994        IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
2995            .AND. j>=icb(il)) THEN ! cld
2996          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
2997            *(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
2998        END IF ! cld
2999      END DO ! cld
3000    END DO ! cld
3001  END DO ! cld
3002
3003  DO i = 1, nl ! cld
3004    DO il = 1, ncum ! cld
3005      IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3006          .AND. sax(il,i)>0.0) THEN ! cld
3007        wa(il, i) = sqrt(2.*sax(il,i)) ! cld
3008      END IF ! cld
3009    END DO ! cld
3010  END DO ! cld
3011
3012  DO i = 1, nl ! cld
3013    DO il = 1, ncum ! cld
3014      IF (wa(il,i)>0.0) &          ! cld
3015        siga(il, i) = mac(il, i)/wa(il, i) & ! cld
3016        *rrd*tvp(il, i)/p(il, i)/100./delta ! cld
3017      siga(il, i) = min(siga(il,i), 1.0) ! cld
3018      ! IM cf. FH
3019      IF (iflag_clw==0) THEN
3020        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) & ! cld
3021          +(1.-siga(il,i))*qcond(il, i) ! cld
3022      ELSE IF (iflag_clw==1) THEN
3023        qcondc(il, i) = qcond(il, i) ! cld
3024      END IF
3025
3026    END DO ! cld
3027  END DO ! cld
3028
3029  RETURN
3030END SUBROUTINE cv30_yield
3031
3032! !RomP >>>
3033SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
3034    d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
3035  IMPLICIT NONE
3036
3037  include "cv30param.h"
3038
3039  ! inputs:
3040  INTEGER ncum, nd, na, nloc, len
3041  REAL ment(nloc, na, na), sij(nloc, na, na)
3042  REAL clw(nloc, nd), elij(nloc, na, na)
3043  REAL ep(nloc, na)
3044  INTEGER icb(nloc), inb(nloc)
3045  REAL vprecip(nloc, nd+1)
3046  ! ouputs:
3047  REAL da(nloc, na), phi(nloc, na, na)
3048  REAL phi2(nloc, na, na)
3049  REAL d1a(nloc, na), dam(nloc, na)
3050  REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
3051  ! variables pour tracer dans precip de l'AA et des mel
3052  ! local variables:
3053  INTEGER i, j, k
3054  REAL epm(nloc, na, na)
3055
3056  ! variables d'Emanuel : du second indice au troisieme
3057  ! --->    tab(i,k,j) -> de l origine k a l arrivee j
3058  ! ment, sij, elij
3059  ! variables personnelles : du troisieme au second indice
3060  ! --->    tab(i,j,k) -> de k a j
3061  ! phi, phi2
3062
3063  ! initialisations
3064  DO j = 1, na
3065    DO i = 1, ncum
3066      da(i, j) = 0.
3067      d1a(i, j) = 0.
3068      dam(i, j) = 0.
3069      eplamm(i, j) = 0.
3070    END DO
3071  END DO
3072  DO k = 1, na
3073    DO j = 1, na
3074      DO i = 1, ncum
3075        epm(i, j, k) = 0.
3076        epmlmmm(i, j, k) = 0.
3077        phi(i, j, k) = 0.
3078        phi2(i, j, k) = 0.
3079      END DO
3080    END DO
3081  END DO
3082
3083  ! fraction deau condensee dans les melanges convertie en precip : epm
3084  ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3085  DO j = 1, na
3086    DO k = 1, j - 1
3087      DO i = 1, ncum
3088        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3089          ! !jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3090          epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3091          ! !
3092          epm(i, j, k) = max(epm(i,j,k), 0.0)
3093        END IF
3094      END DO
3095    END DO
3096  END DO
3097
3098  DO j = 1, na
3099    DO k = 1, na
3100      DO i = 1, ncum
3101        IF (k>=icb(i) .AND. k<=inb(i)) THEN
3102          eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- &
3103            sij(i,j,k))
3104        END IF
3105      END DO
3106    END DO
3107  END DO
3108
3109  DO j = 1, na
3110    DO k = 1, j - 1
3111      DO i = 1, ncum
3112        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3113          epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
3114        END IF
3115      END DO
3116    END DO
3117  END DO
3118
3119  ! matrices pour calculer la tendance des concentrations dans cvltr.F90
3120  DO j = 1, na
3121    DO k = 1, na
3122      DO i = 1, ncum
3123        da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
3124        phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
3125        d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sij(i,k,j))
3126      END DO
3127    END DO
3128  END DO
3129
3130  DO j = 1, na
3131    DO k = 1, j - 1
3132      DO i = 1, ncum
3133        dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, j, k)*(1.-ep(i,k))*(1.- &
3134          sij(i,k,j))
3135        phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
3136      END DO
3137    END DO
3138  END DO
3139
3140  RETURN
3141END SUBROUTINE cv30_tracer
3142! RomP <<<
3143
3144SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
3145    vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
3146    dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
3147    epmlmmm, eplamm, wdtraina, wdtrainm,epmax_diag, iflag1, precip1, vprecip1, evap1, &
3148    ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
3149    dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
3150    elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1,epmax_diag1) ! epmax_cape
3151  IMPLICIT NONE
3152
3153  include "cv30param.h"
3154
3155  ! inputs:
3156  INTEGER len, ncum, nd, ntra, nloc
3157  INTEGER idcum(nloc)
3158  INTEGER iflag(nloc)
3159  INTEGER inb(nloc)
3160  REAL precip(nloc)
3161  REAL vprecip(nloc, nd+1), evap(nloc, nd)
3162  REAL ep(nloc, nd)
3163  REAL sig(nloc, nd), w0(nloc, nd)
3164  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3165  REAL ftra(nloc, nd, ntra)
3166  REAL ma(nloc, nd)
3167  REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3168  REAL qcondc(nloc, nd)
3169  REAL wd(nloc), cape(nloc)
3170  REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
3171  REAL epmax_diag(nloc) ! epmax_cape
3172  ! RomP >>>
3173  REAL phi2(nloc, nd, nd)
3174  REAL d1a(nloc, nd), dam(nloc, nd)
3175  REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
3176  REAL sij(nloc, nd, nd)
3177  REAL elij(nloc, nd, nd), clw(nloc, nd)
3178  REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
3179  ! RomP <<<
3180
3181  ! outputs:
3182  INTEGER iflag1(len)
3183  INTEGER inb1(len)
3184  REAL precip1(len)
3185  REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP
3186  REAL ep1(len, nd) !<<< RomP
3187  REAL sig1(len, nd), w01(len, nd)
3188  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
3189  REAL ftra1(len, nd, ntra)
3190  REAL ma1(len, nd)
3191  REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
3192  REAL qcondc1(nloc, nd)
3193  REAL wd1(nloc), cape1(nloc)
3194  REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
3195  REAL epmax_diag1(len) ! epmax_cape
3196  ! RomP >>>
3197  REAL phi21(len, nd, nd)
3198  REAL d1a1(len, nd), dam1(len, nd)
3199  REAL wdtraina1(len, nd), wdtrainm1(len, nd)
3200  REAL sij1(len, nd, nd)
3201  REAL elij1(len, nd, nd), clw1(len, nd)
3202  REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
3203  ! RomP <<<
3204
3205  ! local variables:
3206  INTEGER i, k, j
3207
3208  DO i = 1, ncum
3209    precip1(idcum(i)) = precip(i)
3210    iflag1(idcum(i)) = iflag(i)
3211    wd1(idcum(i)) = wd(i)
3212    inb1(idcum(i)) = inb(i)
3213    cape1(idcum(i)) = cape(i)
3214    epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
3215  END DO
3216
3217  DO k = 1, nl
3218    DO i = 1, ncum
3219      vprecip1(idcum(i), k) = vprecip(i, k)
3220      evap1(idcum(i), k) = evap(i, k) !<<< RomP
3221      sig1(idcum(i), k) = sig(i, k)
3222      w01(idcum(i), k) = w0(i, k)
3223      ft1(idcum(i), k) = ft(i, k)
3224      fq1(idcum(i), k) = fq(i, k)
3225      fu1(idcum(i), k) = fu(i, k)
3226      fv1(idcum(i), k) = fv(i, k)
3227      ma1(idcum(i), k) = ma(i, k)
3228      upwd1(idcum(i), k) = upwd(i, k)
3229      dnwd1(idcum(i), k) = dnwd(i, k)
3230      dnwd01(idcum(i), k) = dnwd0(i, k)
3231      qcondc1(idcum(i), k) = qcondc(i, k)
3232      da1(idcum(i), k) = da(i, k)
3233      mp1(idcum(i), k) = mp(i, k)
3234      ! RomP >>>
3235      ep1(idcum(i), k) = ep(i, k)
3236      d1a1(idcum(i), k) = d1a(i, k)
3237      dam1(idcum(i), k) = dam(i, k)
3238      clw1(idcum(i), k) = clw(i, k)
3239      eplamm1(idcum(i), k) = eplamm(i, k)
3240      wdtraina1(idcum(i), k) = wdtraina(i, k)
3241      wdtrainm1(idcum(i), k) = wdtrainm(i, k)
3242      ! RomP <<<
3243    END DO
3244  END DO
3245
3246  DO i = 1, ncum
3247    sig1(idcum(i), nd) = sig(i, nd)
3248  END DO
3249
3250
3251  ! do 2100 j=1,ntra
3252  ! do 2110 k=1,nd ! oct3
3253  ! do 2120 i=1,ncum
3254  ! ftra1(idcum(i),k,j)=ftra(i,k,j)
3255  ! 2120     continue
3256  ! 2110    continue
3257  ! 2100   continue
3258  DO j = 1, nd
3259    DO k = 1, nd
3260      DO i = 1, ncum
3261        sij1(idcum(i), k, j) = sij(i, k, j)
3262        phi1(idcum(i), k, j) = phi(i, k, j)
3263        phi21(idcum(i), k, j) = phi2(i, k, j)
3264        elij1(idcum(i), k, j) = elij(i, k, j)
3265        epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
3266      END DO
3267    END DO
3268  END DO
3269
3270  RETURN
3271END SUBROUTINE cv30_uncompress
3272
3273        subroutine cv30_epmax_fn_cape(nloc,ncum,nd &
3274                ,cape,ep,hp,icb,inb,clw,nk,t,h,lv &
3275                ,epmax_diag)
3276        implicit none
3277
3278        ! On fait varier epmax en fn de la cape
3279        ! Il faut donc recalculer ep, et hp qui a déjà été calculé et
3280        ! qui en dépend
3281        ! Toutes les autres variables fn de ep sont calculées plus bas.
3282
3283#include "cvthermo.h"
3284#include "cv30param.h"
3285#include "conema3.h"
3286
3287! inputs:
3288      integer ncum, nd, nloc
3289      integer icb(nloc), inb(nloc)
3290      real cape(nloc)
3291      real clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd)
3292      integer nk(nloc)
3293! inouts:
3294      real ep(nloc,nd)
3295      real hp(nloc,nd)
3296! outputs ou local
3297      real epmax_diag(nloc)
3298! locals
3299      integer i,k   
3300      real hp_bak(nloc,nd)
3301
3302        ! on recalcule ep et hp
3303       
3304        if (coef_epmax_cape.gt.1e-12) then
3305        do i=1,ncum
3306           epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
3307           do k=1,nl
3308                ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
3309                ep(i,k)=amax1(ep(i,k),0.0)
3310                ep(i,k)=amin1(ep(i,k),epmax_diag(i))
3311           enddo
3312        enddo
3313
3314! On recalcule hp:
3315      do k=1,nl
3316        do i=1,ncum
3317          hp_bak(i,k)=hp(i,k)
3318        enddo
3319      enddo
3320      do k=1,nlp
3321        do i=1,ncum
3322          hp(i,k)=h(i,k)
3323        enddo
3324      enddo
3325      do k=minorig+1,nl
3326       do i=1,ncum
3327        if((k.ge.icb(i)).and.(k.le.inb(i)))then
3328          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
3329        endif
3330       enddo
3331      enddo !do k=minorig+1,n
3332!     write(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
3333      do i=1,ncum 
3334       do k=1,nl
3335        if (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) then
3336           write(*,*) 'i,k=',i,k
3337           write(*,*) 'coef_epmax_cape=',coef_epmax_cape
3338           write(*,*) 'epmax_diag(i)=',epmax_diag(i)
3339           write(*,*) 'ep(i,k)=',ep(i,k)
3340           write(*,*) 'hp(i,k)=',hp(i,k)
3341           write(*,*) 'hp_bak(i,k)=',hp_bak(i,k)
3342           write(*,*) 'h(i,k)=',h(i,k)
3343           write(*,*) 'nk(i)=',nk(i)
3344           write(*,*) 'h(i,nk(i))=',h(i,nk(i))
3345           write(*,*) 'lv(i,k)=',lv(i,k)
3346           write(*,*) 't(i,k)=',t(i,k)
3347           write(*,*) 'clw(i,k)=',clw(i,k)
3348           write(*,*) 'cpd,cpv=',cpd,cpv
3349           stop
3350        endif
3351       enddo !do k=1,nl
3352      enddo !do i=1,ncum 
3353      endif !if (coef_epmax_cape.gt.1e-12) then
3354
3355      return
3356      end subroutine cv30_epmax_fn_cape
3357
3358
Note: See TracBrowser for help on using the repository browser.