source: LMDZ5/branches/testing/libf/phylmd/cv30_routines.F90 @ 2435

Last change on this file since 2435 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch
  • 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: 95.8 KB
Line 
1
2! $Id: cv30_routines.F90 2408 2015-12-14 10:43:09Z fairhead $
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  IMPLICIT NONE
842
843  ! ---------------------------------------------------------------------
844  ! Purpose:
845  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
846  ! &
847  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
848  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
849  ! &
850  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
851
852  ! Main differences convect3/convect4:
853  ! - icbs (input) is the first level above LCL (may differ from icb)
854  ! - many minor differences in the iterations
855  ! - condensed water not removed from tvp in convect3
856  ! - vertical profile of buoyancy computed here (use of buoybase)
857  ! - the determination of inb is different
858  ! - no inb1, only inb in output
859  ! ---------------------------------------------------------------------
860
861  include "cvthermo.h"
862  include "cv30param.h"
863  include "conema3.h"
864
865  ! inputs:
866  INTEGER ncum, nd, nloc
867  INTEGER icb(nloc), icbs(nloc), nk(nloc)
868  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
869  REAL p(nloc, nd)
870  REAL tnk(nloc), qnk(nloc), gznk(nloc)
871  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
872  REAL pbase(nloc), buoybase(nloc), plcl(nloc)
873
874  ! outputs:
875  INTEGER inb(nloc)
876  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
877  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
878  REAL buoy(nloc, nd)
879
880  ! local variables:
881  INTEGER i, k
882  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
883  REAL by, defrac, pden
884  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
885  LOGICAL lcape(nloc)
886
887  ! =====================================================================
888  ! --- SOME INITIALIZATIONS
889  ! =====================================================================
890
891  DO k = 1, nl
892    DO i = 1, ncum
893      ep(i, k) = 0.0
894      sigp(i, k) = spfac
895    END DO
896  END DO
897
898  ! =====================================================================
899  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
900  ! =====================================================================
901
902  ! ---       The procedure is to solve the equation.
903  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
904
905  ! ***  Calculate certain parcel quantities, including static energy   ***
906
907
908  DO i = 1, ncum
909    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) & ! debug     &
910                                                  ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
911      +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
912  END DO
913
914
915  ! ***  Find lifted parcel quantities above cloud base    ***
916
917
918  DO k = minorig + 1, nl
919    DO i = 1, ncum
920      ! ori         if(k.ge.(icb(i)+1))then
921      IF (k>=(icbs(i)+1)) THEN ! convect3
922        tg = t(i, k)
923        qg = qs(i, k)
924        ! debug       alv=lv0-clmcpv*(t(i,k)-t0)
925        alv = lv0 - clmcpv*(t(i,k)-273.15)
926
927        ! First iteration.
928
929        ! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
930        s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
931          +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
932        s = 1./s
933        ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
934        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
935        tg = tg + s*(ah0(i)-ahg)
936        ! ori          tg=max(tg,35.0)
937        ! debug        tc=tg-t0
938        tc = tg - 273.15
939        denom = 243.5 + tc
940        denom = max(denom, 1.0) ! convect3
941        ! ori          if(tc.ge.0.0)then
942        es = 6.112*exp(17.67*tc/denom)
943        ! ori          else
944        ! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
945        ! ori          endif
946        qg = eps*es/(p(i,k)-es*(1.-eps))
947
948        ! Second iteration.
949
950        ! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
951        ! ori          s=1./s
952        ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
953        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
954        tg = tg + s*(ah0(i)-ahg)
955        ! ori          tg=max(tg,35.0)
956        ! debug        tc=tg-t0
957        tc = tg - 273.15
958        denom = 243.5 + tc
959        denom = max(denom, 1.0) ! convect3
960        ! ori          if(tc.ge.0.0)then
961        es = 6.112*exp(17.67*tc/denom)
962        ! ori          else
963        ! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
964        ! ori          endif
965        qg = eps*es/(p(i,k)-es*(1.-eps))
966
967        ! debug        alv=lv0-clmcpv*(t(i,k)-t0)
968        alv = lv0 - clmcpv*(t(i,k)-273.15)
969        ! print*,'cpd dans convect2 ',cpd
970        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
971        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
972
973        ! ori c approximation here:
974        ! ori
975        ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
976
977        ! convect3: no approximation:
978        tp(i, k) = (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
979
980        clw(i, k) = qnk(i) - qg
981        clw(i, k) = max(0.0, clw(i,k))
982        rg = qg/(1.-qnk(i))
983        ! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
984        ! convect3: (qg utilise au lieu du vrai mixing ratio rg):
985        tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
986      END IF
987    END DO
988  END DO
989
990  ! =====================================================================
991  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
992  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
993  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
994  ! =====================================================================
995
996  ! ori      do 320 k=minorig+1,nl
997  DO k = 1, nl ! convect3
998    DO i = 1, ncum
999      pden = ptcrit - pbcrit
1000      ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1001      ep(i, k) = amax1(ep(i,k), 0.0)
1002      ep(i, k) = amin1(ep(i,k), epmax)
1003      sigp(i, k) = spfac
1004      ! ori          if(k.ge.(nk(i)+1))then
1005      ! ori            tca=tp(i,k)-t0
1006      ! ori            if(tca.ge.0.0)then
1007      ! ori              elacrit=elcrit
1008      ! ori            else
1009      ! ori              elacrit=elcrit*(1.0-tca/tlcrit)
1010      ! ori            endif
1011      ! ori            elacrit=max(elacrit,0.0)
1012      ! ori            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
1013      ! ori            ep(i,k)=max(ep(i,k),0.0 )
1014      ! ori            ep(i,k)=min(ep(i,k),1.0 )
1015      ! ori            sigp(i,k)=sigs
1016      ! ori          endif
1017    END DO
1018  END DO
1019
1020  ! =====================================================================
1021  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1022  ! --- VIRTUAL TEMPERATURE
1023  ! =====================================================================
1024
1025  ! dans convect3, tvp est calcule en une seule fois, et sans retirer
1026  ! l'eau condensee (~> reversible CAPE)
1027
1028  ! ori      do 340 k=minorig+1,nl
1029  ! ori        do 330 i=1,ncum
1030  ! ori        if(k.ge.(icb(i)+1))then
1031  ! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
1032  ! oric         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1033  ! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1034  ! ori        endif
1035  ! ori 330    continue
1036  ! ori 340  continue
1037
1038  ! ori      do 350 i=1,ncum
1039  ! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1040  ! ori 350  continue
1041
1042  DO i = 1, ncum ! convect3
1043    tp(i, nlp) = tp(i, nl) ! convect3
1044  END DO ! convect3
1045
1046  ! =====================================================================
1047  ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1048  ! =====================================================================
1049
1050  ! -- this is for convect3 only:
1051
1052  ! first estimate of buoyancy:
1053
1054  DO i = 1, ncum
1055    DO k = 1, nl
1056      buoy(i, k) = tvp(i, k) - tv(i, k)
1057    END DO
1058  END DO
1059
1060  ! set buoyancy=buoybase for all levels below base
1061  ! for safety, set buoy(icb)=buoybase
1062
1063  DO i = 1, ncum
1064    DO k = 1, nl
1065      IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1066        buoy(i, k) = buoybase(i)
1067      END IF
1068    END DO
1069    ! IM cf. CRio/JYG 270807   buoy(icb(i),k)=buoybase(i)
1070    buoy(i, icb(i)) = buoybase(i)
1071  END DO
1072
1073  ! -- end convect3
1074
1075  ! =====================================================================
1076  ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1077  ! --- LEVEL OF NEUTRAL BUOYANCY
1078  ! =====================================================================
1079
1080  ! -- this is for convect3 only:
1081
1082  DO i = 1, ncum
1083    inb(i) = nl - 1
1084  END DO
1085
1086  DO i = 1, ncum
1087    DO k = 1, nl - 1
1088      IF ((k>=icb(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1089        inb(i) = min(inb(i), k)
1090      END IF
1091    END DO
1092  END DO
1093
1094  ! -- end convect3
1095
1096  ! ori      do 510 i=1,ncum
1097  ! ori        cape(i)=0.0
1098  ! ori        capem(i)=0.0
1099  ! ori        inb(i)=icb(i)+1
1100  ! ori        inb1(i)=inb(i)
1101  ! ori 510  continue
1102
1103  ! Originial Code
1104
1105  ! do 530 k=minorig+1,nl-1
1106  ! do 520 i=1,ncum
1107  ! if(k.ge.(icb(i)+1))then
1108  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1109  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1110  ! cape(i)=cape(i)+by
1111  ! if(by.ge.0.0)inb1(i)=k+1
1112  ! if(cape(i).gt.0.0)then
1113  ! inb(i)=k+1
1114  ! capem(i)=cape(i)
1115  ! endif
1116  ! endif
1117  ! 520    continue
1118  ! 530  continue
1119  ! do 540 i=1,ncum
1120  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1121  ! cape(i)=capem(i)+byp
1122  ! defrac=capem(i)-cape(i)
1123  ! defrac=max(defrac,0.001)
1124  ! frac(i)=-cape(i)/defrac
1125  ! frac(i)=min(frac(i),1.0)
1126  ! frac(i)=max(frac(i),0.0)
1127  ! 540   continue
1128
1129  ! K Emanuel fix
1130
1131  ! call zilch(byp,ncum)
1132  ! do 530 k=minorig+1,nl-1
1133  ! do 520 i=1,ncum
1134  ! if(k.ge.(icb(i)+1))then
1135  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1136  ! cape(i)=cape(i)+by
1137  ! if(by.ge.0.0)inb1(i)=k+1
1138  ! if(cape(i).gt.0.0)then
1139  ! inb(i)=k+1
1140  ! capem(i)=cape(i)
1141  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1142  ! endif
1143  ! endif
1144  ! 520    continue
1145  ! 530  continue
1146  ! do 540 i=1,ncum
1147  ! inb(i)=max(inb(i),inb1(i))
1148  ! cape(i)=capem(i)+byp(i)
1149  ! defrac=capem(i)-cape(i)
1150  ! defrac=max(defrac,0.001)
1151  ! frac(i)=-cape(i)/defrac
1152  ! frac(i)=min(frac(i),1.0)
1153  ! frac(i)=max(frac(i),0.0)
1154  ! 540   continue
1155
1156  ! J Teixeira fix
1157
1158  ! ori      call zilch(byp,ncum)
1159  ! ori      do 515 i=1,ncum
1160  ! ori        lcape(i)=.true.
1161  ! ori 515  continue
1162  ! ori      do 530 k=minorig+1,nl-1
1163  ! ori        do 520 i=1,ncum
1164  ! ori          if(cape(i).lt.0.0)lcape(i)=.false.
1165  ! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
1166  ! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1167  ! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1168  ! ori            cape(i)=cape(i)+by
1169  ! ori            if(by.ge.0.0)inb1(i)=k+1
1170  ! ori            if(cape(i).gt.0.0)then
1171  ! ori              inb(i)=k+1
1172  ! ori              capem(i)=cape(i)
1173  ! ori            endif
1174  ! ori          endif
1175  ! ori 520    continue
1176  ! ori 530  continue
1177  ! ori      do 540 i=1,ncum
1178  ! ori          cape(i)=capem(i)+byp(i)
1179  ! ori          defrac=capem(i)-cape(i)
1180  ! ori          defrac=max(defrac,0.001)
1181  ! ori          frac(i)=-cape(i)/defrac
1182  ! ori          frac(i)=min(frac(i),1.0)
1183  ! ori          frac(i)=max(frac(i),0.0)
1184  ! ori 540  continue
1185
1186  ! =====================================================================
1187  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1188  ! =====================================================================
1189
1190  ! ym      do i=1,ncum*nlp
1191  ! ym       hp(i,1)=h(i,1)
1192  ! ym      enddo
1193
1194  DO k = 1, nlp
1195    DO i = 1, ncum
1196      hp(i, k) = h(i, k)
1197    END DO
1198  END DO
1199
1200  DO k = minorig + 1, nl
1201    DO i = 1, ncum
1202      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1203        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
1204          )
1205      END IF
1206    END DO
1207  END DO
1208
1209  RETURN
1210END SUBROUTINE cv30_undilute2
1211
1212SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, &
1213    sig, w0, cape, m)
1214  IMPLICIT NONE
1215
1216  ! ===================================================================
1217  ! ---  CLOSURE OF CONVECT3
1218
1219  ! vectorization: S. Bony
1220  ! ===================================================================
1221
1222  include "cvthermo.h"
1223  include "cv30param.h"
1224
1225  ! input:
1226  INTEGER ncum, nd, nloc
1227  INTEGER icb(nloc), inb(nloc)
1228  REAL pbase(nloc)
1229  REAL p(nloc, nd), ph(nloc, nd+1)
1230  REAL tv(nloc, nd), buoy(nloc, nd)
1231
1232  ! input/output:
1233  REAL sig(nloc, nd), w0(nloc, nd)
1234
1235  ! output:
1236  REAL cape(nloc)
1237  REAL m(nloc, nd)
1238
1239  ! local variables:
1240  INTEGER i, j, k, icbmax
1241  REAL deltap, fac, w, amu
1242  REAL dtmin(nloc, nd), sigold(nloc, nd)
1243
1244
1245  ! -------------------------------------------------------
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
1351  DO k = 1, nl
1352    DO i = 1, ncum
1353
1354      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1355
1356        deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1357        cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1)
1358        cape(i) = amax1(0.0, cape(i))
1359        sigold(i, k) = sig(i, k)
1360
1361        ! dtmin(i,k)=100.0
1362        ! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1363        ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1364        ! 97     continue
1365
1366        sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
1367        sig(i, k) = amax1(sig(i,k), 0.0)
1368        sig(i, k) = amin1(sig(i,k), 0.01)
1369        fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1370        w = (1.-beta)*fac*sqrt(cape(i)) + beta*w0(i, k)
1371        amu = 0.5*(sig(i,k)+sigold(i,k))*w
1372        m(i, k) = amu*0.007*p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1373        w0(i, k) = w
1374      END IF
1375
1376    END DO
1377  END DO
1378
1379  DO i = 1, ncum
1380    w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1381    m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/ &
1382      (ph(i,icb(i)+1)-ph(i,icb(i)+2))
1383    sig(i, icb(i)) = sig(i, icb(i)+1)
1384    sig(i, icb(i)-1) = sig(i, icb(i))
1385  END DO
1386
1387
1388  ! !      cape=0.0
1389  ! !      do 98 i=icb+1,inb
1390  ! !         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1391  ! !         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
1392  ! !         dcape=rrd*buoy(i-1)*deltap/p(i-1)
1393  ! !         dlnp=deltap/p(i-1)
1394  ! !         cape=amax1(0.0,cape)
1395  ! !         sigold=sig(i)
1396
1397  ! !         dtmin=100.0
1398  ! !         do 97 j=icb,i-1
1399  ! !            dtmin=amin1(dtmin,buoy(j))
1400  ! !   97    continue
1401
1402  ! !         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
1403  ! !         sig(i)=amax1(sig(i),0.0)
1404  ! !         sig(i)=amin1(sig(i),0.01)
1405  ! !         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1406  ! !         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
1407  ! !         amu=0.5*(sig(i)+sigold)*w
1408  ! !         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
1409  ! !         w0(i)=w
1410  ! !   98 continue
1411  ! !      w0(icb)=0.5*w0(icb+1)
1412  ! !      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1413  ! !      sig(icb)=sig(icb+1)
1414  ! !      sig(icb-1)=sig(icb)
1415
1416  RETURN
1417END SUBROUTINE cv30_closure
1418
1419SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
1420    u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, &
1421    vent, sij, elij, ments, qents, traent)
1422  IMPLICIT NONE
1423
1424  ! ---------------------------------------------------------------------
1425  ! a faire:
1426  ! - changer rr(il,1) -> qnk(il)
1427  ! - vectorisation de la partie normalisation des flux (do 789...)
1428  ! ---------------------------------------------------------------------
1429
1430  include "cvthermo.h"
1431  include "cv30param.h"
1432
1433  ! inputs:
1434  INTEGER ncum, nd, na, ntra, nloc
1435  INTEGER icb(nloc), inb(nloc), nk(nloc)
1436  REAL sig(nloc, nd)
1437  REAL qnk(nloc)
1438  REAL ph(nloc, nd+1)
1439  REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1440  REAL u(nloc, nd), v(nloc, nd)
1441  REAL tra(nloc, nd, ntra) ! input of convect3
1442  REAL lv(nloc, na), h(nloc, na), hp(nloc, na)
1443  REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
1444  REAL m(nloc, na) ! input of convect3
1445
1446  ! outputs:
1447  REAL ment(nloc, na, na), qent(nloc, na, na)
1448  REAL uent(nloc, na, na), vent(nloc, na, na)
1449  REAL sij(nloc, na, na), elij(nloc, na, na)
1450  REAL traent(nloc, nd, nd, ntra)
1451  REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
1452  REAL sigij(nloc, nd, nd)
1453
1454  ! local variables:
1455  INTEGER i, j, k, il, im, jm
1456  INTEGER num1, num2
1457  INTEGER nent(nloc, na)
1458  REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1459  REAL alt, smid, sjmin, sjmax, delp, delm
1460  REAL asij(nloc), smax(nloc), scrit(nloc)
1461  REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1462  REAL wgh
1463  REAL zm(nloc, na)
1464  LOGICAL lwork(nloc)
1465
1466  ! =====================================================================
1467  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1468  ! =====================================================================
1469
1470  ! ori        do 360 i=1,ncum*nlp
1471  DO j = 1, nl
1472    DO i = 1, ncum
1473      nent(i, j) = 0
1474      ! in convect3, m is computed in cv3_closure
1475      ! ori          m(i,1)=0.0
1476    END DO
1477  END DO
1478
1479  ! ori      do 400 k=1,nlp
1480  ! ori       do 390 j=1,nlp
1481  DO j = 1, nl
1482    DO k = 1, nl
1483      DO i = 1, ncum
1484        qent(i, k, j) = rr(i, j)
1485        uent(i, k, j) = u(i, j)
1486        vent(i, k, j) = v(i, j)
1487        elij(i, k, j) = 0.0
1488        ! ym            ment(i,k,j)=0.0
1489        ! ym            sij(i,k,j)=0.0
1490      END DO
1491    END DO
1492  END DO
1493
1494  ! ym
1495  ment(1:ncum, 1:nd, 1:nd) = 0.0
1496  sij(1:ncum, 1:nd, 1:nd) = 0.0
1497
1498  ! do k=1,ntra
1499  ! do j=1,nd  ! instead nlp
1500  ! do i=1,nd ! instead nlp
1501  ! do il=1,ncum
1502  ! traent(il,i,j,k)=tra(il,j,k)
1503  ! enddo
1504  ! enddo
1505  ! enddo
1506  ! enddo
1507  zm(:, :) = 0.
1508
1509  ! =====================================================================
1510  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1511  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1512  ! --- FRACTION (sij)
1513  ! =====================================================================
1514
1515  DO i = minorig + 1, nl
1516
1517    DO j = minorig, nl
1518      DO il = 1, ncum
1519        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
1520            1)) .AND. (j<=inb(il))) THEN
1521
1522          rti = rr(il, 1) - ep(il, i)*clw(il, i)
1523          bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
1524          anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
1525          denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
1526          dei = denom
1527          IF (abs(dei)<0.01) dei = 0.01
1528          sij(il, i, j) = anum/dei
1529          sij(il, i, i) = 1.0
1530          altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
1531          altem = altem/bf2
1532          cwat = clw(il, j)*(1.-ep(il,j))
1533          stemp = sij(il, i, j)
1534          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
1535            anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
1536            denom = denom + lv(il, j)*(rr(il,i)-rti)
1537            IF (abs(denom)<0.01) denom = 0.01
1538            sij(il, i, j) = anum/denom
1539            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - &
1540              rs(il, j)
1541            altem = altem - (bf2-1.)*cwat
1542          END IF
1543          IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
1544            qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti
1545            uent(il, i, j) = sij(il, i, j)*u(il, i) + &
1546              (1.-sij(il,i,j))*u(il, nk(il))
1547            vent(il, i, j) = sij(il, i, j)*v(il, i) + &
1548              (1.-sij(il,i,j))*v(il, nk(il))
1549            ! !!!      do k=1,ntra
1550            ! !!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1551            ! !!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
1552            ! !!!      end do
1553            elij(il, i, j) = altem
1554            elij(il, i, j) = amax1(0.0, elij(il,i,j))
1555            ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
1556            nent(il, i) = nent(il, i) + 1
1557          END IF
1558          sij(il, i, j) = amax1(0.0, sij(il,i,j))
1559          sij(il, i, j) = amin1(1.0, sij(il,i,j))
1560        END IF ! new
1561      END DO
1562    END DO
1563
1564    ! do k=1,ntra
1565    ! do j=minorig,nl
1566    ! do il=1,ncum
1567    ! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1568    ! :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1569    ! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1570    ! :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
1571    ! endif
1572    ! enddo
1573    ! enddo
1574    ! enddo
1575
1576
1577    ! ***   if no air can entrain at level i assume that updraft detrains
1578    ! ***
1579    ! ***   at that level and calculate detrained air flux and properties
1580    ! ***
1581
1582
1583    ! @      do 170 i=icb(il),inb(il)
1584
1585    DO il = 1, ncum
1586      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
1587        ! @      if(nent(il,i).eq.0)then
1588        ment(il, i, i) = m(il, i)
1589        qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i)
1590        uent(il, i, i) = u(il, nk(il))
1591        vent(il, i, i) = v(il, nk(il))
1592        elij(il, i, i) = clw(il, i)
1593        ! MAF      sij(il,i,i)=1.0
1594        sij(il, i, i) = 0.0
1595      END IF
1596    END DO
1597  END DO
1598
1599  ! do j=1,ntra
1600  ! do i=minorig+1,nl
1601  ! do il=1,ncum
1602  ! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1603  ! traent(il,i,i,j)=tra(il,nk(il),j)
1604  ! endif
1605  ! enddo
1606  ! enddo
1607  ! enddo
1608
1609  DO j = minorig, nl
1610    DO i = minorig, nl
1611      DO il = 1, ncum
1612        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
1613            inb(il))) THEN
1614          sigij(il, i, j) = sij(il, i, j)
1615        END IF
1616      END DO
1617    END DO
1618  END DO
1619  ! @      enddo
1620
1621  ! @170   continue
1622
1623  ! =====================================================================
1624  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
1625  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
1626  ! =====================================================================
1627
1628  ! ym      call zilch(asum,ncum*nd)
1629  ! ym      call zilch(bsum,ncum*nd)
1630  ! ym      call zilch(csum,ncum*nd)
1631  CALL zilch(asum, nloc*nd)
1632  CALL zilch(csum, nloc*nd)
1633  CALL zilch(csum, nloc*nd)
1634
1635  DO il = 1, ncum
1636    lwork(il) = .FALSE.
1637  END DO
1638
1639  DO i = minorig + 1, nl
1640
1641    num1 = 0
1642    DO il = 1, ncum
1643      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
1644    END DO
1645    IF (num1<=0) GO TO 789
1646
1647
1648    DO il = 1, ncum
1649      IF (i>=icb(il) .AND. i<=inb(il)) THEN
1650        lwork(il) = (nent(il,i)/=0)
1651        qp = rr(il, 1) - ep(il, i)*clw(il, i)
1652        anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
1653          (cpv-cpd)*t(il, i)*(qp-rr(il,i))
1654        denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
1655          (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
1656        IF (abs(denom)<0.01) denom = 0.01
1657        scrit(il) = anum/denom
1658        alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
1659        IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
1660        smax(il) = 0.0
1661        asij(il) = 0.0
1662      END IF
1663    END DO
1664
1665    DO j = nl, minorig, -1
1666
1667      num2 = 0
1668      DO il = 1, ncum
1669        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1670          il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
1671      END DO
1672      IF (num2<=0) GO TO 175
1673
1674      DO il = 1, ncum
1675        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
1676            il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
1677
1678          IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
1679            wgh = 1.0
1680            IF (j>i) THEN
1681              sjmax = amax1(sij(il,i,j+1), smax(il))
1682              sjmax = amin1(sjmax, scrit(il))
1683              smax(il) = amax1(sij(il,i,j), smax(il))
1684              sjmin = amax1(sij(il,i,j-1), smax(il))
1685              sjmin = amin1(sjmin, scrit(il))
1686              IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
1687              smid = amin1(sij(il,i,j), scrit(il))
1688            ELSE
1689              sjmax = amax1(sij(il,i,j+1), scrit(il))
1690              smid = amax1(sij(il,i,j), scrit(il))
1691              sjmin = 0.0
1692              IF (j>1) sjmin = sij(il, i, j-1)
1693              sjmin = amax1(sjmin, scrit(il))
1694            END IF
1695            delp = abs(sjmax-smid)
1696            delm = abs(sjmin-smid)
1697            asij(il) = asij(il) + wgh*(delp+delm)
1698            ment(il, i, j) = ment(il, i, j)*(delp+delm)*wgh
1699          END IF
1700        END IF
1701      END DO
1702
1703175 END DO
1704
1705    DO il = 1, ncum
1706      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1707        asij(il) = amax1(1.0E-16, asij(il))
1708        asij(il) = 1.0/asij(il)
1709        asum(il, i) = 0.0
1710        bsum(il, i) = 0.0
1711        csum(il, i) = 0.0
1712      END IF
1713    END DO
1714
1715    DO j = minorig, nl
1716      DO il = 1, ncum
1717        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1718            il)-1) .AND. j<=inb(il)) THEN
1719          ment(il, i, j) = ment(il, i, j)*asij(il)
1720        END IF
1721      END DO
1722    END DO
1723
1724    DO j = minorig, nl
1725      DO il = 1, ncum
1726        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1727            il)-1) .AND. j<=inb(il)) THEN
1728          asum(il, i) = asum(il, i) + ment(il, i, j)
1729          ment(il, i, j) = ment(il, i, j)*sig(il, j)
1730          bsum(il, i) = bsum(il, i) + ment(il, i, j)
1731        END IF
1732      END DO
1733    END DO
1734
1735    DO il = 1, ncum
1736      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1737        bsum(il, i) = amax1(bsum(il,i), 1.0E-16)
1738        bsum(il, i) = 1.0/bsum(il, i)
1739      END IF
1740    END DO
1741
1742    DO j = minorig, nl
1743      DO il = 1, ncum
1744        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1745            il)-1) .AND. j<=inb(il)) THEN
1746          ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i)
1747        END IF
1748      END DO
1749    END DO
1750
1751    DO j = minorig, nl
1752      DO il = 1, ncum
1753        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
1754            il)-1) .AND. j<=inb(il)) THEN
1755          csum(il, i) = csum(il, i) + ment(il, i, j)
1756        END IF
1757      END DO
1758    END DO
1759
1760    DO il = 1, ncum
1761      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1762          csum(il,i)<m(il,i)) THEN
1763        nent(il, i) = 0
1764        ment(il, i, i) = m(il, i)
1765        qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i)
1766        uent(il, i, i) = u(il, nk(il))
1767        vent(il, i, i) = v(il, nk(il))
1768        elij(il, i, i) = clw(il, i)
1769        ! MAF        sij(il,i,i)=1.0
1770        sij(il, i, i) = 0.0
1771      END IF
1772    END DO ! il
1773
1774    ! do j=1,ntra
1775    ! do il=1,ncum
1776    ! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1777    ! :     .and. csum(il,i).lt.m(il,i) ) then
1778    ! traent(il,i,i,j)=tra(il,nk(il),j)
1779    ! endif
1780    ! enddo
1781    ! enddo
1782789 END DO
1783
1784  ! MAF: renormalisation de MENT
1785  DO jm = 1, nd
1786    DO im = 1, nd
1787      DO il = 1, ncum
1788        zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
1789      END DO
1790    END DO
1791  END DO
1792
1793  DO jm = 1, nd
1794    DO im = 1, nd
1795      DO il = 1, ncum
1796        IF (zm(il,im)/=0.) THEN
1797          ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
1798        END IF
1799      END DO
1800    END DO
1801  END DO
1802
1803  DO jm = 1, nd
1804    DO im = 1, nd
1805      DO il = 1, ncum
1806        qents(il, im, jm) = qent(il, im, jm)
1807        ments(il, im, jm) = ment(il, im, jm)
1808      END DO
1809    END DO
1810  END DO
1811
1812  RETURN
1813END SUBROUTINE cv30_mixing
1814
1815
1816SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, &
1817    v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, &
1818    mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg
1819    , wdtraina, wdtrainm) ! 26/08/10  RomP-jyg
1820  IMPLICIT NONE
1821
1822
1823  include "cvthermo.h"
1824  include "cv30param.h"
1825  include "cvflag.h"
1826
1827  ! inputs:
1828  INTEGER ncum, nd, na, ntra, nloc
1829  INTEGER icb(nloc), inb(nloc)
1830  REAL delt, plcl(nloc)
1831  REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
1832  REAL u(nloc, nd), v(nloc, nd)
1833  REAL tra(nloc, nd, ntra)
1834  REAL p(nloc, nd), ph(nloc, nd+1)
1835  REAL th(nloc, na), gz(nloc, na)
1836  REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na)
1837  REAL cpn(nloc, na), tv(nloc, na)
1838  REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
1839
1840  ! outputs:
1841  REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
1842  REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
1843  REAL trap(nloc, na, ntra)
1844  REAL b(nloc, na)
1845  ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
1846  ! lascendance adiabatique et des flux melanges Pa et Pm.
1847  ! Distinction des wdtrain
1848  ! Pa = wdtrainA     Pm = wdtrainM
1849  REAL wdtraina(nloc, na), wdtrainm(nloc, na)
1850
1851  ! local variables
1852  INTEGER i, j, k, il, num1
1853  REAL tinv, delti
1854  REAL awat, afac, afac1, afac2, bfac
1855  REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth
1856  REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1857  REAL ampmax
1858  REAL lvcp(nloc, na)
1859  REAL wdtrain(nloc)
1860  LOGICAL lwork(nloc)
1861
1862
1863  ! ------------------------------------------------------
1864
1865  delti = 1./delt
1866  tinv = 1./3.
1867
1868  mp(:, :) = 0.
1869
1870  DO i = 1, nl
1871    DO il = 1, ncum
1872      mp(il, i) = 0.0
1873      rp(il, i) = rr(il, i)
1874      up(il, i) = u(il, i)
1875      vp(il, i) = v(il, i)
1876      wt(il, i) = 0.001
1877      water(il, i) = 0.0
1878      evap(il, i) = 0.0
1879      b(il, i) = 0.0
1880      lvcp(il, i) = lv(il, i)/cpn(il, i)
1881    END DO
1882  END DO
1883
1884  ! do k=1,ntra
1885  ! do i=1,nd
1886  ! do il=1,ncum
1887  ! trap(il,i,k)=tra(il,i,k)
1888  ! enddo
1889  ! enddo
1890  ! enddo
1891  ! ! RomP >>>
1892  DO i = 1, nd
1893    DO il = 1, ncum
1894      wdtraina(il, i) = 0.0
1895      wdtrainm(il, i) = 0.0
1896    END DO
1897  END DO
1898  ! ! RomP <<<
1899
1900  ! ***  check whether ep(inb)=0, if so, skip precipitating    ***
1901  ! ***             downdraft calculation                      ***
1902
1903
1904  DO il = 1, ncum
1905    lwork(il) = .TRUE.
1906    IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE.
1907  END DO
1908
1909  CALL zilch(wdtrain, ncum)
1910
1911  DO i = nl + 1, 1, -1
1912
1913    num1 = 0
1914    DO il = 1, ncum
1915      IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
1916    END DO
1917    IF (num1<=0) GO TO 400
1918
1919
1920    ! ***  integrate liquid water equation to find condensed water   ***
1921    ! ***                and condensed water flux                    ***
1922
1923
1924
1925    ! ***                    begin downdraft loop                    ***
1926
1927
1928
1929    ! ***              calculate detrained precipitation             ***
1930
1931    DO il = 1, ncum
1932      IF (i<=inb(il) .AND. lwork(il)) THEN
1933        IF (cvflag_grav) THEN
1934          wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
1935          wdtraina(il, i) = wdtrain(il)/grav !   Pa  26/08/10   RomP
1936        ELSE
1937          wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
1938          wdtraina(il, i) = wdtrain(il)/10. !   Pa  26/08/10   RomP
1939        END IF
1940      END IF
1941    END DO
1942
1943    IF (i>1) THEN
1944
1945      DO j = 1, i - 1
1946        DO il = 1, ncum
1947          IF (i<=inb(il) .AND. lwork(il)) THEN
1948            awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
1949            awat = amax1(awat, 0.0)
1950            IF (cvflag_grav) THEN
1951              wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
1952            ELSE
1953              wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
1954            END IF
1955          END IF
1956        END DO
1957      END DO
1958      DO il = 1, ncum
1959        IF (cvflag_grav) THEN
1960          wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) !   Pm  26/08/10   RomP
1961        ELSE
1962          wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) !   Pm  26/08/10   RomP
1963        END IF
1964      END DO
1965
1966    END IF
1967
1968
1969    ! ***    find rain water and evaporation using provisional   ***
1970    ! ***              estimates of rp(i)and rp(i-1)             ***
1971
1972
1973    DO il = 1, ncum
1974
1975      IF (i<=inb(il) .AND. lwork(il)) THEN
1976
1977        wt(il, i) = 45.0
1978
1979        IF (i<inb(il)) THEN
1980          rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
1981            i))+gz(il,i+1)-gz(il,i))/lv(il, i)
1982          rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
1983        END IF
1984        rp(il, i) = amax1(rp(il,i), 0.0)
1985        rp(il, i) = amin1(rp(il,i), rs(il,i))
1986        rp(il, inb(il)) = rr(il, inb(il))
1987
1988        IF (i==1) THEN
1989          afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
1990        ELSE
1991          rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
1992            i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
1993          rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
1994          rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
1995          rp(il, i-1) = amax1(rp(il,i-1), 0.0)
1996          afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i) &
1997            )
1998          afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
1999            (1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
2000          afac = 0.5*(afac1+afac2)
2001        END IF
2002        IF (i==inb(il)) afac = 0.0
2003        afac = amax1(afac, 0.0)
2004        bfac = 1./(sigd*wt(il,i))
2005
2006        ! jyg1
2007        ! cc        sigt=1.0
2008        ! cc        if(i.ge.icb)sigt=sigp(i)
2009        ! prise en compte de la variation progressive de sigt dans
2010        ! les couches icb et icb-1:
2011        ! pour plcl<ph(i+1), pr1=0 & pr2=1
2012        ! pour plcl>ph(i),   pr1=1 & pr2=0
2013        ! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2014        ! sur le nuage, et pr2 est la proportion sous la base du
2015        ! nuage.
2016        pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2017        pr1 = max(0., min(1.,pr1))
2018        pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2019        pr2 = max(0., min(1.,pr2))
2020        sigt = sigp(il, i)*pr1 + pr2
2021        ! jyg2
2022
2023        b6 = bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac
2024        c6 = water(il, i+1) + bfac*wdtrain(il) - 50.*sigd*bfac*(ph(il,i)-ph( &
2025          il,i+1))*evap(il, i+1)
2026        IF (c6>0.0) THEN
2027          revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
2028          evap(il, i) = sigt*afac*revap
2029          water(il, i) = revap*revap
2030        ELSE
2031          evap(il, i) = -evap(il, i+1) + 0.02*(wdtrain(il)+sigd*wt(il,i)* &
2032            water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1)))
2033        END IF
2034
2035        ! ***  calculate precipitating downdraft mass flux under     ***
2036        ! ***              hydrostatic approximation                 ***
2037
2038        IF (i/=1) THEN
2039
2040          tevap = amax1(0.0, evap(il,i))
2041          delth = amax1(0.001, (th(il,i)-th(il,i-1)))
2042          IF (cvflag_grav) THEN
2043            mp(il, i) = 100.*ginv*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/ &
2044              delth
2045          ELSE
2046            mp(il, i) = 10.*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth
2047          END IF
2048
2049          ! ***           if hydrostatic assumption fails,             ***
2050          ! ***   solve cubic difference equation for downdraft theta  ***
2051          ! ***  and mass flux from two simultaneous differential eqns ***
2052
2053          amfac = sigd*sigd*70.0*ph(il, i)*(p(il,i-1)-p(il,i))* &
2054            (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2055          amp2 = abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
2056          IF (amp2>(0.1*amfac)) THEN
2057            xf = 100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1))
2058            tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i)* &
2059              sigd*th(il,i))
2060            af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv
2061            bf = 2.*(tinv*mp(il,i+1))**3 + tinv*mp(il, i+1)*xf*tf + &
2062              50.*(p(il,i-1)-p(il,i))*xf*tevap
2063            fac2 = 1.0
2064            IF (bf<0.0) fac2 = -1.0
2065            bf = abs(bf)
2066            ur = 0.25*bf*bf - af*af*af*tinv*tinv*tinv
2067            IF (ur>=0.0) THEN
2068              sru = sqrt(ur)
2069              fac = 1.0
2070              IF ((0.5*bf-sru)<0.0) fac = -1.0
2071              mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
2072                fac*(abs(0.5*bf-sru))**tinv
2073            ELSE
2074              d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2075              IF (fac2<0.0) d = 3.14159 - d
2076              mp(il, i) = mp(il, i+1)*tinv + 2.*sqrt(af*tinv)*cos(d*tinv)
2077            END IF
2078            mp(il, i) = amax1(0.0, mp(il,i))
2079
2080            IF (cvflag_grav) THEN
2081              ! jyg : il y a vraisemblablement une erreur dans la ligne 2
2082              ! suivante:
2083              ! il faut diviser par (mp(il,i)*sigd*grav) et non par
2084              ! (mp(il,i)+sigd*0.1).
2085              ! Et il faut bien revoir les facteurs 100.
2086              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, &
2087                i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2088                )*sigd*th(il,i))
2089            ELSE
2090              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, &
2091                i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i &
2092                )*sigd*th(il,i))
2093            END IF
2094            b(il, i-1) = amax1(b(il,i-1), 0.0)
2095          END IF
2096
2097          ! ***         limit magnitude of mp(i) to meet cfl condition
2098          ! ***
2099
2100          ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
2101          amp2 = 2.0*(ph(il,i-1)-ph(il,i))*delti
2102          ampmax = amin1(ampmax, amp2)
2103          mp(il, i) = amin1(mp(il,i), ampmax)
2104
2105          ! ***      force mp to decrease linearly to zero
2106          ! ***
2107          ! ***       between cloud base and the surface
2108          ! ***
2109
2110          IF (p(il,i)>p(il,icb(il))) THEN
2111            mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ &
2112              (p(il,1)-p(il,icb(il)))
2113          END IF
2114
2115        END IF ! i.eq.1
2116
2117        ! ***       find mixing ratio of precipitating downdraft     ***
2118
2119
2120        IF (i/=inb(il)) THEN
2121
2122          rp(il, i) = rr(il, i)
2123
2124          IF (mp(il,i)>mp(il,i+1)) THEN
2125
2126            IF (cvflag_grav) THEN
2127              rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2128                rr(il, i)*(mp(il,i)-mp(il,i+1)) + 100.*ginv*0.5*sigd*(ph(il,i &
2129                )-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2130            ELSE
2131              rp(il, i) = rp(il, i+1)*mp(il, i+1) + &
2132                rr(il, i)*(mp(il,i)-mp(il,i+1)) + 5.*sigd*(ph(il,i)-ph(il,i+1 &
2133                ))*(evap(il,i+1)+evap(il,i))
2134            END IF
2135            rp(il, i) = rp(il, i)/mp(il, i)
2136            up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+ &
2137              1))
2138            up(il, i) = up(il, i)/mp(il, i)
2139            vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+ &
2140              1))
2141            vp(il, i) = vp(il, i)/mp(il, i)
2142
2143            ! do j=1,ntra
2144            ! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2145            ! testmaf     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2146            ! :            +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2147            ! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2148            ! end do
2149
2150          ELSE
2151
2152            IF (mp(il,i+1)>1.0E-16) THEN
2153              IF (cvflag_grav) THEN
2154                rp(il, i) = rp(il, i+1) + 100.*ginv*0.5*sigd*(ph(il,i)-ph(il, &
2155                  i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2156              ELSE
2157                rp(il, i) = rp(il, i+1) + 5.*sigd*(ph(il,i)-ph(il,i+1))*(evap &
2158                  (il,i+1)+evap(il,i))/mp(il, i+1)
2159              END IF
2160              up(il, i) = up(il, i+1)
2161              vp(il, i) = vp(il, i+1)
2162
2163              ! do j=1,ntra
2164              ! trap(il,i,j)=trap(il,i+1,j)
2165              ! end do
2166
2167            END IF
2168          END IF
2169          rp(il, i) = amin1(rp(il,i), rs(il,i))
2170          rp(il, i) = amax1(rp(il,i), 0.0)
2171
2172        END IF
2173      END IF
2174    END DO
2175
2176400 END DO
2177
2178  RETURN
2179END SUBROUTINE cv30_unsat
2180
2181SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, &
2182    tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, &
2183    wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, &
2184    tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, &
2185    mike, tls, tps, qcondc, wd)
2186  IMPLICIT NONE
2187
2188  include "cvthermo.h"
2189  include "cv30param.h"
2190  include "cvflag.h"
2191  include "conema3.h"
2192
2193  ! inputs:
2194  INTEGER ncum, nd, na, ntra, nloc
2195  INTEGER icb(nloc), inb(nloc)
2196  REAL delt
2197  REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
2198  REAL tra(nloc, nd, ntra), sig(nloc, nd)
2199  REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
2200  REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
2201  REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
2202  REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
2203  REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
2204  REAL water(nloc, na), evap(nloc, na), b(nloc, na)
2205  REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
2206  ! ym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2207  REAL vent(nloc, na, na), elij(nloc, na, na)
2208  INTEGER nent(nloc, na)
2209  REAL traent(nloc, na, na, ntra)
2210  REAL tv(nloc, nd), tvp(nloc, nd)
2211
2212  ! input/output:
2213  INTEGER iflag(nloc)
2214
2215  ! outputs:
2216  REAL precip(nloc)
2217  REAL vprecip(nloc, nd+1)
2218  REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
2219  REAL ftra(nloc, nd, ntra)
2220  REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
2221  REAL dnwd0(nloc, nd), mike(nloc, nd)
2222  REAL tls(nloc, nd), tps(nloc, nd)
2223  REAL qcondc(nloc, nd) ! cld
2224  REAL wd(nloc) ! gust
2225
2226  ! local variables:
2227  INTEGER i, k, il, n, j, num1
2228  REAL rat, awat, delti
2229  REAL ax, bx, cx, dx, ex
2230  REAL cpinv, rdcp, dpinv
2231  REAL lvcp(nloc, na), mke(nloc, na)
2232  REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
2233  ! !!      real up1(nloc), dn1(nloc)
2234  REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
2235  REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2236  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
2237  REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
2238
2239
2240  ! -------------------------------------------------------------
2241
2242  ! initialization:
2243
2244  delti = 1.0/delt
2245
2246  DO il = 1, ncum
2247    precip(il) = 0.0
2248    wd(il) = 0.0 ! gust
2249    vprecip(il, nd+1) = 0.
2250  END DO
2251
2252  DO i = 1, nd
2253    DO il = 1, ncum
2254      vprecip(il, i) = 0.0
2255      ft(il, i) = 0.0
2256      fr(il, i) = 0.0
2257      fu(il, i) = 0.0
2258      fv(il, i) = 0.0
2259      qcondc(il, i) = 0.0 ! cld
2260      qcond(il, i) = 0.0 ! cld
2261      nqcond(il, i) = 0.0 ! cld
2262    END DO
2263  END DO
2264
2265  ! do j=1,ntra
2266  ! do i=1,nd
2267  ! do il=1,ncum
2268  ! ftra(il,i,j)=0.0
2269  ! enddo
2270  ! enddo
2271  ! enddo
2272
2273  DO i = 1, nl
2274    DO il = 1, ncum
2275      lvcp(il, i) = lv(il, i)/cpn(il, i)
2276    END DO
2277  END DO
2278
2279
2280
2281  ! ***  calculate surface precipitation in mm/day     ***
2282
2283  DO il = 1, ncum
2284    IF (ep(il,inb(il))>=0.0001) THEN
2285      IF (cvflag_grav) THEN
2286        precip(il) = wt(il, 1)*sigd*water(il, 1)*86400.*1000./(rowl*grav)
2287      ELSE
2288        precip(il) = wt(il, 1)*sigd*water(il, 1)*8640.
2289      END IF
2290    END IF
2291  END DO
2292
2293  ! ***  CALCULATE VERTICAL PROFILE OF  PRECIPITATIONs IN kg/m2/s  ===
2294
2295  ! MAF rajout pour lessivage
2296  DO k = 1, nl
2297    DO il = 1, ncum
2298      IF (k<=inb(il)) THEN
2299        IF (cvflag_grav) THEN
2300          vprecip(il, k) = wt(il, k)*sigd*water(il, k)/grav
2301        ELSE
2302          vprecip(il, k) = wt(il, k)*sigd*water(il, k)/10.
2303        END IF
2304      END IF
2305    END DO
2306  END DO
2307
2308
2309  ! ***  Calculate downdraft velocity scale    ***
2310  ! ***  NE PAS UTILISER POUR L'INSTANT ***
2311
2312  ! !      do il=1,ncum
2313  ! !        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
2314  ! !     :                                  /(sigd*p(il,icb(il)))
2315  ! !      enddo
2316
2317
2318  ! ***  calculate tendencies of lowest level potential temperature  ***
2319  ! ***                      and mixing ratio                        ***
2320
2321  DO il = 1, ncum
2322    work(il) = 1.0/(ph(il,1)-ph(il,2))
2323    am(il) = 0.0
2324  END DO
2325
2326  DO k = 2, nl
2327    DO il = 1, ncum
2328      IF (k<=inb(il)) THEN
2329        am(il) = am(il) + m(il, k)
2330      END IF
2331    END DO
2332  END DO
2333
2334  DO il = 1, ncum
2335
2336    ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
2337    IF (cvflag_grav) THEN
2338      IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
2339      ft(il, 1) = 0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2340        1))/cpn(il,1))
2341    ELSE
2342      IF ((0.1*work(il)*am(il))>=delti) iflag(il) = 1 !consistency vect
2343      ft(il, 1) = 0.1*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, &
2344        1))/cpn(il,1))
2345    END IF
2346
2347    ft(il, 1) = ft(il, 1) - 0.5*lvcp(il, 1)*sigd*(evap(il,1)+evap(il,2))
2348
2349    IF (cvflag_grav) THEN
2350      ft(il, 1) = ft(il, 1) - 0.009*grav*sigd*mp(il, 2)*t(il, 1)*b(il, 1)* &
2351        work(il)
2352    ELSE
2353      ft(il, 1) = ft(il, 1) - 0.09*sigd*mp(il, 2)*t(il, 1)*b(il, 1)*work(il)
2354    END IF
2355
2356    ft(il, 1) = ft(il, 1) + 0.01*sigd*wt(il, 1)*(cl-cpd)*water(il, 2)*(t(il,2 &
2357      )-t(il,1))*work(il)/cpn(il, 1)
2358
2359    IF (cvflag_grav) THEN
2360      ! jyg1  Correction pour mieux conserver l'eau (conformite avec
2361      ! CONVECT4.3)
2362      ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas
2363      ! evap)
2364      fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + &
2365        sigd*0.5*(evap(il,1)+evap(il,2))
2366      ! +tard     :          +sigd*evap(il,1)
2367
2368      fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
2369
2370      fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il, &
2371        1))+am(il)*(u(il,2)-u(il,1)))
2372      fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
2373        1))+am(il)*(v(il,2)-v(il,1)))
2374    ELSE ! cvflag_grav
2375      fr(il, 1) = 0.1*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + &
2376        sigd*0.5*(evap(il,1)+evap(il,2))
2377      fr(il, 1) = fr(il, 1) + 0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
2378      fu(il, 1) = fu(il, 1) + 0.1*work(il)*(mp(il,2)*(up(il,2)-u(il, &
2379        1))+am(il)*(u(il,2)-u(il,1)))
2380      fv(il, 1) = fv(il, 1) + 0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
2381        1))+am(il)*(v(il,2)-v(il,1)))
2382    END IF ! cvflag_grav
2383
2384  END DO ! il
2385
2386  ! do j=1,ntra
2387  ! do il=1,ncum
2388  ! if (cvflag_grav) then
2389  ! ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
2390  ! :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2391  ! :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
2392  ! else
2393  ! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2394  ! :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2395  ! :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
2396  ! endif
2397  ! enddo
2398  ! enddo
2399
2400  DO j = 2, nl
2401    DO il = 1, ncum
2402      IF (j<=inb(il)) THEN
2403        IF (cvflag_grav) THEN
2404          fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il, &
2405            j,1)-rr(il,1))
2406          fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il, &
2407            j,1)-u(il,1))
2408          fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il, &
2409            j,1)-v(il,1))
2410        ELSE ! cvflag_grav
2411          fr(il, 1) = fr(il, 1) + 0.1*work(il)*ment(il, j, 1)*(qent(il,j,1)- &
2412            rr(il,1))
2413          fu(il, 1) = fu(il, 1) + 0.1*work(il)*ment(il, j, 1)*(uent(il,j,1)-u &
2414            (il,1))
2415          fv(il, 1) = fv(il, 1) + 0.1*work(il)*ment(il, j, 1)*(vent(il,j,1)-v &
2416            (il,1))
2417        END IF ! cvflag_grav
2418      END IF ! j
2419    END DO
2420  END DO
2421
2422  ! do k=1,ntra
2423  ! do j=2,nl
2424  ! do il=1,ncum
2425  ! if (j.le.inb(il)) then
2426
2427  ! if (cvflag_grav) then
2428  ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
2429  ! :                *(traent(il,j,1,k)-tra(il,1,k))
2430  ! else
2431  ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
2432  ! :                *(traent(il,j,1,k)-tra(il,1,k))
2433  ! endif
2434
2435  ! endif
2436  ! enddo
2437  ! enddo
2438  ! enddo
2439
2440
2441  ! ***  calculate tendencies of potential temperature and mixing ratio  ***
2442  ! ***               at levels above the lowest level                   ***
2443
2444  ! ***  first find the net saturated updraft and downdraft mass fluxes  ***
2445  ! ***                      through each level                          ***
2446
2447
2448  DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
2449
2450    num1 = 0
2451    DO il = 1, ncum
2452      IF (i<=inb(il)) num1 = num1 + 1
2453    END DO
2454    IF (num1<=0) GO TO 500
2455
2456    CALL zilch(amp1, ncum)
2457    CALL zilch(ad, ncum)
2458
2459    DO k = i + 1, nl + 1
2460      DO il = 1, ncum
2461        IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN
2462          amp1(il) = amp1(il) + m(il, k)
2463        END IF
2464      END DO
2465    END DO
2466
2467    DO k = 1, i
2468      DO j = i + 1, nl + 1
2469        DO il = 1, ncum
2470          IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
2471            amp1(il) = amp1(il) + ment(il, k, j)
2472          END IF
2473        END DO
2474      END DO
2475    END DO
2476
2477    DO k = 1, i - 1
2478      DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
2479        DO il = 1, ncum
2480          IF (i<=inb(il) .AND. j<=inb(il)) THEN
2481            ad(il) = ad(il) + ment(il, j, k)
2482          END IF
2483        END DO
2484      END DO
2485    END DO
2486
2487    DO il = 1, ncum
2488      IF (i<=inb(il)) THEN
2489        dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2490        cpinv = 1.0/cpn(il, i)
2491
2492        ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
2493        IF (cvflag_grav) THEN
2494          IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2495        ELSE
2496          IF ((0.1*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
2497        END IF
2498
2499        IF (cvflag_grav) THEN
2500          ft(il, i) = 0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
2501            i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
2502            i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( &
2503            il,i)+evap(il,i+1))
2504          rat = cpn(il, i-1)*cpinv
2505          ft(il, i) = ft(il, i) - 0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
2506            -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2507          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h( &
2508            il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2509        ELSE ! cvflag_grav
2510          ft(il, i) = 0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
2511            i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
2512            i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( &
2513            il,i)+evap(il,i+1))
2514          rat = cpn(il, i-1)*cpinv
2515          ft(il, i) = ft(il, i) - 0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i)-mp(il &
2516            ,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2517          ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i)+ &
2518            t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2519        END IF ! cvflag_grav
2520
2521
2522        ft(il, i) = ft(il, i) + 0.01*sigd*wt(il, i)*(cl-cpd)*water(il, i+1)*( &
2523          t(il,i+1)-t(il,i))*dpinv*cpinv
2524
2525        IF (cvflag_grav) THEN
2526          fr(il, i) = 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
2527            i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2528          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
2529            i))-ad(il)*(u(il,i)-u(il,i-1)))
2530          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
2531            i))-ad(il)*(v(il,i)-v(il,i-1)))
2532        ELSE ! cvflag_grav
2533          fr(il, i) = 0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
2534            i))-ad(il)*(rr(il,i)-rr(il,i-1)))
2535          fu(il, i) = fu(il, i) + 0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
2536            i))-ad(il)*(u(il,i)-u(il,i-1)))
2537          fv(il, i) = fv(il, i) + 0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
2538            i))-ad(il)*(v(il,i)-v(il,i-1)))
2539        END IF ! cvflag_grav
2540
2541      END IF ! i
2542    END DO
2543
2544    ! do k=1,ntra
2545    ! do il=1,ncum
2546    ! if (i.le.inb(il)) then
2547    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2548    ! cpinv=1.0/cpn(il,i)
2549    ! if (cvflag_grav) then
2550    ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
2551    ! :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2552    ! :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2553    ! else
2554    ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2555    ! :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2556    ! :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2557    ! endif
2558    ! endif
2559    ! enddo
2560    ! enddo
2561
2562    DO k = 1, i - 1
2563      DO il = 1, ncum
2564        IF (i<=inb(il)) THEN
2565          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2566          cpinv = 1.0/cpn(il, i)
2567
2568          awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
2569          awat = amax1(awat, 0.0)
2570
2571          IF (cvflag_grav) THEN
2572            fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
2573              ,i)-awat-rr(il,i))
2574            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2575              ,i)-u(il,i))
2576            fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
2577              ,i)-v(il,i))
2578          ELSE ! cvflag_grav
2579            fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)- &
2580              awat-rr(il,i))
2581            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2582              ,i)-u(il,i))
2583            fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
2584              il,i))
2585          END IF ! cvflag_grav
2586
2587          ! (saturated updrafts resulting from mixing)        ! cld
2588          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld
2589          nqcond(il, i) = nqcond(il, i) + 1. ! cld
2590        END IF ! i
2591      END DO
2592    END DO
2593
2594    ! do j=1,ntra
2595    ! do k=1,i-1
2596    ! do il=1,ncum
2597    ! if (i.le.inb(il)) then
2598    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2599    ! cpinv=1.0/cpn(il,i)
2600    ! if (cvflag_grav) then
2601    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2602    ! :        *(traent(il,k,i,j)-tra(il,i,j))
2603    ! else
2604    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2605    ! :        *(traent(il,k,i,j)-tra(il,i,j))
2606    ! endif
2607    ! endif
2608    ! enddo
2609    ! enddo
2610    ! enddo
2611
2612    DO k = i, nl + 1
2613      DO il = 1, ncum
2614        IF (i<=inb(il) .AND. k<=inb(il)) THEN
2615          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2616          cpinv = 1.0/cpn(il, i)
2617
2618          IF (cvflag_grav) THEN
2619            fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
2620              ,i)-rr(il,i))
2621            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
2622              ,i)-u(il,i))
2623            fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
2624              ,i)-v(il,i))
2625          ELSE ! cvflag_grav
2626            fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)-rr &
2627              (il,i))
2628            fu(il, i) = fu(il, i) + 0.1*dpinv*ment(il, k, i)*(uent(il,k,i)-u( &
2629              il,i))
2630            fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
2631              il,i))
2632          END IF ! cvflag_grav
2633        END IF ! i and k
2634      END DO
2635    END DO
2636
2637    ! do j=1,ntra
2638    ! do k=i,nl+1
2639    ! do il=1,ncum
2640    ! if (i.le.inb(il) .and. k.le.inb(il)) then
2641    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2642    ! cpinv=1.0/cpn(il,i)
2643    ! if (cvflag_grav) then
2644    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2645    ! :         *(traent(il,k,i,j)-tra(il,i,j))
2646    ! else
2647    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2648    ! :             *(traent(il,k,i,j)-tra(il,i,j))
2649    ! endif
2650    ! endif ! i and k
2651    ! enddo
2652    ! enddo
2653    ! enddo
2654
2655    DO il = 1, ncum
2656      IF (i<=inb(il)) THEN
2657        dpinv = 1.0/(ph(il,i)-ph(il,i+1))
2658        cpinv = 1.0/cpn(il, i)
2659
2660        IF (cvflag_grav) THEN
2661          ! sb: on ne fait pas encore la correction permettant de mieux
2662          ! conserver l'eau:
2663          fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + &
2664            0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il, &
2665            i)-rr(il,i-1)))*dpinv
2666
2667          fu(il, i) = fu(il, i) + 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il, &
2668            i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2669          fv(il, i) = fv(il, i) + 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il, &
2670            i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2671        ELSE ! cvflag_grav
2672          fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + &
2673            0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, &
2674            i-1)))*dpinv
2675          fu(il, i) = fu(il, i) + 0.1*(mp(il,i+1)*(up(il,i+1)-u(il, &
2676            i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2677          fv(il, i) = fv(il, i) + 0.1*(mp(il,i+1)*(vp(il,i+1)-v(il, &
2678            i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2679        END IF ! cvflag_grav
2680
2681      END IF ! i
2682    END DO
2683
2684    ! sb: interface with the cloud parameterization:          ! cld
2685
2686    DO k = i + 1, nl
2687      DO il = 1, ncum
2688        IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
2689          ! (saturated downdrafts resulting from mixing)            ! cld
2690          qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
2691          nqcond(il, i) = nqcond(il, i) + 1. ! cld
2692        END IF ! cld
2693      END DO ! cld
2694    END DO ! cld
2695
2696    ! (particular case: no detraining level is found)         ! cld
2697    DO il = 1, ncum ! cld
2698      IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld
2699        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
2700        nqcond(il, i) = nqcond(il, i) + 1. ! cld
2701      END IF ! cld
2702    END DO ! cld
2703
2704    DO il = 1, ncum ! cld
2705      IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld
2706        qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
2707      END IF ! cld
2708    END DO
2709
2710    ! do j=1,ntra
2711    ! do il=1,ncum
2712    ! if (i.le.inb(il)) then
2713    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2714    ! cpinv=1.0/cpn(il,i)
2715
2716    ! if (cvflag_grav) then
2717    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
2718    ! :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2719    ! :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2720    ! else
2721    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2722    ! :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2723    ! :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2724    ! endif
2725    ! endif ! i
2726    ! enddo
2727    ! enddo
2728
2729500 END DO
2730
2731
2732  ! ***   move the detrainment at level inb down to level inb-1   ***
2733  ! ***        in such a way as to preserve the vertically        ***
2734  ! ***          integrated enthalpy and water tendencies         ***
2735
2736  DO il = 1, ncum
2737
2738    ax = 0.1*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t(il, &
2739      inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
2740      inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2741    ft(il, inb(il)) = ft(il, inb(il)) - ax
2742    ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il &
2743      ))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, &
2744      inb(il))))
2745
2746    bx = 0.1*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb( &
2747      il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2748    fr(il, inb(il)) = fr(il, inb(il)) - bx
2749    fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ &
2750      1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2751
2752    cx = 0.1*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il &
2753      )))/(ph(il,inb(il))-ph(il,inb(il)+1))
2754    fu(il, inb(il)) = fu(il, inb(il)) - cx
2755    fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ &
2756      1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2757
2758    dx = 0.1*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il &
2759      )))/(ph(il,inb(il))-ph(il,inb(il)+1))
2760    fv(il, inb(il)) = fv(il, inb(il)) - dx
2761    fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ &
2762      1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
2763
2764  END DO
2765
2766  ! do j=1,ntra
2767  ! do il=1,ncum
2768  ! ex=0.1*ment(il,inb(il),inb(il))
2769  ! :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2770  ! :      /(ph(il,inb(il))-ph(il,inb(il)+1))
2771  ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2772  ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2773  ! :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2774  ! :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
2775  ! enddo
2776  ! enddo
2777
2778
2779  ! ***    homoginize tendencies below cloud base    ***
2780
2781
2782  DO il = 1, ncum
2783    asum(il) = 0.0
2784    bsum(il) = 0.0
2785    csum(il) = 0.0
2786    dsum(il) = 0.0
2787  END DO
2788
2789  DO i = 1, nl
2790    DO il = 1, ncum
2791      IF (i<=(icb(il)-1)) THEN
2792        asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1))
2793        bsum(il) = bsum(il) + fr(il, i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2794          1)))*(ph(il,i)-ph(il,i+1))
2795        csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
2796          1)))*(ph(il,i)-ph(il,i+1))
2797        dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
2798      END IF
2799    END DO
2800  END DO
2801
2802  ! !!!      do 700 i=1,icb(il)-1
2803  DO i = 1, nl
2804    DO il = 1, ncum
2805      IF (i<=(icb(il)-1)) THEN
2806        ft(il, i) = asum(il)*t(il, i)/(th(il,i)*dsum(il))
2807        fr(il, i) = bsum(il)/csum(il)
2808      END IF
2809    END DO
2810  END DO
2811
2812
2813  ! ***           reset counter and return           ***
2814
2815  DO il = 1, ncum
2816    sig(il, nd) = 2.0
2817  END DO
2818
2819
2820  DO i = 1, nd
2821    DO il = 1, ncum
2822      upwd(il, i) = 0.0
2823      dnwd(il, i) = 0.0
2824    END DO
2825  END DO
2826
2827  DO i = 1, nl
2828    DO il = 1, ncum
2829      dnwd0(il, i) = -mp(il, i)
2830    END DO
2831  END DO
2832  DO i = nl + 1, nd
2833    DO il = 1, ncum
2834      dnwd0(il, i) = 0.
2835    END DO
2836  END DO
2837
2838
2839  DO i = 1, nl
2840    DO il = 1, ncum
2841      IF (i>=icb(il) .AND. i<=inb(il)) THEN
2842        upwd(il, i) = 0.0
2843        dnwd(il, i) = 0.0
2844      END IF
2845    END DO
2846  END DO
2847
2848  DO i = 1, nl
2849    DO k = 1, nl
2850      DO il = 1, ncum
2851        up1(il, k, i) = 0.0
2852        dn1(il, k, i) = 0.0
2853      END DO
2854    END DO
2855  END DO
2856
2857  DO i = 1, nl
2858    DO k = i, nl
2859      DO n = 1, i - 1
2860        DO il = 1, ncum
2861          IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
2862            up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
2863            dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
2864          END IF
2865        END DO
2866      END DO
2867    END DO
2868  END DO
2869
2870  DO i = 2, nl
2871    DO k = i, nl
2872      DO il = 1, ncum
2873        ! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
2874        ! then
2875        IF (i<=inb(il) .AND. k<=inb(il)) THEN
2876          upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
2877          dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
2878        END IF
2879      END DO
2880    END DO
2881  END DO
2882
2883
2884  ! !!!      DO il=1,ncum
2885  ! !!!      do i=icb(il),inb(il)
2886  ! !!!
2887  ! !!!      upwd(il,i)=0.0
2888  ! !!!      dnwd(il,i)=0.0
2889  ! !!!      do k=i,inb(il)
2890  ! !!!      up1=0.0
2891  ! !!!      dn1=0.0
2892  ! !!!      do n=1,i-1
2893  ! !!!      up1=up1+ment(il,n,k)
2894  ! !!!      dn1=dn1-ment(il,k,n)
2895  ! !!!      enddo
2896  ! !!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
2897  ! !!!      dnwd(il,i)=dnwd(il,i)+dn1
2898  ! !!!      enddo
2899  ! !!!      enddo
2900  ! !!!
2901  ! !!!      ENDDO
2902
2903  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2904  ! determination de la variation de flux ascendant entre
2905  ! deux niveau non dilue mike
2906  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2907
2908  DO i = 1, nl
2909    DO il = 1, ncum
2910      mike(il, i) = m(il, i)
2911    END DO
2912  END DO
2913
2914  DO i = nl + 1, nd
2915    DO il = 1, ncum
2916      mike(il, i) = 0.
2917    END DO
2918  END DO
2919
2920  DO i = 1, nd
2921    DO il = 1, ncum
2922      ma(il, i) = 0
2923    END DO
2924  END DO
2925
2926  DO i = 1, nl
2927    DO j = i, nl
2928      DO il = 1, ncum
2929        ma(il, i) = ma(il, i) + m(il, j)
2930      END DO
2931    END DO
2932  END DO
2933
2934  DO i = nl + 1, nd
2935    DO il = 1, ncum
2936      ma(il, i) = 0.
2937    END DO
2938  END DO
2939
2940  DO i = 1, nl
2941    DO il = 1, ncum
2942      IF (i<=(icb(il)-1)) THEN
2943        ma(il, i) = 0
2944      END IF
2945    END DO
2946  END DO
2947
2948  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2949  ! icb represente de niveau ou se trouve la
2950  ! base du nuage , et inb le top du nuage
2951  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2952
2953  DO i = 1, nd
2954    DO il = 1, ncum
2955      mke(il, i) = upwd(il, i) + dnwd(il, i)
2956    END DO
2957  END DO
2958
2959  DO i = 1, nd
2960    DO il = 1, ncum
2961      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il, &
2962        i))+rr(il,i)*cpv)
2963      tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp
2964      tps(il, i) = tp(il, i)
2965    END DO
2966  END DO
2967
2968
2969  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
2970  ! ***           of condensed water         ***            ! cld
2971  ! ! cld
2972
2973  DO i = 1, nd ! cld
2974    DO il = 1, ncum ! cld
2975      mac(il, i) = 0.0 ! cld
2976      wa(il, i) = 0.0 ! cld
2977      siga(il, i) = 0.0 ! cld
2978      sax(il, i) = 0.0 ! cld
2979    END DO ! cld
2980  END DO ! cld
2981
2982  DO i = minorig, nl ! cld
2983    DO k = i + 1, nl + 1 ! cld
2984      DO il = 1, ncum ! cld
2985        IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld
2986          mac(il, i) = mac(il, i) + m(il, k) ! cld
2987        END IF ! cld
2988      END DO ! cld
2989    END DO ! cld
2990  END DO ! cld
2991
2992  DO i = 1, nl ! cld
2993    DO j = 1, i ! cld
2994      DO il = 1, ncum ! cld
2995        IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
2996            .AND. j>=icb(il)) THEN ! cld
2997          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
2998            *(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
2999        END IF ! cld
3000      END DO ! cld
3001    END DO ! cld
3002  END DO ! cld
3003
3004  DO i = 1, nl ! cld
3005    DO il = 1, ncum ! cld
3006      IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
3007          .AND. sax(il,i)>0.0) THEN ! cld
3008        wa(il, i) = sqrt(2.*sax(il,i)) ! cld
3009      END IF ! cld
3010    END DO ! cld
3011  END DO ! cld
3012
3013  DO i = 1, nl ! cld
3014    DO il = 1, ncum ! cld
3015      IF (wa(il,i)>0.0) &          ! cld
3016        siga(il, i) = mac(il, i)/wa(il, i) & ! cld
3017        *rrd*tvp(il, i)/p(il, i)/100./delta ! cld
3018      siga(il, i) = min(siga(il,i), 1.0) ! cld
3019      ! IM cf. FH
3020      IF (iflag_clw==0) THEN
3021        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) & ! cld
3022          +(1.-siga(il,i))*qcond(il, i) ! cld
3023      ELSE IF (iflag_clw==1) THEN
3024        qcondc(il, i) = qcond(il, i) ! cld
3025      END IF
3026
3027    END DO ! cld
3028  END DO ! cld
3029
3030  RETURN
3031END SUBROUTINE cv30_yield
3032
3033! !RomP >>>
3034SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
3035    d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
3036  IMPLICIT NONE
3037
3038  include "cv30param.h"
3039
3040  ! inputs:
3041  INTEGER ncum, nd, na, nloc, len
3042  REAL ment(nloc, na, na), sij(nloc, na, na)
3043  REAL clw(nloc, nd), elij(nloc, na, na)
3044  REAL ep(nloc, na)
3045  INTEGER icb(nloc), inb(nloc)
3046  REAL vprecip(nloc, nd+1)
3047  ! ouputs:
3048  REAL da(nloc, na), phi(nloc, na, na)
3049  REAL phi2(nloc, na, na)
3050  REAL d1a(nloc, na), dam(nloc, na)
3051  REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
3052  ! variables pour tracer dans precip de l'AA et des mel
3053  ! local variables:
3054  INTEGER i, j, k
3055  REAL epm(nloc, na, na)
3056
3057  ! variables d'Emanuel : du second indice au troisieme
3058  ! --->    tab(i,k,j) -> de l origine k a l arrivee j
3059  ! ment, sij, elij
3060  ! variables personnelles : du troisieme au second indice
3061  ! --->    tab(i,j,k) -> de k a j
3062  ! phi, phi2
3063
3064  ! initialisations
3065  DO j = 1, na
3066    DO i = 1, ncum
3067      da(i, j) = 0.
3068      d1a(i, j) = 0.
3069      dam(i, j) = 0.
3070      eplamm(i, j) = 0.
3071    END DO
3072  END DO
3073  DO k = 1, na
3074    DO j = 1, na
3075      DO i = 1, ncum
3076        epm(i, j, k) = 0.
3077        epmlmmm(i, j, k) = 0.
3078        phi(i, j, k) = 0.
3079        phi2(i, j, k) = 0.
3080      END DO
3081    END DO
3082  END DO
3083
3084  ! fraction deau condensee dans les melanges convertie en precip : epm
3085  ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3086  DO j = 1, na
3087    DO k = 1, j - 1
3088      DO i = 1, ncum
3089        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3090          ! !jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3091          epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
3092          ! !
3093          epm(i, j, k) = max(epm(i,j,k), 0.0)
3094        END IF
3095      END DO
3096    END DO
3097  END DO
3098
3099  DO j = 1, na
3100    DO k = 1, na
3101      DO i = 1, ncum
3102        IF (k>=icb(i) .AND. k<=inb(i)) THEN
3103          eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- &
3104            sij(i,j,k))
3105        END IF
3106      END DO
3107    END DO
3108  END DO
3109
3110  DO j = 1, na
3111    DO k = 1, j - 1
3112      DO i = 1, ncum
3113        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3114          epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
3115        END IF
3116      END DO
3117    END DO
3118  END DO
3119
3120  ! matrices pour calculer la tendance des concentrations dans cvltr.F90
3121  DO j = 1, na
3122    DO k = 1, na
3123      DO i = 1, ncum
3124        da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j)
3125        phi(i, j, k) = sij(i, k, j)*ment(i, k, j)
3126        d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sij(i,k,j))
3127      END DO
3128    END DO
3129  END DO
3130
3131  DO j = 1, na
3132    DO k = 1, j - 1
3133      DO i = 1, ncum
3134        dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, j, k)*(1.-ep(i,k))*(1.- &
3135          sij(i,k,j))
3136        phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
3137      END DO
3138    END DO
3139  END DO
3140
3141  RETURN
3142END SUBROUTINE cv30_tracer
3143! RomP <<<
3144
3145SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
3146    vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
3147    dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
3148    epmlmmm, eplamm, wdtraina, wdtrainm, iflag1, precip1, vprecip1, evap1, &
3149    ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
3150    dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
3151    elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1)
3152  IMPLICIT NONE
3153
3154  include "cv30param.h"
3155
3156  ! inputs:
3157  INTEGER len, ncum, nd, ntra, nloc
3158  INTEGER idcum(nloc)
3159  INTEGER iflag(nloc)
3160  INTEGER inb(nloc)
3161  REAL precip(nloc)
3162  REAL vprecip(nloc, nd+1), evap(nloc, nd)
3163  REAL ep(nloc, nd)
3164  REAL sig(nloc, nd), w0(nloc, nd)
3165  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3166  REAL ftra(nloc, nd, ntra)
3167  REAL ma(nloc, nd)
3168  REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3169  REAL qcondc(nloc, nd)
3170  REAL wd(nloc), cape(nloc)
3171  REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
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  ! RomP >>>
3196  REAL phi21(len, nd, nd)
3197  REAL d1a1(len, nd), dam1(len, nd)
3198  REAL wdtraina1(len, nd), wdtrainm1(len, nd)
3199  REAL sij1(len, nd, nd)
3200  REAL elij1(len, nd, nd), clw1(len, nd)
3201  REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
3202  ! RomP <<<
3203
3204  ! local variables:
3205  INTEGER i, k, j
3206
3207  DO i = 1, ncum
3208    precip1(idcum(i)) = precip(i)
3209    iflag1(idcum(i)) = iflag(i)
3210    wd1(idcum(i)) = wd(i)
3211    inb1(idcum(i)) = inb(i)
3212    cape1(idcum(i)) = cape(i)
3213  END DO
3214
3215  DO k = 1, nl
3216    DO i = 1, ncum
3217      vprecip1(idcum(i), k) = vprecip(i, k)
3218      evap1(idcum(i), k) = evap(i, k) !<<< RomP
3219      sig1(idcum(i), k) = sig(i, k)
3220      w01(idcum(i), k) = w0(i, k)
3221      ft1(idcum(i), k) = ft(i, k)
3222      fq1(idcum(i), k) = fq(i, k)
3223      fu1(idcum(i), k) = fu(i, k)
3224      fv1(idcum(i), k) = fv(i, k)
3225      ma1(idcum(i), k) = ma(i, k)
3226      upwd1(idcum(i), k) = upwd(i, k)
3227      dnwd1(idcum(i), k) = dnwd(i, k)
3228      dnwd01(idcum(i), k) = dnwd0(i, k)
3229      qcondc1(idcum(i), k) = qcondc(i, k)
3230      da1(idcum(i), k) = da(i, k)
3231      mp1(idcum(i), k) = mp(i, k)
3232      ! RomP >>>
3233      ep1(idcum(i), k) = ep(i, k)
3234      d1a1(idcum(i), k) = d1a(i, k)
3235      dam1(idcum(i), k) = dam(i, k)
3236      clw1(idcum(i), k) = clw(i, k)
3237      eplamm1(idcum(i), k) = eplamm(i, k)
3238      wdtraina1(idcum(i), k) = wdtraina(i, k)
3239      wdtrainm1(idcum(i), k) = wdtrainm(i, k)
3240      ! RomP <<<
3241    END DO
3242  END DO
3243
3244  DO i = 1, ncum
3245    sig1(idcum(i), nd) = sig(i, nd)
3246  END DO
3247
3248
3249  ! do 2100 j=1,ntra
3250  ! do 2110 k=1,nd ! oct3
3251  ! do 2120 i=1,ncum
3252  ! ftra1(idcum(i),k,j)=ftra(i,k,j)
3253  ! 2120     continue
3254  ! 2110    continue
3255  ! 2100   continue
3256  DO j = 1, nd
3257    DO k = 1, nd
3258      DO i = 1, ncum
3259        sij1(idcum(i), k, j) = sij(i, k, j)
3260        phi1(idcum(i), k, j) = phi(i, k, j)
3261        phi21(idcum(i), k, j) = phi2(i, k, j)
3262        elij1(idcum(i), k, j) = elij(i, k, j)
3263        epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
3264      END DO
3265    END DO
3266  END DO
3267
3268  RETURN
3269END SUBROUTINE cv30_uncompress
3270
Note: See TracBrowser for help on using the repository browser.