source: LMDZ.3.3/trunk/libf/phylmd/convect3.F @ 5306

Last change on this file since 5306 was 207, checked in by lmdz, 24 years ago

petit detail
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.7 KB
RevLine 
[207]1c $Header$
[199]2      SUBROUTINE CONVECT
3     *    (DTIME,  T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
4     *     ND,       NDP1,     NL, NTRA,  DELT,  IFLAG,
5     *     FT, FR, FU,  FV,  FTRA,  PRECIP,
6     *     icb,inb,   upwd,dnwd,dnwd0,SIG, W0,Mike,Mke,
7     *     Ma,MENTS,QENTS,TPS,TLS,SIGIJ,CAPE,TVP,PBASE,BUOYBASE,
8     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
9C
10C    ***  THE PARAMETER NA SHOULD IN GENERAL EQUAL ND   ***
11C
12c#################################################################
13cFleur       Introduction des traceurs dans convect3 le 6 juin 200
14c#################################################################
15#include "dimensions.h"
16#include "dimphy.h"
17      PARAMETER (NA=60)
18
19      integer NTRAC
20      PARAMETER (NTRAC=nqmx-2)
21
22      INTEGER NENT(NA)
23      REAL T1(ND),R1(ND),RS(ND),U(ND),V(ND),TRA(ND,NTRA)
24      REAL P(ND),PH(NDP1)
25      REAL FT(ND),FR(ND),FU(ND),FV(ND),FTRA(ND,NTRA)
26      REAL SIG(ND),W0(ND)
27      REAL UENT(NA,NA),VENT(NA,NA),TRAENT(NA,NA,NTRAC),TRATM(NA)
28      REAL UP(NA),VP(NA),TRAP(NA,NTRAC)
29      REAL M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA)
30      REAL SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA)
31      REAL RP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA)
32      REAL SIGP(NA),B(NA),TP(NA),CPN(NA)
33      REAL LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA)
34      REAL T(NA),RR(NA)
35C
36      REAL BUOY(NA)     !  Lifted parcel buoyancy
37      REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
38C                                      temperature wrt T1 and Q1
39      REAL CLW_NEW(NA),QI(NA)
40C
41      LOGICAL ICE_CONV
42      PARAMETER (ICE_CONV=.TRUE.)
43 
44cccccccccccccccccccccccccccccccccccccccccccccc
45c     declaration des variables a sortir
46ccccccccccccccccccccccccccccccccccccccccccccc
47      real Mke(nd)
48      real Mike(nd)
49      real Ma(nd)
50      real TPS(ND) !temperature dans les ascendances non diluees
51      real TLS(ND) !temperature potentielle
52      real MENTS(nd,nd)
53      real QENTS(nd,nd)
54      REAL SIGIJ(KLEV,KLEV)
55      REAL PBASE ! pressure at the cloud base level
56      REAL BUOYBASE ! buoyancy at the cloud base level
57cccccccccccccccccccccccccccccccccccccccccccccc
58 
59 
60 
61c
62      real dnwd0(nd)  !  precipitation driven unsaturated downdraft flux
63      real dnwd(nd), dn1  ! in-cloud saturated downdraft mass flux
64      real upwd(nd), up1  ! in-cloud saturated updraft mass flux
65C
66C   ***         ASSIGN VALUES OF THERMODYNAMIC CONSTANTS        ***
67C   ***             THESE SHOULD BE CONSISTENT WITH             ***
68C   ***              THOSE USED IN CALLING PROGRAM              ***
69C   ***     NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT  ***
70C
71c sb      CPD=1005.7
72c sb      CPV=1870.0
73c sb      CL=4190.0
74c sb      CPVMCL=CL-CPV
75c sb      RV=461.5
76c sb      RD=287.04
77c sb      EPS=RD/RV
78c sb      ALV0=2.501E6
79ccccccccccccccccccccccc
80c constantes coherentes avec le modele du Centre Europeen
81c sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
82c sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
83c sb      CPD = 3.5 * RD
84c sb      CPV = 4.0 * RV
85c sb      CL = 4218.0
86c sb      CPVMCL=CL-CPV
87c sb      EPS=RD/RV
88c sb      ALV0=2.5008E+06
89cccccccccccccccccccccc
90c on utilise les constantes thermo du Centre Europeen: (SB)
91c
92#include "YOMCST.h"
93c
94       CPD = RCPD
95       CPV = RCPV
96       CL = RCW
97       CPVMCL = CL-CPV
98       EPS = RD/RV
99       ALV0 = RLVTT
100c
101       NK = 1 ! origin level of the lifted parcel
102c
103cccccccccccccccccccccc
104C
105C           ***  INITIALIZE OUTPUT ARRAYS AND PARAMETERS  ***
106C
107      DO 5 I=1,ND
108         FT(I)=0.0
109         FR(I)=0.0
110         FU(I)=0.0
111         FV(I)=0.0
112         DO 4 J=1,NTRA
113          FTRA(I,J)=0.0
114    4    CONTINUE
115         T(I)=T1(I)
116         RR(I)=R1(I)
117    5 CONTINUE
118      DO 7 I=1,NL
119         RDCP=(RD*(1.-RR(I))+RR(I)*RV)/ (CPD*(1.-RR(I))+RR(I)*CPV)
120         TH(I)=T(I)*(1000.0/P(I))**RDCP
121    7 CONTINUE
122C
123*************************************************************
124**    CALCUL DES TEMPERATURES POTENTIELLES A SORTIR
125*************************************************************
126      do i=1,ND
127      RDCP=(RD*(1.-RR(I))+RR(I)*RV)/ (CPD*(1.-RR(I))+RR(I)*CPV)
128 
129      TLS(i)=T(I)*(1000.0/P(I))**RDCP
130      enddo
131 
132 
133 
134 
135************************************************************
136 
137 
138      PRECIP=0.0
139      IFLAG=1
140C
141C   ***                    SPECIFY PARAMETERS                        ***
142C   ***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE   ***
143C   ***       PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO         ***
144C   ***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP.      ***
145C   ***            EFFICIENCY IS ASSUMED TO BE UNITY                 ***
146C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
147C   ***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE      ***
148C   ***                        OF CLOUD                              ***
149C   ***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF    ***
150C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
151C   ***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY)    ***
152C   ***           (BETA MUST BE LESS THAN OR EQUAL TO 1)             ***
153C   ***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE    ***
154C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
155C   ***                     IT MUST BE LESS THAN 0                   ***
156C
157      PBCRIT=150.0
158      PTCRIT=500.0
159      SIGD=0.01
160      SPFAC=0.15
161c sb:
162      EPMAX=0.993 ! precip efficiency less than unity
163c      EPMAX=1. ! precip efficiency less than unity
164C
165Cjyg
166CCC      BETA=0.96
167C  Beta is now expressed as a function of the characteristic time
168C  of the convective process.
169CCC        Old value : TAU = 15000.   !(for dtime = 600.s)
170CCC        Other value (inducing little change) :TAU = 8000.
171      TAU = 8000.
172      BETA = 1.-DTIME/TAU
173Cjyg
174CCC      ALPHA=1.0
175      ALPHA=1.5E-3*DTIME/TAU
176C        Increase alpha in order to compensate W decrease
177      ALPHA = ALPHA*1.5
178C
179Cjyg (voir CONVECT 3)
180CCC      DTCRIT=-0.2
181      DTCRIT=-2.
182Cgf&jyg
183CCC     DT pour l'overshoot.
184      DTOVSH = -0.2
185 
186C
187C           ***        INCREMENT THE COUNTER       ***
188C
189      SIG(ND)=SIG(ND)+1.0
190      SIG(ND)=AMIN1(SIG(ND),12.1)
191C
192C           ***    IF NOPT IS AN INTEGER OTHER THAN 0, CONVECT     ***
193C           ***     RETURNS ARRAYS T AND R THAT MAY HAVE BEEN      ***
194C           ***  ALTERED BY DRY ADIABATIC ADJUSTMENT; OTHERWISE    ***
195C           ***        THE RETURNED ARRAYS ARE UNALTERED.          ***
196C
197      NOPT=0
198C
199C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
200C
201C     ***  DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A  ***
202C     ***                BOUNDARY LAYER SCHEME !!!               ***
203C
204c test sb      DO 30 I=NL-1,1,-1
205c test sb         JN=0
206c test sb         DO 10 J=I+1,NL
207c test sb   10    IF(TH(J).LT.TH(I))JN=J
208c test sb         IF(JN.EQ.0)GOTO 30
209c test sb         AHM=0.0
210c test sb         RM=0.0
211c test sb         UM=0.0
212c test sb         VM=0.0
213c test sb         DO K=1,NTRA
214c test sb          TRATM(K)=0.0
215c test sb         END DO
216c test sb         DO 15 J=I,JN
217c test sb          AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1))
218c test sb          RM=RM+RR(J)*(PH(J)-PH(J+1))
219c test sb          UM=UM+U(J)*(PH(J)-PH(J+1))
220c test sb          VM=VM+V(J)*(PH(J)-PH(J+1))
221c test sb          DO K=1,NTRA
222c test sb           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
223c test sb          END DO
224c test sb   15    CONTINUE
225c test sb         DPHINV=1./(PH(I)-PH(JN+1))
226c test sb         RM=RM*DPHINV
227c test sb         UM=UM*DPHINV
228c test sb         VM=VM*DPHINV
229c test sb         DO K=1,NTRA
230c test sb          TRATM(K)=TRATM(K)*DPHINV
231c test sb         END DO
232c test sb         A2=0.0
233c test sb         DO 20 J=I,JN
234c test sb            RR(J)=RM
235c test sb          U(J)=UM
236c test sb          V(J)=VM
237c test sb          DO K=1,NTRA
238c test sb           TRA(J,K)=TRATM(K)
239c test sb          END DO
240c test sb            RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV)
241c test sb            X=(0.001*P(J))**RDCP
242c test sb            T(J)=X
243c test sb            A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1))
244c test sb   20    CONTINUE
245c test sb         DO 25 J=I,JN
246c test sb            TH(J)=AHM/A2
247c test sb            T(J)=T(J)*TH(J)
248c test sb   25    CONTINUE
249c test sb   30 CONTINUE
250C
251C   ***   RESET INPUT ARRAYS IF NOPT .NE. 0   ***
252C
253      IF (NOPT.NE.0)THEN
254         DO 35 I=1,ND
255            T1(I)=T(I)
256            R1(I)=RR(I)
257   35    CONTINUE
258      END IF
259C
260C  *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY
261C
262      GZ(1)=0.0
263      CPN(1)=CPD*(1.-RR(1))+RR(1)*CPV
264      H(1)=T(1)*CPN(1)
265      DO 40 I=2,NL
266        TVX=T(I)*(1.+RR(I)/EPS-RR(I))
267        TVY=T(I-1)*(1.+RR(I-1)/EPS-RR(I-1))
268        GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(P(I-1)-P(I))/PH(I)
269        CPN(I)=CPD*(1.-RR(I))+CPV*RR(I)
270        H(I)=T(I)*CPN(I)+GZ(I)
271   40 CONTINUE
272C
273C   ***  CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT LOWEST MODEL LEVEL ***
274C   ***       (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980)     ***
275C
276      IF (T(1).LT.250.0.OR.RR(1).LE.0.0)THEN
277         IFLAG=0
278c sb3d         print*,'je suis passe par 366'
279         RETURN
280      END IF
281
282cjyg1 Utilisation de la subroutine CLIFT
283CC      RH=RR(1)/RS(1)
284CC      CHI=T(1)/(1669.0-122.0*RH-T(1))
285CC      PLCL=P(1)*(RH**CHI)
286      CALL CLIFT(P(1),T(1),RR(1),RS(1),PLCL,DPLCLDT,DPLCLDR)
287cjyg2
288c sb3d      PRINT *,' em_plcl,p1,t1,r1,rs1,rh '
289c sb3d     $        ,PLCL,P(1),T(1),RR(1),RS(1),RH
290c
291      IF (PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN
292         IFLAG=2
293         RETURN
294      END IF
295Cjyg1
296C     Essais de modification de ICB
297C
298C   ***  CALCULATE FIRST LEVEL ABOVE LCL (=ICB)  ***
299C
300CC      ICB=NL-1
301CC      DO 50 I=2,NL-1
302CC         IF(P(I).LT.PLCL)THEN
303CC            ICB=MIN(ICB,I)   ! ICB sup ou egal a 2
304CC         END IF
305CC   50 CONTINUE
306CC      IF(ICB.EQ.(NL-1))THEN
307CC         IFLAG=3
308CC         RETURN
309CC      END IF
310C
311C   *** CALCULATE LAYER CONTAINING LCL (=ICB)   ***
312C
313      ICB=NL-1
314c sb      DO 50 I=2,NL-1
315      DO 50 I=3,NL-1 ! modif sb pour que ICB soit sup/egal a 2
316C   la modification consiste a comparer PLCL a PH et non a P:
317C   ICB est defini par :  PH(ICB)<PLCL<PH(ICB-!)
318         IF(PH(I).LT.PLCL)THEN
319            ICB=MIN(ICB,I)
320         END IF
321   50 CONTINUE
322      IF(ICB.EQ.(NL-1))THEN
323         IFLAG=3
324         RETURN
325      END IF
326      ICB = ICB-1 ! ICB sup ou egal a 2
327Cjyg2
328C
329C
330 
331C   *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL      ***
332C   ***  TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC             ***
333C   ***                   LIQUID WATER CONTENT                             ***
334C
335 
336cjyg1
337c make sure that "Cloud base" seen by TLIFT is actually the
338c fisrt level where adiabatic ascent is saturated
339       IF (PLCL .GT. P(ICB)) THEN
340c sb        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,TVP,TP,CLW,ND,NL)
341        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK,TVP,TP,CLW,ND,NL
342     :            ,DTVPDT1,DTVPDQ1)
343       ELSE
344c sb        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,TVP,TP,CLW,ND,NL)
345        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,NK,TVP,TP,CLW,ND,NL
346     :            ,DTVPDT1,DTVPDQ1)
347       ENDIF
348cjyg2
349 
350******************************************************************************
351****     SORTIE DE LA TEMPERATURE DE L ASCENDANCE NON DILUE
352******************************************************************************
353        do i=1,ND
354        TPS(i)=TP(i)
355        enddo
356 
357 
358******************************************************************************
359 
360C
361C   ***  SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF   ***
362C   ***          PRECIPITATION FALLING OUTSIDE OF CLOUD           ***
363C   ***      THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)     ***
364C
365      DO 55 I=1,NL
366         PDEN=PTCRIT-PBCRIT
367c
368cjyg
369ccc         EP(I)=(P(ICB)-P(I)-PBCRIT)/PDEN
370c sb         EP(I)=(PLCL-P(I)-PBCRIT)/PDEN
371         EP(I)=(PLCL-P(I)-PBCRIT)/PDEN * EPMAX ! sb
372c
373         EP(I)=AMAX1(EP(I),0.0)
374c sb         EP(I)=AMIN1(EP(I),1.0)
375         EP(I)=AMIN1(EP(I),EPMAX) ! sb
376         SIGP(I)=SPFAC
377C
378C   ***       CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL     ***
379C   ***                    VIRTUAL TEMPERATURE                    ***
380C
381         TV(I)=T(I)*(1.+RR(I)/EPS-RR(I))
382Ccd1
383C    . Keep all liquid water in lifted parcel (-> adiabatic CAPE)
384C
385ccc    TVP(I)=TVP(I)-TP(I)*(RR(1)-EP(I)*CLW(I))
386c!!! sb         TVP(I)=TVP(I)-TP(I)*RR(1) ! calcule dans tlift
387Ccd2
388C
389C   ***       Calculate first estimate of buoyancy
390C
391         BUOY(I) = TVP(I) - TV(I)
392   55 CONTINUE
393C
394C   ***   Set Cloud Base Buoyancy at (Plcl+DPbase) level buoyancy
395C
396      DPBASE = -40.   !That is 400m above LCL
397      PBASE = PLCL + DPBASE
398      TVPBASE = TVP(ICB  )*(PBASE -P(ICB+1))/(P(ICB)-P(ICB+1))
399     $         +TVP(ICB+1)*(P(ICB)-PBASE   )/(P(ICB)-P(ICB+1))
400      TVBASE = TV(ICB  )*(PBASE -P(ICB+1))/(P(ICB)-P(ICB+1))
401     $        +TV(ICB+1)*(P(ICB)-PBASE   )/(P(ICB)-P(ICB+1))
402C
403      BUOYBASE = TVPBASE - TVBASE
404C
405CC       Set buoyancy = BUOYBASE for all levels below BASE.
406CC       For safety, set : BUOY(ICB) = BUOYBASE
407      DO I = ICB,NL
408        IF (P(I) .GE. PBASE) THEN
409          BUOY(I) = BUOYBASE
410        ENDIF
411      ENDDO
412      BUOY(ICB) = BUOYBASE
413C
414c sb3d      print *,'buoybase,tvp_tv,tvpbase,tvbase,pbase,plcl'
415c sb3d     $,        buoybase,tvp(icb)-tv(icb),tvpbase,tvbase,pbase,plcl
416c sb3d      print *,'TVP ',(tvp(i),i=1,nl)
417c sb3d      print *,'TV ',(tv(i),i=1,nl)
418c sb3d      print *, 'P ',(p(i),i=1,nl)
419c sb3d      print *,'ICB ',icb
420C
421C   ***   MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE  ***
422C   ***    AND CLOUD BASE, AND THAT LIFTED AIR IS POSITIVELY BUOYANT  ***
423C   ***                         AT CLOUD BASE                         ***
424C   ***       IF NOT, RETURN TO CALLING PROGRAM AFTER RESETTING       ***
425C   ***                        SIG(I) AND W0(I)                       ***
426C
427Cjyg
428CCC      TDIF=TVP(ICB)-TV(ICB)
429      TDIF = BUOY(ICB)
430      ATH1=TH(1)
431Cjyg
432CCC      ATH=TH(ICB-1)-1.0
433      ATH=TH(ICB-1)-5.0
434c      ATH=0.                          ! ajout sb
435c      IF (ICB.GT.1) ATH=TH(ICB-1)-5.0 ! modif sb
436      IF(TDIF.LT.DTCRIT.OR.ATH.GT.ATH1)THEN
437         DO 60 I=1,NL
438            SIG(I)=BETA*SIG(I)-2.*ALPHA*TDIF*TDIF
439            SIG(I)=AMAX1(SIG(I),0.0)
440            W0(I)=BETA*W0(I)
441   60    CONTINUE
442         IFLAG=0
443         RETURN
444      END IF
445C
446 
447 
448C   ***  IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY ***
449C   ***        NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS       ***
450C
451      DO 70 I=1,NL
452         HP(I)=H(I)
453         WT(I)=0.001
454         RP(I)=RR(I)
455         UP(I)=U(I)
456         VP(I)=V(I)
457         DO 71 J=1,NTRA
458          TRAP(I,J)=TRA(I,J)
459   71    CONTINUE
460         NENT(I)=0
461         WATER(I)=0.0
462         EVAP(I)=0.0
463         B(I)=0.0
464         MP(I)=0.0
465         M(I)=0.0
466         LV(I)=ALV0-CPVMCL*(T(I)-273.15)
467         LVCP(I)=LV(I)/CPN(I)
468         DO 70 J=1,NL
469            QENT(I,J)=RR(J)
470            ELIJ(I,J)=0.0
471            MENT(I,J)=0.0
472            SIJ(I,J)=0.0
473          UENT(I,J)=U(J)
474          VENT(I,J)=V(J)
475          DO 70 K=1,NTRA
476           TRAENT(I,J,K)=TRA(J,K)
477   70 CONTINUE
478 
479      DELTI=1.0/DELT
480C
481C  ***  FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S       ***
482C  ***                LEVEL OF NEUTRAL BUOYANCY                   ***
483C
484      INB=NL-1
485      DO 80 I=ICB,NL-1
486Cjyg
487CCC         IF((TVP(I)-TV(I)).LT.DTCRIT)THEN
488         IF(BUOY(I).LT.DTOVSH)THEN
489            INB=MIN(INB,I)
490         END IF
491   80 CONTINUE
492 
493 
494 
495 
496C
497C   ***          RESET SIG(I) AND W0(I) FOR I>INB AND I<ICB       ***
498C
499      IF(INB.LT.(NL-1))THEN
500         DO 85 I=INB+1,NL-1
501Cjyg
502CCC            SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(INB)-TVP(INB))*
503CCC     1              ABS(TV(INB)-TVP(INB))
504            SIG(I)=BETA*SIG(I)+2.*ALPHA*BUOY(INB)*
505     1              ABS(BUOY(INB))
506            SIG(I)=AMAX1(SIG(I),0.0)
507            W0(I)=BETA*W0(I)
508   85    CONTINUE
509      END IF
510      DO 87 I=1,ICB
511Cjyg
512CCC         SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(ICB)-TVP(ICB))*
513CCC     1           (TV(ICB)-TVP(ICB))
514         SIG(I)=BETA*SIG(I)-2.*ALPHA*BUOY(ICB)*BUOY(ICB)
515         SIG(I)=AMAX1(SIG(I),0.0)
516         W0(I)=BETA*W0(I)
517   87 CONTINUE
518C
519C   ***    RESET FRACTIONAL AREAS OF UPDRAFTS AND W0 AT INITIAL TIME    ***
520C   ***           AND AFTER 10 TIME STEPS OF NO CONVECTION              ***
521C
522 
523      IF(SIG(ND).LT.1.5.OR.SIG(ND).GT.12.0)THEN
524         DO 90 I=1,NL-1
525            SIG(I)=0.0
526            W0(I)=0.0
527   90    CONTINUE
528      END IF
529C
530C   ***   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL   ***
531C
532      DO 95 I=ICB,INB
533         HP(I)=H(1)+(LV(I)+(CPD-CPV)*T(I))*EP(I)*CLW(I)
534   95 CONTINUE
535C
536C   ***  CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE),  ***
537C   ***     VERTICAL VELOCITY (W), FRACTIONAL AREA COVERED BY     ***
538C   ***     UNDILUTE UPDRAFT (SIG),  AND UPDRAFT MASS FLUX (M)    ***
539C
540      CAPE=0.0
541C
542      DO 98 I=ICB+1,INB
543Cjyg1
544CCC         CAPE=CAPE+RD*(TVP(I-1)-TV(I-1))*(PH(I-1)-PH(I))/P(I-1)
545CCC         DCAPE=RD*BUOY(I-1)*(PH(I-1)-PH(I))/P(I-1)
546CCC         DLNP=(PH(I-1)-PH(I))/P(I-1)
547C          The interval on which CAPE is computed starts at PBASE :
548         DELTAP = MIN(PBASE,PH(I-1))-MIN(PBASE,PH(I))
549         CAPE=CAPE+RD*BUOY(I-1)*DELTAP/P(I-1)
550         DCAPE=RD*BUOY(I-1)*DELTAP/P(I-1)
551         DLNP=DELTAP/P(I-1)
552Cjyg2
553c sb3d         print *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape
554         CAPE=AMAX1(0.0,CAPE)
555C
556         SIGOLD=SIG(I)
557         DTMIN=100.0
558         DO 97 J=ICB,I-1
559Cjyg
560CCC            DTMIN=AMIN1(DTMIN,(TVP(J)-TV(J)))
561            DTMIN=AMIN1(DTMIN,BUOY(J))
562   97    CONTINUE
563c sb3d     print *, 'DTMIN, BETA, ALPHA, SIG = ',DTMIN,BETA,ALPHA,SIG(I)
564         SIG(I)=BETA*SIG(I)+ALPHA*DTMIN*ABS(DTMIN)
565         SIG(I)=AMAX1(SIG(I),0.0)
566         SIG(I)=AMIN1(SIG(I),0.01)
567         FAC=AMIN1(((DTCRIT-DTMIN)/DTCRIT),1.0)
568Cjyg
569CC    Essais de reduction de la vitesse
570CC         FAC = FAC*.5
571C
572         W=(1.-BETA)*FAC*SQRT(CAPE)+BETA*W0(I)
573         AMU=0.5*(SIG(I)+SIGOLD)*W
574         M(I)=AMU*0.007*P(I)*(PH(I)-PH(I+1))/TV(I)
575         W0(I)=W
576   98 CONTINUE
577      W0(ICB)=0.5*W0(ICB+1)
578      M(ICB)=0.5*M(ICB+1)*(PH(ICB)-PH(ICB+1))/(PH(ICB+1)-PH(ICB+2))
579      SIG(ICB)=SIG(ICB+1)
580      SIG(ICB-1)=SIG(ICB)
581cjyg1
582c sb3d      print *, 'Cloud base, c. top, CAPE',ICB,INB,cape
583c sb3d      print *, 'SIG ',(sig(i),i=1,inb)
584c sb3d      print *, 'W ',(w0(i),i=1,inb)
585c sb3d      print *, 'M ',(m(i), i=1,inb)
586c sb3d      print *, 'Dt1 ',(tvp(i)-tv(i),i=1,inb)
587c sb3d      print *, 'Dt_vrai ',(buoy(i),i=1,inb)
588Cjyg2
589C
590C   ***  CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING  ***
591C   ***     RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING     ***
592C   ***                        FRACTION (SIJ)                          ***
593C
594 
595 
596      DO 170 I=ICB,INB
597         RTI=RR(1)-EP(I)*CLW(I)
598         DO 160 J=ICB-1,INB
599            BF2=1.+LV(J)*LV(J)*RS(J)/(RV*T(J)*T(J)*CPD)
600            ANUM=H(J)-HP(I)+(CPV-CPD)*T(J)*(RTI-RR(J))
601            DENOM=H(I)-HP(I)+(CPD-CPV)*(RR(I)-RTI)*T(J)
602            DEI=DENOM
603            IF(ABS(DEI).LT.0.01)DEI=0.01
604            SIJ(I,J)=ANUM/DEI
605            SIJ(I,I)=1.0
606            ALTEM=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI-RS(J)
607            ALTEM=ALTEM/BF2
608            CWAT=CLW(J)*(1.-EP(J))
609            STEMP=SIJ(I,J)
610            IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR.
611     1      ALTEM.GT.CWAT).AND.J.GT.I)THEN
612            ANUM=ANUM-LV(J)*(RTI-RS(J)-CWAT*BF2)
613            DENOM=DENOM+LV(J)*(RR(I)-RTI)
614            IF(ABS(DENOM).LT.0.01)DENOM=0.01
615            SIJ(I,J)=ANUM/DENOM
616            ALTEM=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI-RS(J)
617            ALTEM=ALTEM-(BF2-1.)*CWAT
618            END IF
619 
620 
621            IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.95)THEN
622               QENT(I,J)=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI
623               UENT(I,J)=SIJ(I,J)*U(I)+(1.-SIJ(I,J))*U(NK)
624               VENT(I,J)=SIJ(I,J)*V(I)+(1.-SIJ(I,J))*V(NK)
625               DO K=1,NTRA
626               TRAENT(I,J,K)=SIJ(I,J)*TRA(I,K)+(1.-SIJ(I,J))*
627     1         TRA(NK,K)
628               END DO
629               ELIJ(I,J)=ALTEM
630               ELIJ(I,J)=AMAX1(0.0,ELIJ(I,J))
631               MENT(I,J)=M(I)/(1.-SIJ(I,J))
632               NENT(I)=NENT(I)+1
633            END IF
634            SIJ(I,J)=AMAX1(0.0,SIJ(I,J))
635            SIJ(I,J)=AMIN1(1.0,SIJ(I,J))
636  160    CONTINUE
637C
638C   ***   IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS  ***
639C   ***   AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES  ***
640C
641         IF(NENT(I).EQ.0)THEN
642            MENT(I,I)=M(I)
643            QENT(I,I)=RR(NK)-EP(I)*CLW(I)
644           UENT(I,I)=U(NK)
645           VENT(I,I)=V(NK)
646           DO J=1,NTRA
647            TRAENT(I,I,J)=TRA(NK,J)
648           END DO
649            ELIJ(I,I)=CLW(I)
650            SIJ(I,I)=1.0
651         END IF
652C
653         DO J = ICB-1,INB
654           SIGIJ(I,J) = SIJ(I,J)
655         ENDDO
656C       
657  170 CONTINUE
658C
659C   ***  NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL  ***
660C   ***              PROBABILITIES OF MIXING                     ***
661C
662 
663      DO 200 I=ICB,INB
664      IF(NENT(I).NE.0)THEN
665       QP=RR(1)-EP(I)*CLW(I)
666       ANUM=H(I)-HP(I)-LV(I)*(QP-RS(I))+(CPV-CPD)*T(I)*
667     1    (QP-RR(I))
668       DENOM=H(I)-HP(I)+LV(I)*(RR(I)-QP)+
669     1    (CPD-CPV)*T(I)*(RR(I)-QP)
670       IF(ABS(DENOM).LT.0.01)DENOM=0.01
671       SCRIT=ANUM/DENOM
672       ALT=QP-RS(I)+SCRIT*(RR(I)-QP)
673       IF(SCRIT.LE.0.0.OR.ALT.LE.0.0)SCRIT=1.0
674       SMAX=0.0
675       ASIJ=0.0
676        DO 175 J=INB,ICB-1,-1
677        IF(SIJ(I,J).GT.1.0E-16.AND.SIJ(I,J).LT.0.95)THEN
678         WGH=1.0
679         IF(J.GT.I)THEN
680          SJMAX=AMAX1(SIJ(I,J+1),SMAX)
681          SJMAX=AMIN1(SJMAX,SCRIT)
682          SMAX=AMAX1(SIJ(I,J),SMAX)
683          SJMIN=AMAX1(SIJ(I,J-1),SMAX)
684          SJMIN=AMIN1(SJMIN,SCRIT)
685          IF(SIJ(I,J).LT.(SMAX-1.0E-16))WGH=0.0
686          SMID=AMIN1(SIJ(I,J),SCRIT)
687         ELSE
688          SJMAX=AMAX1(SIJ(I,J+1),SCRIT)
689          SMID=AMAX1(SIJ(I,J),SCRIT)
690          SJMIN=0.0
691          IF(J.GT.1)SJMIN=SIJ(I,J-1)
692          SJMIN=AMAX1(SJMIN,SCRIT)
693         END IF
694         DELP=ABS(SJMAX-SMID)
695         DELM=ABS(SJMIN-SMID)
696         ASIJ=ASIJ+WGH*(DELP+DELM)
697         MENT(I,J)=MENT(I,J)*(DELP+DELM)*WGH
698        END IF
699  175       CONTINUE
700       ASIJ=AMAX1(1.0E-16,ASIJ)
701       ASIJ=1.0/ASIJ
702       DO 180 J=ICB-1,INB
703        MENT(I,J)=MENT(I,J)*ASIJ
704  180    CONTINUE
705       ASUM=0.0
706       BSUM=0.0
707       DO 190 J=ICB-1,INB
708        ASUM=ASUM+MENT(I,J)
709        MENT(I,J)=MENT(I,J)*SIG(J)
710        BSUM=BSUM+MENT(I,J)
711  190       CONTINUE
712       BSUM=AMAX1(BSUM,1.0E-16)
713       BSUM=1.0/BSUM
714       DO 195 J=ICB-1,INB
715        MENT(I,J)=MENT(I,J)*ASUM*BSUM   
716  195       CONTINUE
717       CSUM=0.0
718       DO 197 J=ICB-1,INB
719        CSUM=CSUM+MENT(I,J)
720  197       CONTINUE
721 
722       IF(CSUM.LT.M(I))THEN
723        NENT(I)=0
724        MENT(I,I)=M(I)
725        QENT(I,I)=RR(1)-EP(I)*CLW(I)
726          UENT(I,I)=U(NK)
727          VENT(I,I)=V(NK)
728          DO J=1,NTRA
729           TRAENT(I,I,J)=TRA(NK,J)
730          END DO
731        ELIJ(I,I)=CLW(I)
732        SIJ(I,I)=1.0
733       END IF
734      END IF
735  200      CONTINUE
736 
737 
738 
739***************************************************************
740**       CALCUL DES MENTS(I,J) ET DES QENTS(I,J)
741**************************************************************
742 
743         DO im=1,nd
744         do jm=1,nd
745 
746         QENTS(im,jm)=QENT(im,jm)
747         MENTS(im,jm)=MENT(im,jm)
748         enddo
749         enddo
750 
751***********************************************************
752 
753 
754 
755C
756C   ***  CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING    ***
757C   ***             DOWNDRAFT CALCULATION                      ***
758C
759        IF(EP(INB).LT.0.0001)GOTO 405
760C
761C   ***  INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER   ***
762C   ***                AND CONDENSED WATER FLUX                    ***
763C
764        WFLUX=0.0
765        TINV=1./3.
766C
767C    ***                    BEGIN DOWNDRAFT LOOP                    ***
768C
769        DO 400 I=INB,1,-1
770C
771C    ***              CALCULATE DETRAINED PRECIPITATION             ***
772C
773 
774 
775        WDTRAIN=10.0*EP(I)*M(I)*CLW(I)
776        IF(I.GT.1)THEN
777         DO 320 J=1,I-1
778       AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I)
779       AWAT=AMAX1(AWAT,0.0)
780  320    WDTRAIN=WDTRAIN+10.0*AWAT*MENT(J,I)
781        END IF
782C
783C    ***    FIND RAIN WATER AND EVAPORATION USING PROVISIONAL   ***
784C    ***              ESTIMATES OF RP(I)AND RP(I-1)             ***
785C
786 
787 
788        WT(I)=45.0
789      IF(I.LT.INB)THEN
790       RP(I)=RP(I+1)+(CPD*(T(I+1)-T(I))+GZ(I+1)-GZ(I))/LV(I)
791       RP(I)=0.5*(RP(I)+RR(I))
792      END IF
793      RP(I)=AMAX1(RP(I),0.0)
794      RP(I)=AMIN1(RP(I),RS(I))
795      RP(INB)=RR(INB)
796      IF(I.EQ.1)THEN
797       AFAC=P(1)*(RS(1)-RP(1))/(1.0E4+2000.0*P(1)*RS(1))
798      ELSE
799       RP(I-1)=RP(I)+(CPD*(T(I)-T(I-1))+GZ(I)-GZ(I-1))/LV(I)
800       RP(I-1)=0.5*(RP(I-1)+RR(I-1))
801       RP(I-1)=AMIN1(RP(I-1),RS(I-1))
802       RP(I-1)=AMAX1(RP(I-1),0.0)
803       AFAC1=P(I)*(RS(I)-RP(I))/(1.0E4+2000.0*P(I)*RS(I))
804       AFAC2=P(I-1)*(RS(I-1)-RP(I-1))/(1.0E4+
805     1    2000.0*P(I-1)*RS(I-1))
806       AFAC=0.5*(AFAC1+AFAC2)
807      END IF
808      IF(I.EQ.INB)AFAC=0.0
809        AFAC=AMAX1(AFAC,0.0)
810        BFAC=1./(SIGD*WT(I))
811C
812Cjyg1
813CCC        SIGT=1.0
814CCC        IF(I.GE.ICB)SIGT=SIGP(I)
815C Prise en compte de la variation progressive de SIGT dans
816C les couches ICB et ICB-1:
817C       Pour PLCL<PH(I+1), PR1=0 & PR2=1
818C       Pour PLCL>PH(I),   PR1=1 & PR2=0
819C       Pour PH(I+1)<PLCL<PH(I), PR1 est la proportion a cheval
820C    sur le nuage, et PR2 est la proportion sous la base du
821C    nuage.
822         PR1 =(PLCL-PH(I+1))/(PH(I)-PH(I+1))
823         PR1 = MAX(0.,MIN(1.,PR1))
824         PR2 = (PH(I)-PLCL)/(PH(I)-PH(I+1))
825         PR2 = MAX(0.,MIN(1.,PR2))
826         SIGT = SIGP(I)*PR1 + PR2
827c sb3d         print *,'i,sigt,pr1,pr2', i,sigt,pr1,pr2
828Cjyg2
829C
830        B6=BFAC*50.*SIGD*(PH(I)-PH(I+1))*SIGT*AFAC
831        C6=WATER(I+1)+BFAC*WDTRAIN-50.*SIGD*BFAC*
832     1   (PH(I)-PH(I+1))*EVAP(I+1)
833      IF(C6.GT.0.0)THEN
834         REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6))
835         EVAP(I)=SIGT*AFAC*REVAP
836         WATER(I)=REVAP*REVAP
837      ELSE
838       EVAP(I)=-EVAP(I+1)+0.02*(WDTRAIN+SIGD*WT(I)*
839     1    WATER(I+1))/(SIGD*(PH(I)-PH(I+1)))
840      END IF
841 
842 
843C
844C    ***  CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER     ***
845C    ***              HYDROSTATIC APPROXIMATION                 ***
846C
847        IF(I.EQ.1)GOTO 360
848      TEVAP=AMAX1(0.0,EVAP(I))
849      DELTH=AMAX1(0.001,(TH(I)-TH(I-1)))
850      MP(I)=10.*LVCP(I)*SIGD*TEVAP*(P(I-1)-P(I))/DELTH
851C
852C    ***           IF HYDROSTATIC ASSUMPTION FAILS,             ***
853C    ***   SOLVE CUBIC DIFFERENCE EQUATION FOR DOWNDRAFT THETA  ***
854C    ***  AND MASS FLUX FROM TWO SIMULTANEOUS DIFFERENTIAL EQNS ***
855C
856      AMFAC=SIGD*SIGD*70.0*PH(I)*(P(I-1)-P(I))*
857     1   (TH(I)-TH(I-1))/(TV(I)*TH(I))
858      AMP2=ABS(MP(I+1)*MP(I+1)-MP(I)*MP(I))
859      IF(AMP2.GT.(0.1*AMFAC))THEN
860         XF=100.0*SIGD*SIGD*SIGD*(PH(I)-PH(I+1))
861         TF=B(I)-5.0*(TH(I)-TH(I-1))*T(I)/(LVCP(I)*SIGD*TH(I))
862         AF=XF*TF+MP(I+1)*MP(I+1)*TINV
863         BF=2.*(TINV*MP(I+1))**3+TINV*MP(I+1)*XF*TF+50.*
864     1    (P(I-1)-P(I))*XF*TEVAP
865         FAC2=1.0
866         IF(BF.LT.0.0)FAC2=-1.0
867         BF=ABS(BF)
868         UR=0.25*BF*BF-AF*AF*AF*TINV*TINV*TINV
869         IF(UR.GE.0.0)THEN
870          SRU=SQRT(UR)
871          FAC=1.0
872          IF((0.5*BF-SRU).LT.0.0)FAC=-1.0
873          MP(I)=MP(I+1)*TINV+(0.5*BF+SRU)**TINV+
874     1     FAC*(ABS(0.5*BF-SRU))**TINV
875         ELSE
876          D=ATAN(2.*SQRT(-UR)/(BF+1.0E-28))
877          IF(FAC2.LT.0.0)D=3.14159-D
878          MP(I)=MP(I+1)*TINV+2.*SQRT(AF*TINV)*COS(D*TINV)
879         END IF
880         MP(I)=AMAX1(0.0,MP(I))
881         B(I-1)=B(I)+100.0*(P(I-1)-P(I))*TEVAP/(MP(I)+SIGD*0.1)-
882     1    10.0*(TH(I)-TH(I-1))*T(I)/(LVCP(I)*SIGD*TH(I))
883         B(I-1)=AMAX1(B(I-1),0.0)
884      END IF
885 
886 
887C
888C   ***         LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION      ***
889C
890      AMPMAX=2.0*(PH(I)-PH(I+1))*DELTI
891      AMP2=2.0*(PH(I-1)-PH(I))*DELTI
892      AMPMAX=AMIN1(AMPMAX,AMP2)
893      MP(I)=AMIN1(MP(I),AMPMAX)
894C
895C    ***      FORCE MP TO DECREASE LINEARLY TO ZERO                 ***
896C    ***       BETWEEN CLOUD BASE AND THE SURFACE                   ***
897C
898          IF(P(I).GT.P(ICB))THEN
899           MP(I)=MP(ICB)*(P(1)-P(I))/(P(1)-P(ICB))
900          END IF
901  360   CONTINUE
902C
903C    ***       FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT     ***
904C
905        IF(I.EQ.INB)GOTO 400
906      RP(I)=RR(I)
907        IF(MP(I).GT.MP(I+1))THEN
908        RP(I)=RP(I+1)*MP(I+1)+RR(I)*(MP(I)-MP(I+1))+
909     1       5.*SIGD*(PH(I)-PH(I+1))*(EVAP(I+1)+EVAP(I))
910        RP(I)=RP(I)/MP(I)
911          UP(I)=UP(I+1)*MP(I+1)+U(I)*(MP(I)-MP(I+1))
912         UP(I)=UP(I)/MP(I)
913          VP(I)=VP(I+1)*MP(I+1)+V(I)*(MP(I)-MP(I+1))
914         VP(I)=VP(I)/MP(I)
915          DO J=1,NTRA
916           TRAP(I,J)=TRAP(I+1,J)*MP(I+1)+
917     s     TRAP(I,J)*(MP(I)-MP(I+1))
918           TRAP(I,J)=TRAP(I,J)/MP(I)
919          END DO
920        ELSE
921        IF(MP(I+1).GT.1.0E-16)THEN
922           RP(I)=RP(I+1)+5.0*SIGD*(PH(I)-PH(I+1))*(EVAP(I+1)+
923     1      EVAP(I))/MP(I+1)
924            UP(I)=UP(I+1)
925            VP(I)=VP(I+1)
926            DO J=1,NTRA
927             TRAP(I,J)=TRAP(I+1,J)
928            END DO
929        END IF
930        END IF
931      RP(I)=AMIN1(RP(I),RS(I))
932      RP(I)=AMAX1(RP(I),0.0)
933  400   CONTINUE
934C
935C   ***  CALCULATE SURFACE PRECIPITATION IN MM/DAY     ***
936C
937        PRECIP=WT(1)*SIGD*WATER(1)*8640.0
938  405   CONTINUE
939C
940C   ***  CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE  ***
941C   ***                      AND MIXING RATIO                        ***
942C
943      DPINV=1.0/(PH(1)-PH(2))
944        AM=0.0
945        DO 410 K=2,INB
946  410   AM=AM+M(K)
947      IF((0.1*DPINV*AM).GE.DELTI)IFLAG=4
948      FT(1)=0.1*DPINV*AM*(T(2)-T(1)+(GZ(2)-GZ(1))/CPN(1))
949        FT(1)=FT(1)-0.5*LVCP(1)*SIGD*(EVAP(1)+EVAP(2))
950        FT(1)=FT(1)-0.09*SIGD*MP(2)*T(1)*B(1)*DPINV
951      FT(1)=FT(1)+0.01*SIGD*WT(1)*(CL-CPD)*WATER(2)*(T(2)-
952     1   T(1))*DPINV/CPN(1)
953        FR(1)=0.1*MP(2)*(RP(2)-RR(1))*
954     1    DPINV+SIGD*0.5*(EVAP(1)+EVAP(2))
955        FR(1)=FR(1)+0.1*AM*(RR(2)-RR(1))*DPINV
956        FU(1)=FU(1)+0.1*DPINV*(MP(2)*(UP(2)-U(1))+AM*(U(2)-U(1)))
957        FV(1)=FV(1)+0.1*DPINV*(MP(2)*(VP(2)-V(1))+AM*(V(2)-V(1)))
958        DO J=1,NTRA
959         FTRA(1,J)=FTRA(1,J)+0.1*DPINV*(MP(2)*(TRAP(2,J)-TRA(1,J))+
960     1    AM*(TRA(2,J)-TRA(1,J)))
961        END DO
962        AMDE=0.0
963        DO 415 J=2,INB
964         FR(1)=FR(1)+0.1*DPINV*MENT(J,1)*(QENT(J,1)-RR(1))
965         FU(1)=FU(1)+0.1*DPINV*MENT(J,1)*(UENT(J,1)-U(1))
966         FV(1)=FV(1)+0.1*DPINV*MENT(J,1)*(VENT(J,1)-V(1))
967         DO K=1,NTRA
968          FTRA(1,K)=FTRA(1,K)+0.1*DPINV*MENT(J,1)*(TRAENT(J,1,K)-
969     1     TRA(1,K))
970         END DO
971  415      CONTINUE
972C
973C   ***  CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO  ***
974C   ***               AT LEVELS ABOVE THE LOWEST LEVEL                   ***
975C
976C   ***  FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES  ***
977C   ***                      THROUGH EACH LEVEL                          ***
978C
979 
980 
981        DO 500 I=2,INB
982        DPINV=1.0/(PH(I)-PH(I+1))
983      CPINV=1.0/CPN(I)
984        AMP1=0.0
985        DO 440 K=I+1,INB+1
986  440   AMP1=AMP1+M(K)
987        DO 450 K=1,I
988        DO 450 J=I+1,INB+1
989         AMP1=AMP1+MENT(K,J)
990  450   CONTINUE
991      IF((0.1*DPINV*AMP1).GE.DELTI)IFLAG=4
992        AD=0.0
993        DO 470 K=1,I-1
994        DO 470 J=I,INB
995  470   AD=AD+MENT(J,K)
996      FT(I)=0.1*DPINV*(AMP1*(T(I+1)-T(I)+(GZ(I+1)-GZ(I))*
997     1   CPINV)-AD*(T(I)-T(I-1)+(GZ(I)-GZ(I-1))*CPINV))
998     2   -0.5*SIGD*LVCP(I)*(EVAP(I)+EVAP(I+1))
999      RAT=CPN(I-1)*CPINV
1000        FT(I)=FT(I)-0.09*SIGD*(MP(I+1)*T(I)*
1001     1    B(I)-MP(I)*T(I-1)*RAT*B(I-1))*DPINV
1002      FT(I)=FT(I)+0.1*DPINV*MENT(I,I)*(HP(I)-H(I)+
1003     1    T(I)*(CPV-CPD)*(RR(I)-QENT(I,I)))*CPINV
1004      FT(I)=FT(I)+0.01*SIGD*WT(I)*(CL-CPD)*WATER(I+1)*
1005     1    (T(I+1)-T(I))*DPINV*CPINV
1006        FR(I)=0.1*DPINV*(AMP1*(RR(I+1)-RR(I))-
1007     1    AD*(RR(I)-RR(I-1)))
1008        FU(I)=FU(I)+0.1*DPINV*(AMP1*(U(I+1)-U(I))-
1009     1    AD*(U(I)-U(I-1)))
1010        FV(I)=FV(I)+0.1*DPINV*(AMP1*(V(I+1)-V(I))-
1011     1    AD*(V(I)-V(I-1)))
1012        DO K=1,NTRA
1013         FTRA(I,K)=FTRA(I,K)+0.1*DPINV*(AMP1*(TRA(I+1,K)-
1014     1    TRA(I,K))-AD*(TRA(I,K)-TRA(I-1,K)))
1015        END DO
1016        DO 480 K=1,I-1
1017       AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I)
1018       AWAT=AMAX1(AWAT,0.0)
1019         FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-AWAT
1020     1    -RR(I))
1021         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
1022         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
1023         DO J=1,NTRA
1024          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
1025     1     TRA(I,J))
1026         END DO
1027  480   CONTINUE
1028      DO 490 K=I,INB
1029       FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-RR(I))
1030         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
1031         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
1032         DO J=1,NTRA
1033          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
1034     1     TRA(I,J))
1035         END DO
1036  490      CONTINUE
1037        FR(I)=FR(I)+0.5*SIGD*(EVAP(I)+EVAP(I+1))+0.1*(MP(I+1)*
1038     1    (RP(I+1)-RR(I))-MP(I)*(RP(I)-RR(I-1)))*DPINV
1039        FU(I)=FU(I)+0.1*(MP(I+1)*(UP(I+1)-U(I))-MP(I)*
1040     1    (UP(I)-U(I-1)))*DPINV
1041        FV(I)=FV(I)+0.1*(MP(I+1)*(VP(I+1)-V(I))-MP(I)*
1042     1    (VP(I)-V(I-1)))*DPINV
1043        DO J=1,NTRA
1044         FTRA(I,J)=FTRA(I,J)+0.1*DPINV*(MP(I+1)*(TRAP(I+1,J)-TRA(I,J))-
1045     1    MP(I)*(TRAP(I,J)-TRAP(I-1,J)))
1046        END DO
1047  500   CONTINUE
1048 
1049 
1050 
1051C
1052C   ***   MOVE THE DETRAINMENT AT LEVEL INB DOWN TO LEVEL INB-1   ***
1053C   ***        IN SUCH A WAY AS TO PRESERVE THE VERTICALLY        ***
1054C   ***          INTEGRATED ENTHALPY AND WATER TENDENCIES         ***
1055C
1056      AX=0.1*MENT(INB,INB)*(HP(INB)-H(INB)+T(INB)*
1057     1    (CPV-CPD)*(RR(INB)-QENT(INB,INB)))/(CPN(INB)*
1058     2    (PH(INB)-PH(INB+1)))
1059      FT(INB)=FT(INB)-AX
1060      FT(INB-1)=FT(INB-1)+AX*CPN(INB)*(PH(INB)-PH(INB+1))/
1061     1    (CPN(INB-1)*(PH(INB-1)-PH(INB)))
1062      BX=0.1*MENT(INB,INB)*(QENT(INB,INB)-RR(INB))/
1063     1    (PH(INB)-PH(INB+1))
1064      FR(INB)=FR(INB)-BX
1065      FR(INB-1)=FR(INB-1)+BX*(PH(INB)-PH(INB+1))/
1066     1    (PH(INB-1)-PH(INB))
1067      CX=0.1*MENT(INB,INB)*(UENT(INB,INB)-U(INB))/
1068     1    (PH(INB)-PH(INB+1))
1069      FU(INB)=FU(INB)-CX
1070      FU(INB-1)=FU(INB-1)+CX*(PH(INB)-PH(INB+1))/
1071     1    (PH(INB-1)-PH(INB))
1072      DX=0.1*MENT(INB,INB)*(VENT(INB,INB)-V(INB))/
1073     1    (PH(INB)-PH(INB+1))
1074      FV(INB)=FV(INB)-DX
1075      FV(INB-1)=FV(INB-1)+DX*(PH(INB)-PH(INB+1))/
1076     1    (PH(INB-1)-PH(INB))
1077      DO J=1,NTRA
1078      EX=0.1*MENT(INB,INB)*(TRAENT(INB,INB,J)
1079     1    -TRA(INB,J))/(PH(INB)-PH(INB+1))
1080      FTRA(INB,J)=FTRA(INB,J)-EX
1081      FTRA(INB-1,J)=FTRA(INB-1,J)+EX*
1082     1     (PH(INB)-PH(INB+1))/(PH(INB-1)-PH(INB))
1083      ENDDO   
1084C
1085C   ***    HOMOGINIZE TENDENCIES BELOW CLOUD BASE    ***
1086C
1087      ASUM=0.0
1088      BSUM=0.0
1089      CSUM=0.0
1090        DSUM=0.0
1091      DO 650 I=1,ICB-1
1092       ASUM=ASUM+FT(I)*(PH(I)-PH(I+1))
1093         BSUM=BSUM+FR(I)*(LV(I)+(CL-CPD)*(T(I)-T(1)))*
1094     1    (PH(I)-PH(I+1))
1095       CSUM=CSUM+(LV(I)+(CL-CPD)*(T(I)-T(1)))*(PH(I)-PH(I+1))
1096       DSUM=DSUM+T(I)*(PH(I)-PH(I+1))/TH(I)
1097  650      CONTINUE
1098      DO 700 I=1,ICB-1
1099       FT(I)=ASUM*T(I)/(TH(I)*DSUM)
1100       FR(I)=BSUM/CSUM
1101  700      CONTINUE
1102C
1103C   ***           RESET COUNTER AND RETURN           ***
1104C
1105      SIG(ND)=2.0
1106c
1107c
1108      do i = 1, nd
1109         upwd(i) = 0.0
1110         dnwd(i) = 0.0
1111c sb       dnwd0(i) = - mp(i)
1112      enddo
1113c
1114      do i = 1, nl
1115       dnwd0(i) = - mp(i)
1116      enddo
1117      do i = nl+1, nd
1118       dnwd0(i) = 0.
1119      enddo
1120c
1121      do i = icb, inb
1122         upwd(i) = 0.0
1123         dnwd(i) = 0.0
1124
1125         do k =i, inb
1126            up1=0.0
1127            dn1=0.0
1128            do n = 1, i-1
1129               up1 = up1 + ment(n,k)
1130               dn1 = dn1 - ment(k,n)
1131            enddo
1132            upwd(i) = upwd(i)+ m(k) + up1
1133            dnwd(i) = dnwd(i) + dn1
1134         enddo
1135        enddo
1136 
1137cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1138c        DETERMINATION DE LA VARIATION DE FLUX ASCENDANT ENTRE
1139C        DEUX NIVEAU NON DILUE Mike
1140cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1141 
1142 
1143c sb      do i=1,ND
1144c sb      Mike(i)=M(i)
1145c sb      enddo
1146 
1147      do i = 1, NL
1148       Mike(i) = M(i)
1149      enddo
1150      do i = NL+1, ND
1151       Mike(i) = 0.
1152      enddo
1153 
1154      do i=1,nd
1155      Ma(i)=0
1156      enddo
1157 
1158c sb      do i=1,nd
1159c sb      do j=i,nd
1160c sb      Ma(i)=Ma(i)+M(j)
1161c sb      enddo
1162c sb      enddo
1163
1164      do i = 1, NL
1165      do j = i, NL
1166       Ma(i) = Ma(i) + M(j)
1167      enddo
1168      enddo
1169c
1170      do i = NL+1, ND
1171       Ma(i) = 0.
1172      enddo
1173c
1174      do i=1,ICB-1
1175      Ma(i)=0
1176      enddo
1177 
1178 
1179 
1180ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1181C        ICB REPRESENTE DE NIVEAU OU SE TROUVE LA
1182c        BASE DU NUAGE , ET INB LE TOP DU NUAGE
1183cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1184 
1185 
1186       do i=1,ND
1187       Mke(i)=upwd(i)+dnwd(i)
1188       enddo
1189c$$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1190c$$$         call writeg1d(1,klev,ma,'ma  ','ma  ')
1191c$$$          call writeg1d(1,klev,upwd,'upwd  ','upwd  ')
1192c$$$          call writeg1d(1,klev,dnwd,'dnwd  ','dnwd  ')
1193c$$$          call writeg1d(1,klev,dnwd0,'dnwd0  ','dnwd0  ')
1194c$$$          call writeg1d(1,klev,tvp,'tvp  ','tvp  ')
1195c$$$          call writeg1d(1,klev,tra(1:klev,3),'tra3  ','tra3  ')
1196c$$$          call writeg1d(1,klev,tra(1:klev,4),'tra4  ','tra4  ')
1197c$$$          call writeg1d(1,klev,tra(1:klev,5),'tra5  ','tra5  ')
1198c$$$          call writeg1d(1,klev,tra(1:klev,6),'tra6  ','tra6  ')
1199c$$$          call writeg1d(1,klev,tra(1:klev,7),'tra7  ','tra7  ')
1200c$$$          call writeg1d(1,klev,tra(1:klev,8),'tra8  ','tra8  ')
1201c$$$          call writeg1d(1,klev,tra(1:klev,9),'tra9  ','tra9  ')
1202c$$$          call writeg1d(1,klev,tra(1:klev,10),'tra10','tra10')
1203c$$$          call writeg1d(1,klev,tra(1:klev,11),'tra11','tra11')
1204c$$$          call writeg1d(1,klev,tra(1:klev,12),'tra12','tra12')
1205c$$$          call writeg1d(1,klev,tra(1:klev,13),'tra13','tra13')
1206c$$$          call writeg1d(1,klev,tra(1:klev,14),'tra14','tra14')
1207c$$$          call writeg1d(1,klev,tra(1:klev,15),'tra15','tra15')
1208c$$$          call writeg1d(1,klev,tra(1:klev,16),'tra16','tra16')
1209c$$$          call writeg1d(1,klev,tra(1:klev,17),'tra17','tra17')
1210c$$$          call writeg1d(1,klev,tra(1:klev,18),'tra18','tra18')
1211c$$$          call writeg1d(1,klev,tra(1:klev,19),'tra19','tra19')
1212c$$$          call writeg1d(1,klev,tra(1:klev,20),'tra20','tra20 ')
1213c$$$          call writeg1d(1,klev,trap(1:klev,1),'trp1','trp1')
1214c$$$          call writeg1d(1,klev,trap(1:klev,2),'trp2','trp2')
1215c$$$          call writeg1d(1,klev,trap(1:klev,3),'trp3','trp3')
1216c$$$          call writeg1d(1,klev,trap(1:klev,4),'trp4','trp4')
1217c$$$          call writeg1d(1,klev,trap(1:klev,5),'trp5','trp5')
1218c$$$          call writeg1d(1,klev,trap(1:klev,10),'trp10','trp10')
1219c$$$          call writeg1d(1,klev,trap(1:klev,12),'trp12','trp12')
1220c$$$          call writeg1d(1,klev,trap(1:klev,15),'trp15','trp15')
1221c$$$          call writeg1d(1,klev,trap(1:klev,20),'trp20','trp20')
1222c$$$          call writeg1d(1,klev,ftra(1:klev,1),'ftr1  ','ftr1  ')
1223c$$$          call writeg1d(1,klev,ftra(1:klev,2),'ftr2  ','ftr2  ')
1224c$$$          call writeg1d(1,klev,ftra(1:klev,3),'ftr3  ','ftr3  ')
1225c$$$          call writeg1d(1,klev,ftra(1:klev,4),'ftr4  ','ftr4  ')
1226c$$$          call writeg1d(1,klev,ftra(1:klev,5),'ftr5  ','ftr5  ')
1227c$$$          call writeg1d(1,klev,ftra(1:klev,6),'ftr6  ','ftr6  ')
1228c$$$          call writeg1d(1,klev,ftra(1:klev,7),'ftr7  ','ftr7  ')
1229c$$$          call writeg1d(1,klev,ftra(1:klev,8),'ftr8  ','ftr8  ')
1230c$$$          call writeg1d(1,klev,ftra(1:klev,9),'ftr9  ','ftr9  ')
1231c$$$          call writeg1d(1,klev,ftra(1:klev,10),'ftr10','ftr10')
1232c$$$          call writeg1d(1,klev,ftra(1:klev,11),'ftr11','ftr11')
1233c$$$          call writeg1d(1,klev,ftra(1:klev,12),'ftr12','ftr12')
1234c$$$          call writeg1d(1,klev,ftra(1:klev,13),'ftr13','ftr13')
1235c$$$          call writeg1d(1,klev,ftra(1:klev,14),'ftr14','ftr14')
1236c$$$          call writeg1d(1,klev,ftra(1:klev,15),'ftr15','ftr15')
1237c$$$          call writeg1d(1,klev,ftra(1:klev,16),'ftr16','ftr16')
1238c$$$          call writeg1d(1,klev,ftra(1:klev,17),'ftr17','ftr17')
1239c$$$          call writeg1d(1,klev,ftra(1:klev,18),'ftr18','ftr18')
1240c$$$          call writeg1d(1,klev,ftra(1:klev,19),'ftr19','ftr19')
1241c$$$          call writeg1d(1,klev,ftra(1:klev,20),'ftr20','ftr20 ')
1242c$$$          call writeg1d(1,klev,mp,'mp  ','mp ')
1243c$$$          call writeg1d(1,klev,Mke,'Mke  ','Mke ')
1244
1245 
1246 
1247ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1248c
1249c
1250        RETURN
1251        END
1252C ---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.