/[lmdze]/trunk/libf/phylmd/convect3.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/convect3.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 10 months ago) by guez
File size: 45918 byte(s)
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

Deleted local variables "dtime", "co2_ppm_etat0", "solaire_etat0",
"length", "tabcntr0" in "physiq". Replaced all references to "dtime"
by references to "pdtphys".

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

  ViewVC Help
Powered by ViewVC 1.1.21