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

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago

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