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