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

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

Inclusion du nouveau schema de nuages de SB. FH
IM/LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.4 KB
RevLine 
[207]1c $Header$
[373]2      SUBROUTINE CONVECT3       
3     *    (DTIME,EPMAX,ok_adj,
4     *     T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
[199]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,
[373]9cccc     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
10     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR,   ! sbl
11     *     FT2,FR2,FU2,FV2,WD,QCOND,QCONDC)   ! sbl
[199]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)
[373]24      REAL DELTAC              ! cld
25      PARAMETER (DELTAC=0.01)  ! cld
[199]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
[373]41      REAL FT2(ND),FR2(ND),FU2(ND),FV2(ND) ! added sbl
42      REAL U1(ND),V1(ND) ! added sbl
43C
[199]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
[373]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
[199]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
[373]124
125         FT2(I)=0.0
126         FR2(I)=0.0
127         FU2(I)=0.0
128         FV2(I)=0.0
129
[199]130         DO 4 J=1,NTRA
131          FTRA(I,J)=0.0
132    4    CONTINUE
[373]133
134         QCONDC(I)=0.0  ! cld
135         QCOND(I)=0.0   ! cld
136         NQCOND(I)=0.0  ! cld
137
[199]138         T(I)=T1(I)
139         RR(I)=R1(I)
[373]140         U1(I)=U(I) ! added sbl
141         V1(I)=V(I) ! added sbl
[199]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
[373]164      WD=0.0 ! sb
[199]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:
[373]188c     EPMAX=0.993 ! precip efficiency less than unity
[199]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
[373]224c!      NOPT=1 ! sbl
[199]225C
226C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
227C
228C     ***  DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A  ***
229C     ***                BOUNDARY LAYER SCHEME !!!               ***
230C
[373]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
[199]281C
[373]282C   ***   RESET INPUT ARRAYS IF ok_adj 0   ***
[199]283C
[373]284      IF (ok_adj)THEN
[199]285         DO 35 I=1,ND
[373]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
[199]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
[373]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
[199]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
[373]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
[199]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
[373]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
[199]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)
[373]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
[199]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***********************************************************
[373]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
[199]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
[373]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
[199]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))*
1043     1    DPINV+SIGD*0.5*(EVAP(1)+EVAP(2))
1044        FR(1)=FR(1)+0.1*AM*(RR(2)-RR(1))*DPINV
1045        FU(1)=FU(1)+0.1*DPINV*(MP(2)*(UP(2)-U(1))+AM*(U(2)-U(1)))
1046        FV(1)=FV(1)+0.1*DPINV*(MP(2)*(VP(2)-V(1))+AM*(V(2)-V(1)))
1047        DO J=1,NTRA
1048         FTRA(1,J)=FTRA(1,J)+0.1*DPINV*(MP(2)*(TRAP(2,J)-TRA(1,J))+
1049     1    AM*(TRA(2,J)-TRA(1,J)))
1050        END DO
1051        AMDE=0.0
1052        DO 415 J=2,INB
1053         FR(1)=FR(1)+0.1*DPINV*MENT(J,1)*(QENT(J,1)-RR(1))
1054         FU(1)=FU(1)+0.1*DPINV*MENT(J,1)*(UENT(J,1)-U(1))
1055         FV(1)=FV(1)+0.1*DPINV*MENT(J,1)*(VENT(J,1)-V(1))
1056         DO K=1,NTRA
1057          FTRA(1,K)=FTRA(1,K)+0.1*DPINV*MENT(J,1)*(TRAENT(J,1,K)-
1058     1     TRA(1,K))
1059         END DO
1060  415      CONTINUE
1061C
1062C   ***  CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO  ***
1063C   ***               AT LEVELS ABOVE THE LOWEST LEVEL                   ***
1064C
1065C   ***  FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES  ***
1066C   ***                      THROUGH EACH LEVEL                          ***
1067C
1068 
1069 
1070        DO 500 I=2,INB
1071        DPINV=1.0/(PH(I)-PH(I+1))
1072      CPINV=1.0/CPN(I)
1073        AMP1=0.0
1074        DO 440 K=I+1,INB+1
1075  440   AMP1=AMP1+M(K)
1076        DO 450 K=1,I
1077        DO 450 J=I+1,INB+1
1078         AMP1=AMP1+MENT(K,J)
1079  450   CONTINUE
1080      IF((0.1*DPINV*AMP1).GE.DELTI)IFLAG=4
1081        AD=0.0
1082        DO 470 K=1,I-1
1083        DO 470 J=I,INB
1084  470   AD=AD+MENT(J,K)
1085      FT(I)=0.1*DPINV*(AMP1*(T(I+1)-T(I)+(GZ(I+1)-GZ(I))*
1086     1   CPINV)-AD*(T(I)-T(I-1)+(GZ(I)-GZ(I-1))*CPINV))
1087     2   -0.5*SIGD*LVCP(I)*(EVAP(I)+EVAP(I+1))
1088      RAT=CPN(I-1)*CPINV
1089        FT(I)=FT(I)-0.09*SIGD*(MP(I+1)*T(I)*
1090     1    B(I)-MP(I)*T(I-1)*RAT*B(I-1))*DPINV
1091      FT(I)=FT(I)+0.1*DPINV*MENT(I,I)*(HP(I)-H(I)+
1092     1    T(I)*(CPV-CPD)*(RR(I)-QENT(I,I)))*CPINV
1093      FT(I)=FT(I)+0.01*SIGD*WT(I)*(CL-CPD)*WATER(I+1)*
1094     1    (T(I+1)-T(I))*DPINV*CPINV
1095        FR(I)=0.1*DPINV*(AMP1*(RR(I+1)-RR(I))-
1096     1    AD*(RR(I)-RR(I-1)))
1097        FU(I)=FU(I)+0.1*DPINV*(AMP1*(U(I+1)-U(I))-
1098     1    AD*(U(I)-U(I-1)))
1099        FV(I)=FV(I)+0.1*DPINV*(AMP1*(V(I+1)-V(I))-
1100     1    AD*(V(I)-V(I-1)))
1101        DO K=1,NTRA
1102         FTRA(I,K)=FTRA(I,K)+0.1*DPINV*(AMP1*(TRA(I+1,K)-
1103     1    TRA(I,K))-AD*(TRA(I,K)-TRA(I-1,K)))
1104        END DO
1105        DO 480 K=1,I-1
1106       AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I)
1107       AWAT=AMAX1(AWAT,0.0)
1108         FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-AWAT
1109     1    -RR(I))
1110         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
1111         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
[373]1112C (saturated updrafts resulting from mixing)      ! cld   
1113         QCOND(I)=QCOND(I)+(ELIJ(K,I)-AWAT)       ! cld
1114         NQCOND(I)=NQCOND(I)+1.                   ! cld
[199]1115         DO J=1,NTRA
1116          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
1117     1     TRA(I,J))
1118         END DO
1119  480   CONTINUE
1120      DO 490 K=I,INB
1121       FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-RR(I))
1122         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
1123         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
1124         DO J=1,NTRA
1125          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
1126     1     TRA(I,J))
1127         END DO
1128  490      CONTINUE
1129        FR(I)=FR(I)+0.5*SIGD*(EVAP(I)+EVAP(I+1))+0.1*(MP(I+1)*
1130     1    (RP(I+1)-RR(I))-MP(I)*(RP(I)-RR(I-1)))*DPINV
1131        FU(I)=FU(I)+0.1*(MP(I+1)*(UP(I+1)-U(I))-MP(I)*
1132     1    (UP(I)-U(I-1)))*DPINV
1133        FV(I)=FV(I)+0.1*(MP(I+1)*(VP(I+1)-V(I))-MP(I)*
1134     1    (VP(I)-V(I-1)))*DPINV
1135        DO J=1,NTRA
1136         FTRA(I,J)=FTRA(I,J)+0.1*DPINV*(MP(I+1)*(TRAP(I+1,J)-TRA(I,J))-
1137     1    MP(I)*(TRAP(I,J)-TRAP(I-1,J)))
1138        END DO
[373]1139C (saturated downdrafts resulting from mixing)    ! cld
1140        DO K=I+1,INB                              ! cld
1141         QCOND(I)=QCOND(I)+ELIJ(K,I)              ! cld
1142         NQCOND(I)=NQCOND(I)+1.                   ! cld
1143        ENDDO                                     ! cld
1144C (particular case: no detraining level is found) ! cld
1145        IF (NENT(I).EQ.0) THEN                    ! cld
1146         QCOND(I)=QCOND(I)+(1-EP(I))*CLW(I)       ! cld
1147         NQCOND(I)=NQCOND(I)+1.                   ! cld
1148        ENDIF                                     ! cld
1149        IF (NQCOND(I).NE.0.) THEN                 ! cld
1150         QCOND(I)=QCOND(I)/NQCOND(I)              ! cld
1151        ENDIF                                     ! cld
[199]1152  500   CONTINUE
1153 
1154 
1155 
1156C
1157C   ***   MOVE THE DETRAINMENT AT LEVEL INB DOWN TO LEVEL INB-1   ***
1158C   ***        IN SUCH A WAY AS TO PRESERVE THE VERTICALLY        ***
1159C   ***          INTEGRATED ENTHALPY AND WATER TENDENCIES         ***
1160C
[373]1161c test sb:
1162c@      write(*,*) '--------------------------------------------'
1163c@      write(*,*) 'inb,ft,hp,h,t,rr,qent,ment,water,waterp,wt,mp,b'
1164c@      write(*,*) inb,ft(inb),hp(inb),h(inb)
1165c@     :   ,t(inb),rr(inb),qent(inb,inb)
1166c@     :   ,ment(inb,inb),water(inb)
1167c@     :   ,water(inb+1),wt(inb),mp(inb),b(inb)
1168c@      write(*,*) '--------------------------------------------'
1169c fin test sb:
1170
[199]1171      AX=0.1*MENT(INB,INB)*(HP(INB)-H(INB)+T(INB)*
1172     1    (CPV-CPD)*(RR(INB)-QENT(INB,INB)))/(CPN(INB)*
1173     2    (PH(INB)-PH(INB+1)))
1174      FT(INB)=FT(INB)-AX
1175      FT(INB-1)=FT(INB-1)+AX*CPN(INB)*(PH(INB)-PH(INB+1))/
1176     1    (CPN(INB-1)*(PH(INB-1)-PH(INB)))
1177      BX=0.1*MENT(INB,INB)*(QENT(INB,INB)-RR(INB))/
1178     1    (PH(INB)-PH(INB+1))
1179      FR(INB)=FR(INB)-BX
1180      FR(INB-1)=FR(INB-1)+BX*(PH(INB)-PH(INB+1))/
1181     1    (PH(INB-1)-PH(INB))
1182      CX=0.1*MENT(INB,INB)*(UENT(INB,INB)-U(INB))/
1183     1    (PH(INB)-PH(INB+1))
1184      FU(INB)=FU(INB)-CX
1185      FU(INB-1)=FU(INB-1)+CX*(PH(INB)-PH(INB+1))/
1186     1    (PH(INB-1)-PH(INB))
1187      DX=0.1*MENT(INB,INB)*(VENT(INB,INB)-V(INB))/
1188     1    (PH(INB)-PH(INB+1))
1189      FV(INB)=FV(INB)-DX
1190      FV(INB-1)=FV(INB-1)+DX*(PH(INB)-PH(INB+1))/
1191     1    (PH(INB-1)-PH(INB))
1192      DO J=1,NTRA
1193      EX=0.1*MENT(INB,INB)*(TRAENT(INB,INB,J)
1194     1    -TRA(INB,J))/(PH(INB)-PH(INB+1))
1195      FTRA(INB,J)=FTRA(INB,J)-EX
1196      FTRA(INB-1,J)=FTRA(INB-1,J)+EX*
1197     1     (PH(INB)-PH(INB+1))/(PH(INB-1)-PH(INB))
1198      ENDDO   
1199C
1200C   ***    HOMOGINIZE TENDENCIES BELOW CLOUD BASE    ***
1201C
1202      ASUM=0.0
1203      BSUM=0.0
1204      CSUM=0.0
1205        DSUM=0.0
1206      DO 650 I=1,ICB-1
1207       ASUM=ASUM+FT(I)*(PH(I)-PH(I+1))
1208         BSUM=BSUM+FR(I)*(LV(I)+(CL-CPD)*(T(I)-T(1)))*
1209     1    (PH(I)-PH(I+1))
1210       CSUM=CSUM+(LV(I)+(CL-CPD)*(T(I)-T(1)))*(PH(I)-PH(I+1))
1211       DSUM=DSUM+T(I)*(PH(I)-PH(I+1))/TH(I)
1212  650      CONTINUE
1213      DO 700 I=1,ICB-1
1214       FT(I)=ASUM*T(I)/(TH(I)*DSUM)
1215       FR(I)=BSUM/CSUM
1216  700      CONTINUE
1217C
1218C   ***           RESET COUNTER AND RETURN           ***
1219C
1220      SIG(ND)=2.0
1221c
1222c
1223      do i = 1, nd
1224         upwd(i) = 0.0
1225         dnwd(i) = 0.0
1226c sb       dnwd0(i) = - mp(i)
1227      enddo
1228c
1229      do i = 1, nl
1230       dnwd0(i) = - mp(i)
1231      enddo
1232      do i = nl+1, nd
1233       dnwd0(i) = 0.
1234      enddo
1235c
1236      do i = icb, inb
1237         upwd(i) = 0.0
1238         dnwd(i) = 0.0
1239
1240         do k =i, inb
1241            up1=0.0
1242            dn1=0.0
1243            do n = 1, i-1
1244               up1 = up1 + ment(n,k)
1245               dn1 = dn1 - ment(k,n)
1246            enddo
1247            upwd(i) = upwd(i)+ m(k) + up1
1248            dnwd(i) = dnwd(i) + dn1
1249         enddo
1250        enddo
1251 
1252cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1253c        DETERMINATION DE LA VARIATION DE FLUX ASCENDANT ENTRE
1254C        DEUX NIVEAU NON DILUE Mike
1255cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1256 
1257 
1258c sb      do i=1,ND
1259c sb      Mike(i)=M(i)
1260c sb      enddo
1261 
1262      do i = 1, NL
1263       Mike(i) = M(i)
1264      enddo
1265      do i = NL+1, ND
1266       Mike(i) = 0.
1267      enddo
1268 
1269      do i=1,nd
1270      Ma(i)=0
1271      enddo
1272 
1273c sb      do i=1,nd
1274c sb      do j=i,nd
1275c sb      Ma(i)=Ma(i)+M(j)
1276c sb      enddo
1277c sb      enddo
1278
1279      do i = 1, NL
1280      do j = i, NL
1281       Ma(i) = Ma(i) + M(j)
1282      enddo
1283      enddo
1284c
1285      do i = NL+1, ND
1286       Ma(i) = 0.
1287      enddo
1288c
1289      do i=1,ICB-1
1290      Ma(i)=0
1291      enddo
1292 
1293 
1294 
1295ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1296C        ICB REPRESENTE DE NIVEAU OU SE TROUVE LA
1297c        BASE DU NUAGE , ET INB LE TOP DU NUAGE
1298cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1299 
1300 
1301       do i=1,ND
1302       Mke(i)=upwd(i)+dnwd(i)
1303       enddo
[373]1304
1305C
1306C   *** Diagnose the in-cloud mixing ratio   ***              ! cld
1307C   ***           of condensed water         ***              ! cld
1308C                                                             ! cld
1309       DO I=1,ND                                              ! cld
1310        MAA(I)=0.0                                            ! cld
1311        WA(I)=0.0                                             ! cld
1312        SIGA(I)=0.0                                           ! cld
1313       ENDDO                                                  ! cld
1314       DO I=NK,INB                                            ! cld
1315       DO K=I+1,INB+1                                         ! cld
1316        MAA(I)=MAA(I)+M(K)                                    ! cld
1317       ENDDO                                                  ! cld
1318       ENDDO                                                  ! cld
1319       DO I=ICB,INB-1                                         ! cld
1320        AXC(I)=0.                                             ! cld
1321        DO J=ICB,I                                            ! cld
1322         AXC(I)=AXC(I)+RD*(TVP(J)-TV(J))*(PH(J)-PH(J+1))/P(J) ! cld
1323        ENDDO                                                 ! cld
1324        IF (AXC(I).GT.0.0) THEN                               ! cld
1325         WA(I)=SQRT(2.*AXC(I))                                ! cld
1326        ENDIF                                                 ! cld
1327       ENDDO                                                  ! cld
1328       DO I=1,NL                                              ! cld
1329        IF (WA(I).GT.0.0)                                     ! cld
1330     1    SIGA(I)=MAA(I)/WA(I)*RD*TVP(I)/P(I)/100./DELTAC     ! cld
1331        SIGA(I) = MIN(SIGA(I),1.0)                            ! cld
1332        QCONDC(I)=SIGA(I)*CLW(I)*(1.-EP(I))                   ! cld
1333     1          + (1.-SIGA(I))*QCOND(I)                       ! cld
1334       ENDDO                                                  ! cld
1335
1336
[199]1337c$$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1338c$$$         call writeg1d(1,klev,ma,'ma  ','ma  ')
1339c$$$          call writeg1d(1,klev,upwd,'upwd  ','upwd  ')
1340c$$$          call writeg1d(1,klev,dnwd,'dnwd  ','dnwd  ')
1341c$$$          call writeg1d(1,klev,dnwd0,'dnwd0  ','dnwd0  ')
1342c$$$          call writeg1d(1,klev,tvp,'tvp  ','tvp  ')
1343c$$$          call writeg1d(1,klev,tra(1:klev,3),'tra3  ','tra3  ')
1344c$$$          call writeg1d(1,klev,tra(1:klev,4),'tra4  ','tra4  ')
1345c$$$          call writeg1d(1,klev,tra(1:klev,5),'tra5  ','tra5  ')
1346c$$$          call writeg1d(1,klev,tra(1:klev,6),'tra6  ','tra6  ')
1347c$$$          call writeg1d(1,klev,tra(1:klev,7),'tra7  ','tra7  ')
1348c$$$          call writeg1d(1,klev,tra(1:klev,8),'tra8  ','tra8  ')
1349c$$$          call writeg1d(1,klev,tra(1:klev,9),'tra9  ','tra9  ')
1350c$$$          call writeg1d(1,klev,tra(1:klev,10),'tra10','tra10')
1351c$$$          call writeg1d(1,klev,tra(1:klev,11),'tra11','tra11')
1352c$$$          call writeg1d(1,klev,tra(1:klev,12),'tra12','tra12')
1353c$$$          call writeg1d(1,klev,tra(1:klev,13),'tra13','tra13')
1354c$$$          call writeg1d(1,klev,tra(1:klev,14),'tra14','tra14')
1355c$$$          call writeg1d(1,klev,tra(1:klev,15),'tra15','tra15')
1356c$$$          call writeg1d(1,klev,tra(1:klev,16),'tra16','tra16')
1357c$$$          call writeg1d(1,klev,tra(1:klev,17),'tra17','tra17')
1358c$$$          call writeg1d(1,klev,tra(1:klev,18),'tra18','tra18')
1359c$$$          call writeg1d(1,klev,tra(1:klev,19),'tra19','tra19')
1360c$$$          call writeg1d(1,klev,tra(1:klev,20),'tra20','tra20 ')
1361c$$$          call writeg1d(1,klev,trap(1:klev,1),'trp1','trp1')
1362c$$$          call writeg1d(1,klev,trap(1:klev,2),'trp2','trp2')
1363c$$$          call writeg1d(1,klev,trap(1:klev,3),'trp3','trp3')
1364c$$$          call writeg1d(1,klev,trap(1:klev,4),'trp4','trp4')
1365c$$$          call writeg1d(1,klev,trap(1:klev,5),'trp5','trp5')
1366c$$$          call writeg1d(1,klev,trap(1:klev,10),'trp10','trp10')
1367c$$$          call writeg1d(1,klev,trap(1:klev,12),'trp12','trp12')
1368c$$$          call writeg1d(1,klev,trap(1:klev,15),'trp15','trp15')
1369c$$$          call writeg1d(1,klev,trap(1:klev,20),'trp20','trp20')
1370c$$$          call writeg1d(1,klev,ftra(1:klev,1),'ftr1  ','ftr1  ')
1371c$$$          call writeg1d(1,klev,ftra(1:klev,2),'ftr2  ','ftr2  ')
1372c$$$          call writeg1d(1,klev,ftra(1:klev,3),'ftr3  ','ftr3  ')
1373c$$$          call writeg1d(1,klev,ftra(1:klev,4),'ftr4  ','ftr4  ')
1374c$$$          call writeg1d(1,klev,ftra(1:klev,5),'ftr5  ','ftr5  ')
1375c$$$          call writeg1d(1,klev,ftra(1:klev,6),'ftr6  ','ftr6  ')
1376c$$$          call writeg1d(1,klev,ftra(1:klev,7),'ftr7  ','ftr7  ')
1377c$$$          call writeg1d(1,klev,ftra(1:klev,8),'ftr8  ','ftr8  ')
1378c$$$          call writeg1d(1,klev,ftra(1:klev,9),'ftr9  ','ftr9  ')
1379c$$$          call writeg1d(1,klev,ftra(1:klev,10),'ftr10','ftr10')
1380c$$$          call writeg1d(1,klev,ftra(1:klev,11),'ftr11','ftr11')
1381c$$$          call writeg1d(1,klev,ftra(1:klev,12),'ftr12','ftr12')
1382c$$$          call writeg1d(1,klev,ftra(1:klev,13),'ftr13','ftr13')
1383c$$$          call writeg1d(1,klev,ftra(1:klev,14),'ftr14','ftr14')
1384c$$$          call writeg1d(1,klev,ftra(1:klev,15),'ftr15','ftr15')
1385c$$$          call writeg1d(1,klev,ftra(1:klev,16),'ftr16','ftr16')
1386c$$$          call writeg1d(1,klev,ftra(1:klev,17),'ftr17','ftr17')
1387c$$$          call writeg1d(1,klev,ftra(1:klev,18),'ftr18','ftr18')
1388c$$$          call writeg1d(1,klev,ftra(1:klev,19),'ftr19','ftr19')
1389c$$$          call writeg1d(1,klev,ftra(1:klev,20),'ftr20','ftr20 ')
1390c$$$          call writeg1d(1,klev,mp,'mp  ','mp ')
1391c$$$          call writeg1d(1,klev,Mke,'Mke  ','Mke ')
1392
1393 
1394 
1395ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1396c
1397c
1398        RETURN
1399        END
1400C ---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.