source: LMDZ6/branches/Amaury_dev/libf/phylmd/cv30_routines.F90 @ 5140

Last change on this file since 5140 was 5140, checked in by abarral, 3 months ago

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