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