source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/cv3_routines.F90 @ 3069

Last change on this file since 3069 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

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