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

Last change on this file since 5133 was 5117, checked in by abarral, 5 months ago

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