source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_cv30.f90 @ 5441

Last change on this file since 5441 was 5159, checked in by abarral, 5 months ago

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