source: LMDZ4/trunk/libf/phylmd/cv3_routines.F @ 862

Last change on this file since 862 was 829, checked in by lmdzadmin, 17 years ago

Correction bug indices buoy dans cv3_undilute2
CRio/JYG/IM

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