source: LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F @ 408

Last change on this file since 408 was 408, checked in by lmdzadmin, 22 years ago

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