/[lmdze]/trunk/phylmd/convect3.f90
ViewVC logotype

Annotation of /trunk/phylmd/convect3.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 81 - (hide annotations)
Wed Mar 5 14:38:41 2014 UTC (10 years, 2 months ago) by guez
File size: 43344 byte(s) 
 Converted to free source form files which were still in fixed source
form. The conversion was done using the polish mode of the NAG Fortran
Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

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