source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cv30_routines.F90 @ 3831

Last change on this file since 3831 was 3831, checked in by ymipsl, 10 years ago

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