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

Last change on this file since 491 was 433, checked in by lmdzadmin, 21 years ago

Convergence avec la version de Ionela dec 2002

YOMCST.? : suppression RI0 (IM)
albedo.F : facteur 1.2 sur le nouveau calcul (IM)
clesphys.h : rajout de différentes ctes (concentration des gaz) (IM)
clmain.F : separation des flux LW, SW (JLD)

remplace qsurf par yqsol (IM)

conf_phys.F90 : rajout de différentes ctes (gaz + orbite) (IM)
convect3.F : DPINV+SIGD*0.5*(EVAP(1)+EVAP(2)) (SBL)
cv3_routines.F:
cvparam3.h : compatibilite avec conema3 TEMPORAIRE (FH)
phyetat0.F : lecture de co2_ppm et solaire pour tests de coherence
phyredem.F : co2_ppm et solaire passé en common
physiq.F : separation flux LW, SW

rajout diagnostiques (slp, w500)
suppression iflag_con = 4
clwcon0=qcondc (FH)
position dU "ENDIF ! ok_cvl"

radlwsw.F : passage des concentrations gaz dans un common (IM)

PEMIS(i) = 1.0 (JLD pour cohérence ORCHIDEE)

stdlevvar.F90 :
suphec.F : suppression init. des ctes orbitales (IM)

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