source: LMDZ4/trunk/libf/phylmd/convect3.F @ 1085

Last change on this file since 1085 was 766, checked in by Laurent Fairhead, 18 years ago

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