source: trunk/libf/phylmd/convect3.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

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