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

Last change on this file since 5111 was 5111, checked in by abarral, 4 months ago

Put abort_physic into a module
Remove -g option from makelmdz_fcm, since that option is linked to a header file that isn't included anywhere.
(lint) light lint on traversed files
  • 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.4 KB
Line 
1! $Id: cv30_routines.F90 5111 2024-07-24 10:17:33Z 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 print_control_mod, 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  ! endif
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  ! endif
1106  ! endif
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  ! endif
1133  ! endif
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    ! endif
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  ! endif
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    ! endif
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  ! endif
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
2406  ! if (cvflag_grav) then
2407  ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
2408  ! :                *(traent(il,j,1,k)-tra(il,1,k))
2409  ! else
2410  ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
2411  ! :                *(traent(il,j,1,k)-tra(il,1,k))
2412  ! endif
2413
2414  ! endif
2415  ! enddo
2416  ! enddo
2417  ! enddo
2418
2419
2420  ! ***  calculate tendencies of potential temperature and mixing ratio  ***
2421  ! ***               at levels above the lowest level                   ***
2422
2423  ! ***  first find the net saturated updraft and downdraft mass fluxes  ***
2424  ! ***                      through each level                          ***
2425
2426  DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
2427
2428    num1 = 0
2429    DO il = 1, ncum
2430      IF (i<=inb(il)) num1 = num1 + 1
2431    END DO
2432    IF (num1<=0) GO TO 500
2433
2434    CALL zilch(amp1, ncum)
2435    CALL zilch(ad, ncum)
2436
2437    DO k = i + 1, nl + 1
2438      DO il = 1, ncum
2439        IF (i<=inb(il) .AND. k<=(inb(il) + 1)) THEN
2440          amp1(il) = amp1(il) + m(il, k)
2441        END IF
2442      END DO
2443    END DO
2444
2445    DO k = 1, i
2446      DO j = i + 1, nl + 1
2447        DO il = 1, ncum
2448          IF (i<=inb(il) .AND. j<=(inb(il) + 1)) THEN
2449            amp1(il) = amp1(il) + ment(il, k, j)
2450          END IF
2451        END DO
2452      END DO
2453    END DO
2454
2455    DO k = 1, i - 1
2456      DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
2457        DO il = 1, ncum
2458          IF (i<=inb(il) .AND. j<=inb(il)) THEN
2459            ad(il) = ad(il) + ment(il, j, k)
2460          END IF
2461        END DO
2462      END DO
2463    END DO
2464
2465    DO il = 1, ncum
2466      IF (i<=inb(il)) THEN
2467        dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
2468        cpinv = 1.0 / cpn(il, i)
2469
2470        ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
2471        IF (cvflag_grav) THEN
2472          IF ((0.01 * grav * dpinv * amp1(il))>=delti) iflag(il) = 1 ! vecto
2473        ELSE
2474          IF ((0.1 * dpinv * amp1(il))>=delti) iflag(il) = 1 ! vecto
2475        END IF
2476
2477        IF (cvflag_grav) THEN
2478          ft(il, i) = 0.01 * grav * dpinv * (amp1(il) * (t(il, i + 1) - t(il, &
2479                  i) + (gz(il, i + 1) - gz(il, i)) * cpinv) - ad(il) * (t(il, i) - t(il, &
2480                  i - 1) + (gz(il, i) - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) * (evap(&
2481                  il, i) + evap(il, i + 1))
2482          rat = cpn(il, i - 1) * cpinv
2483          ft(il, i) = ft(il, i) - 0.009 * grav * sigd * (mp(il, i + 1) * t(il, i) * b(il, i) &
2484                  - mp(il, i) * t(il, i - 1) * rat * b(il, i - 1)) * dpinv
2485          ft(il, i) = ft(il, i) + 0.01 * grav * dpinv * ment(il, i, i) * (hp(il, i) - h(&
2486                  il, i) + t(il, i) * (cpv - cpd) * (rr(il, i) - qent(il, i, i))) * cpinv
2487        ELSE ! cvflag_grav
2488          ft(il, i) = 0.1 * dpinv * (amp1(il) * (t(il, i + 1) - t(il, &
2489                  i) + (gz(il, i + 1) - gz(il, i)) * cpinv) - ad(il) * (t(il, i) - t(il, &
2490                  i - 1) + (gz(il, i) - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) * (evap(&
2491                  il, i) + evap(il, i + 1))
2492          rat = cpn(il, i - 1) * cpinv
2493          ft(il, i) = ft(il, i) - 0.09 * sigd * (mp(il, i + 1) * t(il, i) * b(il, i) - mp(il &
2494                  , i) * t(il, i - 1) * rat * b(il, i - 1)) * dpinv
2495          ft(il, i) = ft(il, i) + 0.1 * dpinv * ment(il, i, i) * (hp(il, i) - h(il, i) + &
2496                  t(il, i) * (cpv - cpd) * (rr(il, i) - qent(il, i, i))) * cpinv
2497        END IF ! cvflag_grav
2498
2499        ft(il, i) = ft(il, i) + 0.01 * sigd * wt(il, i) * (cl - cpd) * water(il, i + 1) * (&
2500                t(il, i + 1) - t(il, i)) * dpinv * cpinv
2501
2502        IF (cvflag_grav) THEN
2503          fr(il, i) = 0.01 * grav * dpinv * (amp1(il) * (rr(il, i + 1) - rr(il, &
2504                  i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
2505          fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * (amp1(il) * (u(il, i + 1) - u(il, &
2506                  i)) - ad(il) * (u(il, i) - u(il, i - 1)))
2507          fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * (amp1(il) * (v(il, i + 1) - v(il, &
2508                  i)) - ad(il) * (v(il, i) - v(il, i - 1)))
2509        ELSE ! cvflag_grav
2510          fr(il, i) = 0.1 * dpinv * (amp1(il) * (rr(il, i + 1) - rr(il, &
2511                  i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
2512          fu(il, i) = fu(il, i) + 0.1 * dpinv * (amp1(il) * (u(il, i + 1) - u(il, &
2513                  i)) - ad(il) * (u(il, i) - u(il, i - 1)))
2514          fv(il, i) = fv(il, i) + 0.1 * dpinv * (amp1(il) * (v(il, i + 1) - v(il, &
2515                  i)) - ad(il) * (v(il, i) - v(il, i - 1)))
2516        END IF ! cvflag_grav
2517
2518      END IF ! i
2519    END DO
2520
2521    ! do k=1,ntra
2522    ! do il=1,ncum
2523    ! if (i.le.inb(il)) then
2524    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2525    ! cpinv=1.0/cpn(il,i)
2526    ! if (cvflag_grav) then
2527    ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
2528    ! :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2529    ! :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2530    ! else
2531    ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2532    ! :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2533    ! :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2534    ! endif
2535    ! endif
2536    ! enddo
2537    ! enddo
2538
2539    DO k = 1, i - 1
2540      DO il = 1, ncum
2541        IF (i<=inb(il)) THEN
2542          dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
2543          cpinv = 1.0 / cpn(il, i)
2544
2545          awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
2546          awat = amax1(awat, 0.0)
2547
2548          IF (cvflag_grav) THEN
2549            fr(il, i) = fr(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (qent(il, k &
2550                    , i) - awat - rr(il, i))
2551            fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (uent(il, k &
2552                    , i) - u(il, i))
2553            fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (vent(il, k &
2554                    , i) - v(il, i))
2555          ELSE ! cvflag_grav
2556            fr(il, i) = fr(il, i) + 0.1 * dpinv * ment(il, k, i) * (qent(il, k, i) - &
2557                    awat - rr(il, i))
2558            fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (uent(il, k &
2559                    , i) - u(il, i))
2560            fv(il, i) = fv(il, i) + 0.1 * dpinv * ment(il, k, i) * (vent(il, k, i) - v(&
2561                    il, i))
2562          END IF ! cvflag_grav
2563
2564          ! (saturated updrafts resulting from mixing)        ! cld
2565          qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat) ! cld
2566          nqcond(il, i) = nqcond(il, i) + 1. ! cld
2567        END IF ! i
2568      END DO
2569    END DO
2570
2571    ! do j=1,ntra
2572    ! do k=1,i-1
2573    ! do il=1,ncum
2574    ! if (i.le.inb(il)) then
2575    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2576    ! cpinv=1.0/cpn(il,i)
2577    ! if (cvflag_grav) then
2578    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2579    ! :        *(traent(il,k,i,j)-tra(il,i,j))
2580    ! else
2581    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2582    ! :        *(traent(il,k,i,j)-tra(il,i,j))
2583    ! endif
2584    ! endif
2585    ! enddo
2586    ! enddo
2587    ! enddo
2588
2589    DO k = i, nl + 1
2590      DO il = 1, ncum
2591        IF (i<=inb(il) .AND. k<=inb(il)) THEN
2592          dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
2593          cpinv = 1.0 / cpn(il, i)
2594
2595          IF (cvflag_grav) THEN
2596            fr(il, i) = fr(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (qent(il, k &
2597                    , i) - rr(il, i))
2598            fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (uent(il, k &
2599                    , i) - u(il, i))
2600            fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * ment(il, k, i) * (vent(il, k &
2601                    , i) - v(il, i))
2602          ELSE ! cvflag_grav
2603            fr(il, i) = fr(il, i) + 0.1 * dpinv * ment(il, k, i) * (qent(il, k, i) - rr &
2604                    (il, i))
2605            fu(il, i) = fu(il, i) + 0.1 * dpinv * ment(il, k, i) * (uent(il, k, i) - u(&
2606                    il, i))
2607            fv(il, i) = fv(il, i) + 0.1 * dpinv * ment(il, k, i) * (vent(il, k, i) - v(&
2608                    il, i))
2609          END IF ! cvflag_grav
2610        END IF ! i and k
2611      END DO
2612    END DO
2613
2614    ! do j=1,ntra
2615    ! do k=i,nl+1
2616    ! do il=1,ncum
2617    ! if (i.le.inb(il) .and. k.le.inb(il)) then
2618    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2619    ! cpinv=1.0/cpn(il,i)
2620    ! if (cvflag_grav) then
2621    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2622    ! :         *(traent(il,k,i,j)-tra(il,i,j))
2623    ! else
2624    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2625    ! :             *(traent(il,k,i,j)-tra(il,i,j))
2626    ! endif
2627    ! endif ! i and k
2628    ! enddo
2629    ! enddo
2630    ! enddo
2631
2632    DO il = 1, ncum
2633      IF (i<=inb(il)) THEN
2634        dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
2635        cpinv = 1.0 / cpn(il, i)
2636
2637        IF (cvflag_grav) THEN
2638          ! sb: on ne fait pas encore la correction permettant de mieux
2639          ! conserver l'eau:
2640          fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) + evap(il, i + 1)) + &
2641                  0.01 * grav * (mp(il, i + 1) * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, &
2642                          i) - rr(il, i - 1))) * dpinv
2643
2644          fu(il, i) = fu(il, i) + 0.01 * grav * (mp(il, i + 1) * (up(il, i + 1) - u(il, &
2645                  i)) - mp(il, i) * (up(il, i) - u(il, i - 1))) * dpinv
2646          fv(il, i) = fv(il, i) + 0.01 * grav * (mp(il, i + 1) * (vp(il, i + 1) - v(il, &
2647                  i)) - mp(il, i) * (vp(il, i) - v(il, i - 1))) * dpinv
2648        ELSE ! cvflag_grav
2649          fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) + evap(il, i + 1)) + &
2650                  0.1 * (mp(il, i + 1) * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) - rr(il, &
2651                          i - 1))) * dpinv
2652          fu(il, i) = fu(il, i) + 0.1 * (mp(il, i + 1) * (up(il, i + 1) - u(il, &
2653                  i)) - mp(il, i) * (up(il, i) - u(il, i - 1))) * dpinv
2654          fv(il, i) = fv(il, i) + 0.1 * (mp(il, i + 1) * (vp(il, i + 1) - v(il, &
2655                  i)) - mp(il, i) * (vp(il, i) - v(il, i - 1))) * dpinv
2656        END IF ! cvflag_grav
2657
2658      END IF ! i
2659    END DO
2660
2661    ! sb: interface with the cloud parameterization:          ! cld
2662
2663    DO k = i + 1, nl
2664      DO il = 1, ncum
2665        IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld
2666          ! (saturated downdrafts resulting from mixing)            ! cld
2667          qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
2668          nqcond(il, i) = nqcond(il, i) + 1. ! cld
2669        END IF ! cld
2670      END DO ! cld
2671    END DO ! cld
2672
2673    ! (particular case: no detraining level is found)         ! cld
2674    DO il = 1, ncum ! cld
2675      IF (i<=inb(il) .AND. nent(il, i)==0) THEN ! cld
2676        qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i) ! cld
2677        nqcond(il, i) = nqcond(il, i) + 1. ! cld
2678      END IF ! cld
2679    END DO ! cld
2680
2681    DO il = 1, ncum ! cld
2682      IF (i<=inb(il) .AND. nqcond(il, i)/=0.) THEN ! cld
2683        qcond(il, i) = qcond(il, i) / nqcond(il, i) ! cld
2684      END IF ! cld
2685    END DO
2686
2687    ! do j=1,ntra
2688    ! do il=1,ncum
2689    ! if (i.le.inb(il)) then
2690    ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
2691    ! cpinv=1.0/cpn(il,i)
2692
2693    ! if (cvflag_grav) then
2694    ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
2695    ! :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2696    ! :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2697    ! else
2698    ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2699    ! :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2700    ! :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2701    ! endif
2702    ! endif ! i
2703    ! enddo
2704    ! enddo
2705
2706  500 END DO
2707
2708
2709  ! ***   move the detrainment at level inb down to level inb-1   ***
2710  ! ***        in such a way as to preserve the vertically        ***
2711  ! ***          integrated enthalpy and water tendencies         ***
2712
2713  DO il = 1, ncum
2714
2715    ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) - h(il, inb(il)) + t(il, &
2716            inb(il)) * (cpv - cpd) * (rr(il, inb(il)) - qent(il, inb(il), &
2717            inb(il)))) / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1)))
2718    ft(il, inb(il)) = ft(il, inb(il)) - ax
2719    ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) * (ph(il, inb(il &
2720            )) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) * (ph(il, inb(il) - 1) - ph(il, &
2721            inb(il))))
2722
2723    bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) - rr(il, inb(&
2724            il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
2725    fr(il, inb(il)) = fr(il, inb(il)) - bx
2726    fr(il, inb(il) - 1) = fr(il, inb(il) - 1) + bx * (ph(il, inb(il)) - ph(il, inb(il) + &
2727            1)) / (ph(il, inb(il) - 1) - ph(il, inb(il)))
2728
2729    cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) - u(il, inb(il &
2730            ))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
2731    fu(il, inb(il)) = fu(il, inb(il)) - cx
2732    fu(il, inb(il) - 1) = fu(il, inb(il) - 1) + cx * (ph(il, inb(il)) - ph(il, inb(il) + &
2733            1)) / (ph(il, inb(il) - 1) - ph(il, inb(il)))
2734
2735    dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) - v(il, inb(il &
2736            ))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
2737    fv(il, inb(il)) = fv(il, inb(il)) - dx
2738    fv(il, inb(il) - 1) = fv(il, inb(il) - 1) + dx * (ph(il, inb(il)) - ph(il, inb(il) + &
2739            1)) / (ph(il, inb(il) - 1) - ph(il, inb(il)))
2740
2741  END DO
2742
2743  ! do j=1,ntra
2744  ! do il=1,ncum
2745  ! ex=0.1*ment(il,inb(il),inb(il))
2746  ! :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2747  ! :      /(ph(il,inb(il))-ph(il,inb(il)+1))
2748  ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2749  ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2750  ! :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2751  ! :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
2752  ! enddo
2753  ! enddo
2754
2755
2756  ! ***    homoginize tendencies below cloud base    ***
2757
2758  DO il = 1, ncum
2759    asum(il) = 0.0
2760    bsum(il) = 0.0
2761    csum(il) = 0.0
2762    dsum(il) = 0.0
2763  END DO
2764
2765  DO i = 1, nl
2766    DO il = 1, ncum
2767      IF (i<=(icb(il) - 1)) THEN
2768        asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
2769        bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (cl - cpd) * (t(il, i) - t(il, &
2770                1))) * (ph(il, i) - ph(il, i + 1))
2771        csum(il) = csum(il) + (lv(il, i) + (cl - cpd) * (t(il, i) - t(il, &
2772                1))) * (ph(il, i) - ph(il, i + 1))
2773        dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) / th(il, i)
2774      END IF
2775    END DO
2776  END DO
2777
2778  ! !!!      do 700 i=1,icb(il)-1
2779  DO i = 1, nl
2780    DO il = 1, ncum
2781      IF (i<=(icb(il) - 1)) THEN
2782        ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
2783        fr(il, i) = bsum(il) / csum(il)
2784      END IF
2785    END DO
2786  END DO
2787
2788
2789  ! ***           reset counter and return           ***
2790
2791  DO il = 1, ncum
2792    sig(il, nd) = 2.0
2793  END DO
2794
2795  DO i = 1, nd
2796    DO il = 1, ncum
2797      upwd(il, i) = 0.0
2798      dnwd(il, i) = 0.0
2799    END DO
2800  END DO
2801
2802  DO i = 1, nl
2803    DO il = 1, ncum
2804      dnwd0(il, i) = -mp(il, i)
2805    END DO
2806  END DO
2807  DO i = nl + 1, nd
2808    DO il = 1, ncum
2809      dnwd0(il, i) = 0.
2810    END DO
2811  END DO
2812
2813  DO i = 1, nl
2814    DO il = 1, ncum
2815      IF (i>=icb(il) .AND. i<=inb(il)) THEN
2816        upwd(il, i) = 0.0
2817        dnwd(il, i) = 0.0
2818      END IF
2819    END DO
2820  END DO
2821
2822  DO i = 1, nl
2823    DO k = 1, nl
2824      DO il = 1, ncum
2825        up1(il, k, i) = 0.0
2826        dn1(il, k, i) = 0.0
2827      END DO
2828    END DO
2829  END DO
2830
2831  DO i = 1, nl
2832    DO k = i, nl
2833      DO n = 1, i - 1
2834        DO il = 1, ncum
2835          IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
2836            up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
2837            dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
2838          END IF
2839        END DO
2840      END DO
2841    END DO
2842  END DO
2843
2844  DO i = 2, nl
2845    DO k = i, nl
2846      DO il = 1, ncum
2847        ! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
2848        ! then
2849        IF (i<=inb(il) .AND. k<=inb(il)) THEN
2850          upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
2851          dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
2852        END IF
2853      END DO
2854    END DO
2855  END DO
2856
2857
2858  ! !!!      DO il=1,ncum
2859  ! !!!      do i=icb(il),inb(il)
2860  ! !!!
2861  ! !!!      upwd(il,i)=0.0
2862  ! !!!      dnwd(il,i)=0.0
2863  ! !!!      do k=i,inb(il)
2864  ! !!!      up1=0.0
2865  ! !!!      dn1=0.0
2866  ! !!!      do n=1,i-1
2867  ! !!!      up1=up1+ment(il,n,k)
2868  ! !!!      dn1=dn1-ment(il,k,n)
2869  ! !!!      enddo
2870  ! !!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
2871  ! !!!      dnwd(il,i)=dnwd(il,i)+dn1
2872  ! !!!      enddo
2873  ! !!!      enddo
2874  ! !!!
2875  ! !!!      ENDDO
2876
2877  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2878  ! determination de la variation de flux ascendant entre
2879  ! deux niveau non dilue mike
2880  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2881
2882  DO i = 1, nl
2883    DO il = 1, ncum
2884      mike(il, i) = m(il, i)
2885    END DO
2886  END DO
2887
2888  DO i = nl + 1, nd
2889    DO il = 1, ncum
2890      mike(il, i) = 0.
2891    END DO
2892  END DO
2893
2894  DO i = 1, nd
2895    DO il = 1, ncum
2896      ma(il, i) = 0
2897    END DO
2898  END DO
2899
2900  DO i = 1, nl
2901    DO j = i, nl
2902      DO il = 1, ncum
2903        ma(il, i) = ma(il, i) + m(il, j)
2904      END DO
2905    END DO
2906  END DO
2907
2908  DO i = nl + 1, nd
2909    DO il = 1, ncum
2910      ma(il, i) = 0.
2911    END DO
2912  END DO
2913
2914  DO i = 1, nl
2915    DO il = 1, ncum
2916      IF (i<=(icb(il) - 1)) THEN
2917        ma(il, i) = 0
2918      END IF
2919    END DO
2920  END DO
2921
2922  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2923  ! icb represente de niveau ou se trouve la
2924  ! base du nuage , et inb le top du nuage
2925  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2926
2927  DO i = 1, nd
2928    DO il = 1, ncum
2929      mke(il, i) = upwd(il, i) + dnwd(il, i)
2930    END DO
2931  END DO
2932
2933  DO i = 1, nd
2934    DO il = 1, ncum
2935      rdcp = (rrd * (1. - rr(il, i)) - rr(il, i) * rrv) / (cpd * (1. - rr(il, &
2936              i)) + rr(il, i) * cpv)
2937      tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp
2938      tps(il, i) = tp(il, i)
2939    END DO
2940  END DO
2941
2942
2943  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
2944  ! ***           of condensed water         ***            ! cld
2945  ! ! cld
2946
2947  DO i = 1, nd ! cld
2948    DO il = 1, ncum ! cld
2949      mac(il, i) = 0.0 ! cld
2950      wa(il, i) = 0.0 ! cld
2951      siga(il, i) = 0.0 ! cld
2952      sax(il, i) = 0.0 ! cld
2953    END DO ! cld
2954  END DO ! cld
2955
2956  DO i = minorig, nl ! cld
2957    DO k = i + 1, nl + 1 ! cld
2958      DO il = 1, ncum ! cld
2959        IF (i<=inb(il) .AND. k<=(inb(il) + 1)) THEN ! cld
2960          mac(il, i) = mac(il, i) + m(il, k) ! cld
2961        END IF ! cld
2962      END DO ! cld
2963    END DO ! cld
2964  END DO ! cld
2965
2966  DO i = 1, nl ! cld
2967    DO j = 1, i ! cld
2968      DO il = 1, ncum ! cld
2969        IF (i>=icb(il) .AND. i<=(inb(il) - 1) & ! cld
2970                .AND. j>=icb(il)) THEN ! cld
2971          sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) & ! cld
2972                  * (ph(il, j) - ph(il, j + 1)) / p(il, j) ! cld
2973        END IF ! cld
2974      END DO ! cld
2975    END DO ! cld
2976  END DO ! cld
2977
2978  DO i = 1, nl ! cld
2979    DO il = 1, ncum ! cld
2980      IF (i>=icb(il) .AND. i<=(inb(il) - 1) & ! cld
2981              .AND. sax(il, i)>0.0) THEN ! cld
2982        wa(il, i) = sqrt(2. * sax(il, i)) ! cld
2983      END IF ! cld
2984    END DO ! cld
2985  END DO ! cld
2986
2987  DO i = 1, nl ! cld
2988    DO il = 1, ncum ! cld
2989      IF (wa(il, i)>0.0) &          ! cld
2990              siga(il, i) = mac(il, i) / wa(il, i) & ! cld
2991                      * rrd * tvp(il, i) / p(il, i) / 100. / delta ! cld
2992      siga(il, i) = min(siga(il, i), 1.0) ! cld
2993      ! IM cf. FH
2994      IF (iflag_clw==0) THEN
2995        qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) & ! cld
2996                + (1. - siga(il, i)) * qcond(il, i) ! cld
2997      ELSE IF (iflag_clw==1) THEN
2998        qcondc(il, i) = qcond(il, i) ! cld
2999      END IF
3000
3001    END DO ! cld
3002  END DO ! cld
3003
3004END SUBROUTINE cv30_yield
3005
3006! !RomP >>>
3007SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, &
3008        d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
3009  IMPLICIT NONE
3010
3011  include "cv30param.h"
3012
3013  ! inputs:
3014  INTEGER ncum, nd, na, nloc, len
3015  REAL ment(nloc, na, na), sij(nloc, na, na)
3016  REAL clw(nloc, nd), elij(nloc, na, na)
3017  REAL ep(nloc, na)
3018  INTEGER icb(nloc), inb(nloc)
3019  REAL vprecip(nloc, nd + 1)
3020  ! ouputs:
3021  REAL da(nloc, na), phi(nloc, na, na)
3022  REAL phi2(nloc, na, na)
3023  REAL d1a(nloc, na), dam(nloc, na)
3024  REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
3025  ! variables pour tracer dans precip de l'AA et des mel
3026  ! local variables:
3027  INTEGER i, j, k, nam1
3028  REAL epm(nloc, na, na)
3029
3030  nam1 = na - 1 ! Introduced because ep is not defined for j=na
3031  ! variables d'Emanuel : du second indice au troisieme
3032  ! --->    tab(i,k,j) -> de l origine k a l arrivee j
3033  ! ment, sij, elij
3034  ! variables personnelles : du troisieme au second indice
3035  ! --->    tab(i,j,k) -> de k a j
3036  ! phi, phi2
3037
3038  ! initialisations
3039  DO j = 1, na
3040    DO i = 1, ncum
3041      da(i, j) = 0.
3042      d1a(i, j) = 0.
3043      dam(i, j) = 0.
3044      eplamm(i, j) = 0.
3045    END DO
3046  END DO
3047  DO k = 1, na
3048    DO j = 1, na
3049      DO i = 1, ncum
3050        epm(i, j, k) = 0.
3051        epmlmmm(i, j, k) = 0.
3052        phi(i, j, k) = 0.
3053        phi2(i, j, k) = 0.
3054      END DO
3055    END DO
3056  END DO
3057
3058  ! fraction deau condensee dans les melanges convertie en precip : epm
3059  ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3060  DO j = 1, nam1
3061    DO k = 1, j - 1
3062      DO i = 1, ncum
3063        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3064          ! !jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3065          epm(i, j, k) = 1. - (1. - ep(i, j)) * clw(i, j) / max(elij(i, k, j), 1.E-16)
3066          ! !
3067          epm(i, j, k) = max(epm(i, j, k), 0.0)
3068        END IF
3069      END DO
3070    END DO
3071  END DO
3072
3073  DO j = 1, nam1
3074    DO k = 1, nam1
3075      DO i = 1, ncum
3076        IF (k>=icb(i) .AND. k<=inb(i)) THEN
3077          eplamm(i, j) = eplamm(i, j) + ep(i, j) * clw(i, j) * ment(i, j, k) * (1. - &
3078                  sij(i, j, k))
3079        END IF
3080      END DO
3081    END DO
3082  END DO
3083
3084  DO j = 1, nam1
3085    DO k = 1, j - 1
3086      DO i = 1, ncum
3087        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
3088          epmlmmm(i, j, k) = epm(i, j, k) * elij(i, k, j) * ment(i, k, j)
3089        END IF
3090      END DO
3091    END DO
3092  END DO
3093
3094  ! matrices pour calculer la tendance des concentrations dans cvltr.F90
3095  DO j = 1, nam1
3096    DO k = 1, nam1
3097      DO i = 1, ncum
3098        da(i, j) = da(i, j) + (1. - sij(i, k, j)) * ment(i, k, j)
3099        phi(i, j, k) = sij(i, k, j) * ment(i, k, j)
3100        d1a(i, j) = d1a(i, j) + ment(i, k, j) * ep(i, k) * (1. - sij(i, k, j))
3101      END DO
3102    END DO
3103  END DO
3104
3105  DO j = 1, nam1
3106    DO k = 1, j - 1
3107      DO i = 1, ncum
3108        dam(i, j) = dam(i, j) + ment(i, k, j) * epm(i, j, k) * (1. - ep(i, k)) * (1. - &
3109                sij(i, k, j))
3110        phi2(i, j, k) = phi(i, j, k) * epm(i, j, k)
3111      END DO
3112    END DO
3113  END DO
3114
3115END SUBROUTINE cv30_tracer
3116! RomP <<<
3117
3118SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
3119        vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, &
3120        dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, &
3121        epmlmmm, eplamm, wdtraina, wdtrainm, epmax_diag, iflag1, precip1, vprecip1, evap1, &
3122        ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, &
3123        dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, &
3124        elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1, epmax_diag1) ! epmax_cape
3125  IMPLICIT NONE
3126
3127  include "cv30param.h"
3128
3129  ! inputs:
3130  INTEGER len, ncum, nd, ntra, nloc
3131  INTEGER idcum(nloc)
3132  INTEGER iflag(nloc)
3133  INTEGER inb(nloc)
3134  REAL precip(nloc)
3135  REAL vprecip(nloc, nd + 1), evap(nloc, nd)
3136  REAL ep(nloc, nd)
3137  REAL sig(nloc, nd), w0(nloc, nd)
3138  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
3139  REAL ftra(nloc, nd, ntra)
3140  REAL ma(nloc, nd)
3141  REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
3142  REAL qcondc(nloc, nd)
3143  REAL wd(nloc), cape(nloc)
3144  REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd)
3145  REAL epmax_diag(nloc) ! epmax_cape
3146  ! RomP >>>
3147  REAL phi2(nloc, nd, nd)
3148  REAL d1a(nloc, nd), dam(nloc, nd)
3149  REAL wdtraina(nloc, nd), wdtrainm(nloc, nd)
3150  REAL sij(nloc, nd, nd)
3151  REAL elij(nloc, nd, nd), clw(nloc, nd)
3152  REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd)
3153  ! RomP <<<
3154
3155  ! outputs:
3156  INTEGER iflag1(len)
3157  INTEGER inb1(len)
3158  REAL precip1(len)
3159  REAL vprecip1(len, nd + 1), evap1(len, nd) !<<< RomP
3160  REAL ep1(len, nd) !<<< RomP
3161  REAL sig1(len, nd), w01(len, nd)
3162  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
3163  REAL ftra1(len, nd, ntra)
3164  REAL ma1(len, nd)
3165  REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
3166  REAL qcondc1(nloc, nd)
3167  REAL wd1(nloc), cape1(nloc)
3168  REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd)
3169  REAL epmax_diag1(len) ! epmax_cape
3170  ! RomP >>>
3171  REAL phi21(len, nd, nd)
3172  REAL d1a1(len, nd), dam1(len, nd)
3173  REAL wdtraina1(len, nd), wdtrainm1(len, nd)
3174  REAL sij1(len, nd, nd)
3175  REAL elij1(len, nd, nd), clw1(len, nd)
3176  REAL epmlmmm1(len, nd, nd), eplamm1(len, nd)
3177  ! RomP <<<
3178
3179  ! local variables:
3180  INTEGER i, k, j
3181
3182  DO i = 1, ncum
3183    precip1(idcum(i)) = precip(i)
3184    iflag1(idcum(i)) = iflag(i)
3185    wd1(idcum(i)) = wd(i)
3186    inb1(idcum(i)) = inb(i)
3187    cape1(idcum(i)) = cape(i)
3188    epmax_diag1(idcum(i)) = epmax_diag(i) ! epmax_cape
3189  END DO
3190
3191  DO k = 1, nl
3192    DO i = 1, ncum
3193      vprecip1(idcum(i), k) = vprecip(i, k)
3194      evap1(idcum(i), k) = evap(i, k) !<<< RomP
3195      sig1(idcum(i), k) = sig(i, k)
3196      w01(idcum(i), k) = w0(i, k)
3197      ft1(idcum(i), k) = ft(i, k)
3198      fq1(idcum(i), k) = fq(i, k)
3199      fu1(idcum(i), k) = fu(i, k)
3200      fv1(idcum(i), k) = fv(i, k)
3201      ma1(idcum(i), k) = ma(i, k)
3202      upwd1(idcum(i), k) = upwd(i, k)
3203      dnwd1(idcum(i), k) = dnwd(i, k)
3204      dnwd01(idcum(i), k) = dnwd0(i, k)
3205      qcondc1(idcum(i), k) = qcondc(i, k)
3206      da1(idcum(i), k) = da(i, k)
3207      mp1(idcum(i), k) = mp(i, k)
3208      ! RomP >>>
3209      ep1(idcum(i), k) = ep(i, k)
3210      d1a1(idcum(i), k) = d1a(i, k)
3211      dam1(idcum(i), k) = dam(i, k)
3212      clw1(idcum(i), k) = clw(i, k)
3213      eplamm1(idcum(i), k) = eplamm(i, k)
3214      wdtraina1(idcum(i), k) = wdtraina(i, k)
3215      wdtrainm1(idcum(i), k) = wdtrainm(i, k)
3216      ! RomP <<<
3217    END DO
3218  END DO
3219
3220  DO i = 1, ncum
3221    sig1(idcum(i), nd) = sig(i, nd)
3222  END DO
3223
3224
3225  ! do 2100 j=1,ntra
3226  ! do 2110 k=1,nd ! oct3
3227  ! do 2120 i=1,ncum
3228  ! ftra1(idcum(i),k,j)=ftra(i,k,j)
3229  ! 2120     continue
3230  ! 2110    continue
3231  ! 2100   continue
3232  DO j = 1, nd
3233    DO k = 1, nd
3234      DO i = 1, ncum
3235        sij1(idcum(i), k, j) = sij(i, k, j)
3236        phi1(idcum(i), k, j) = phi(i, k, j)
3237        phi21(idcum(i), k, j) = phi2(i, k, j)
3238        elij1(idcum(i), k, j) = elij(i, k, j)
3239        epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j)
3240      END DO
3241    END DO
3242  END DO
3243
3244END SUBROUTINE cv30_uncompress
3245
3246SUBROUTINE cv30_epmax_fn_cape(nloc, ncum, nd &
3247        , cape, ep, hp, icb, inb, clw, nk, t, h, lv &
3248        , epmax_diag)
3249  USE lmdz_abort_physic, ONLY: abort_physic
3250  implicit none
3251
3252  ! On fait varier epmax en fn de la cape
3253  ! Il faut donc recalculer ep, et hp qui a déjà été calculé et
3254  ! qui en dépend
3255  ! Toutes les autres variables fn de ep sont calculées plus bas.
3256
3257  INCLUDE "cvthermo.h"
3258  INCLUDE "cv30param.h"
3259  INCLUDE "conema3.h"
3260
3261  ! inputs:
3262  integer ncum, nd, nloc
3263  integer icb(nloc), inb(nloc)
3264  real cape(nloc)
3265  real clw(nloc, nd), lv(nloc, nd), t(nloc, nd), h(nloc, nd)
3266  integer nk(nloc)
3267  ! inouts:
3268  real ep(nloc, nd)
3269  real hp(nloc, nd)
3270  ! outputs ou local
3271  real epmax_diag(nloc)
3272  ! locals
3273  integer i, k
3274  real hp_bak(nloc, nd)
3275  CHARACTER (LEN = 20) :: modname = 'cv30_epmax_fn_cape'
3276  CHARACTER (LEN = 80) :: abort_message
3277
3278  ! on recalcule ep et hp
3279
3280  if (coef_epmax_cape>1e-12) then
3281    do i = 1, ncum
3282      epmax_diag(i) = epmax - coef_epmax_cape * sqrt(cape(i))
3283      do k = 1, nl
3284        ep(i, k) = ep(i, k) / epmax * epmax_diag(i)
3285        ep(i, k) = amax1(ep(i, k), 0.0)
3286        ep(i, k) = amin1(ep(i, k), epmax_diag(i))
3287      enddo
3288    enddo
3289
3290    ! On recalcule hp:
3291    do k = 1, nl
3292      do i = 1, ncum
3293        hp_bak(i, k) = hp(i, k)
3294      enddo
3295    enddo
3296    do k = 1, nlp
3297      do i = 1, ncum
3298        hp(i, k) = h(i, k)
3299      enddo
3300    enddo
3301    do k = minorig + 1, nl
3302      do i = 1, ncum
3303        if((k>=icb(i)).and.(k<=inb(i)))then
3304          hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k)
3305        endif
3306      enddo
3307    enddo !do k=minorig+1,n
3308    !     write(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20)
3309    do i = 1, ncum
3310      do k = 1, nl
3311        if (abs(hp_bak(i, k) - hp(i, k))>0.01) then
3312          write(*, *) 'i,k=', i, k
3313          write(*, *) 'coef_epmax_cape=', coef_epmax_cape
3314          write(*, *) 'epmax_diag(i)=', epmax_diag(i)
3315          write(*, *) 'ep(i,k)=', ep(i, k)
3316          write(*, *) 'hp(i,k)=', hp(i, k)
3317          write(*, *) 'hp_bak(i,k)=', hp_bak(i, k)
3318          write(*, *) 'h(i,k)=', h(i, k)
3319          write(*, *) 'nk(i)=', nk(i)
3320          write(*, *) 'h(i,nk(i))=', h(i, nk(i))
3321          write(*, *) 'lv(i,k)=', lv(i, k)
3322          write(*, *) 't(i,k)=', t(i, k)
3323          write(*, *) 'clw(i,k)=', clw(i, k)
3324          write(*, *) 'cpd,cpv=', cpd, cpv
3325          CALL abort_physic(modname, abort_message, 1)
3326        endif
3327      enddo !do k=1,nl
3328    enddo !do i=1,ncum
3329  endif !if (coef_epmax_cape.gt.1e-12) then
3330
3331END SUBROUTINE  cv30_epmax_fn_cape
3332
3333
Note: See TracBrowser for help on using the repository browser.