source: LMDZ5/trunk/libf/phylmd/convect3.F90 @ 2320

Last change on this file since 2320 was 2320, checked in by Ehouarn Millour, 9 years ago

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