source: LMDZ4/branches/IPSL-CM4_IPCC_patches/libf/phylmd/convect3.F @ 5306

Last change on this file since 5306 was 524, checked in by lmdzadmin, 20 years ago

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