source: LMDZ5/branches/testing/libf/phylmd/convect3.F90 @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago

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