source: LMDZ6/branches/Amaury_dev/libf/phylmd/convect3.F90 @ 5218

Last change on this file since 5218 was 5160, checked in by abarral, 5 months ago

Put .h into modules

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