source: LMDZ5/branches/LMDZ5_AR5/libf/phylmd/cv3_routines.F @ 1634

Last change on this file since 1634 was 1554, checked in by Laurent Fairhead, 13 years ago

Changes in some default valeus of parameters


Modifications des valeurs par défaut de certains paramètres

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