source: LMDZ6/branches/contrails/libf/phylmd/cv30_routines_mod.f90 @ 5456

Last change on this file since 5456 was 5346, checked in by fhourdin, 4 weeks ago

Debut de replaysation de la convection profonde.

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