source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cv3p1_closure.F90 @ 3809

Last change on this file since 3809 was 3809, checked in by ymipsl, 10 years ago

Add LMDZ in aquaplanet configuration
YM
File size: 17.8 KB
Line 
1
2! $Id: cv3p1_closure.F90 2224 2015-03-11 11:23:40Z crio $
3
4SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
5    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
6    iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, &
7    wbeff)
8
9
10  ! **************************************************************
11  ! *
12  ! CV3P1_CLOSURE                                               *
13  ! Ale & Alp Closure of Convect3              *
14  ! *
15  ! written by   :   Kerry Emanuel                              *
16  ! vectorization:   S. Bony                                    *
17  ! modified by :    Jean-Yves Grandpeix, 18/06/2003, 19.32.10  *
18  ! Julie Frohwirth,     14/10/2005  17.44.22  *
19  ! **************************************************************
20
21  IMPLICIT NONE
22
23  include "cvthermo.h"
24  include "cv3param.h"
25  include "YOMCST2.h"
26  include "YOMCST.h"
27  include "conema3.h"
28  include "iniprint.h"
29
30  ! input:
31  INTEGER ncum, nd, nloc
32  INTEGER icb(nloc), inb(nloc)
33  REAL pbase(nloc), plcl(nloc)
34  REAL p(nloc, nd), ph(nloc, nd+1)
35  REAL tv(nloc, nd), tvp(nloc, nd), buoy(nloc, nd)
36  REAL supmax(nloc, nd)
37  LOGICAL ok_inhib ! enable convection inhibition by dryness
38  REAL ale(nloc), alp(nloc)
39  REAL omega(nloc,nd)
40
41  ! input/output:
42  REAL sig(nloc, nd), w0(nloc, nd), ptop2(nloc)
43
44  ! output:
45  REAL cape(nloc), cin(nloc)
46  REAL m(nloc, nd)
47  REAL plim1(nloc), plim2(nloc)
48  REAL asupmax(nloc, nd), supmax0(nloc)
49  REAL asupmaxmin(nloc)
50  REAL cbmf(nloc), plfc(nloc)
51  REAL wbeff(nloc)
52  INTEGER iflag(nloc)
53
54  ! local variables:
55  INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
56  REAL deltap, fac, w, amu
57  REAL rhodp
58  REAL pbmxup
59  REAL dtmin(nloc, nd), sigold(nloc, nd)
60  REAL coefmix(nloc, nd)
61  REAL pzero(nloc), ptop2old(nloc)
62  REAL cina(nloc), cinb(nloc)
63  INTEGER ibeg(nloc)
64  INTEGER nsupmax(nloc)
65  REAL supcrit, temp(nloc, nd)
66  REAL p1(nloc), pmin(nloc)
67  REAL asupmax0(nloc)
68  LOGICAL ok(nloc)
69  REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd)
70  REAL wb2(nloc)
71  REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc)
72  REAL cbmflast(nloc)
73  REAL coef(nloc)
74  REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc)
75  REAL theta(nloc), bb(nloc)
76  REAL term1, term2, term3
77  REAL alp2(nloc) ! Alp with offset
78
79  REAL sigmax
80  PARAMETER (sigmax=0.1)
81
82  CHARACTER (LEN=20) :: modname = 'cv3p1_closure'
83  CHARACTER (LEN=80) :: abort_message
84
85  ! print *,' -> cv3p1_closure, Ale ',ale(1)
86
87
88  ! -------------------------------------------------------
89  ! -- Initialization
90  ! -------------------------------------------------------
91
92
93
94  DO il = 1, ncum
95    alp2(il) = max(alp(il), 1.E-5)
96    ! IM
97    alp2(il) = max(alp(il), 1.E-12)
98  END DO
99
100  pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
101  ! exist (if any)
102
103  IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, &
104    ncum, nd, icb(nloc), inb(nloc), nl
105  DO k = 1, nl
106    DO il = 1, ncum
107      m(il, k) = 0.0
108    END DO
109  END DO
110
111  ! -------------------------------------------------------
112  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
113  ! -------------------------------------------------------
114
115  ! update sig and w0 above LNB:
116
117  DO k = 1, nl - 1
118    DO il = 1, ncum
119      IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
120        sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il &
121          ,inb(il)))
122        sig(il, k) = amax1(sig(il,k), 0.0)
123        w0(il, k) = beta*w0(il, k)
124      END IF
125    END DO
126  END DO
127
128  ! if(prt.level.GE.20) print*,'cv3p1_param apres 100'
129  ! compute icbmax:
130
131  icbmax = 2
132  DO il = 1, ncum
133    icbmax = max(icbmax, icb(il))
134  END DO
135  ! if(prt.level.GE.20) print*,'cv3p1_param apres 200'
136
137  ! update sig and w0 below cloud base:
138
139  DO k = 1, icbmax
140    DO il = 1, ncum
141      IF (k<=icb(il)) THEN
142        sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il, &
143          icb(il))
144        sig(il, k) = amax1(sig(il,k), 0.0)
145        w0(il, k) = beta*w0(il, k)
146      END IF
147    END DO
148  END DO
149  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300'
150  ! -------------------------------------------------------------
151  ! -- Reset fractional areas of updrafts and w0 at initial time
152  ! -- and after 10 time steps of no convection
153  ! -------------------------------------------------------------
154
155  DO k = 1, nl - 1
156    DO il = 1, ncum
157      IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
158        sig(il, k) = 0.0
159        w0(il, k) = 0.0
160      END IF
161    END DO
162  END DO
163  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400'
164
165  ! -------------------------------------------------------------
166  ! jyg1
167  ! --  Calculate adiabatic ascent top pressure (ptop)
168  ! -------------------------------------------------------------
169
170
171  ! c 1. Start at first level where precipitations form
172  DO il = 1, ncum
173    pzero(il) = plcl(il) - pbcrit
174  END DO
175
176  ! c 2. Add offset
177  DO il = 1, ncum
178    pzero(il) = pzero(il) - pbmxup
179  END DO
180  DO il = 1, ncum
181    ptop2old(il) = ptop2(il)
182  END DO
183
184  DO il = 1, ncum
185    ! CR:c est quoi ce 300??
186    p1(il) = pzero(il) - 300.
187  END DO
188
189  ! compute asupmax=abs(supmax) up to lnm+1
190
191  DO il = 1, ncum
192    ok(il) = .TRUE.
193    nsupmax(il) = inb(il)
194  END DO
195
196  DO i = 1, nl
197    DO il = 1, ncum
198      IF (i>icb(il) .AND. i<=inb(il)) THEN
199        IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
200          nsupmax(il) = i
201          ok(il) = .FALSE.
202        END IF ! end IF (P(i) ...  )
203      END IF ! end IF (icb+1 le i le inb)
204    END DO
205  END DO
206
207  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.'
208  DO i = 1, nl
209    DO il = 1, ncum
210      asupmax(il, i) = abs(supmax(il,i))
211    END DO
212  END DO
213
214
215  DO il = 1, ncum
216    asupmaxmin(il) = 10.
217    pmin(il) = 100.
218    ! IM ??
219    asupmax0(il) = 0.
220  END DO
221
222  ! c 3.  Compute in which level is Pzero
223
224  ! IM bug      i0 = 18
225  DO il = 1, ncum
226    i0(il) = nl
227  END DO
228
229  DO i = 1, nl
230    DO il = 1, ncum
231      IF (i>icb(il) .AND. i<=inb(il)) THEN
232        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
233          IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
234            i0(il) = i
235          END IF
236        END IF
237      END IF
238    END DO
239  END DO
240  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 3.'
241
242  ! c 4.  Compute asupmax at Pzero
243
244  DO i = 1, nl
245    DO il = 1, ncum
246      IF (i>icb(il) .AND. i<=inb(il)) THEN
247        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
248          asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( &
249            pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, &
250            i0(il)-1))
251        END IF
252      END IF
253    END DO
254  END DO
255
256
257  DO i = 1, nl
258    DO il = 1, ncum
259      IF (p(il,i)==pzero(il)) THEN
260        asupmax(i, il) = asupmax0(il)
261      END IF
262    END DO
263  END DO
264  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.'
265
266  ! c 5. Compute asupmaxmin, minimum of asupmax
267
268  DO i = 1, nl
269    DO il = 1, ncum
270      IF (i>icb(il) .AND. i<=inb(il)) THEN
271        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
272          IF (asupmax(il,i)<asupmaxmin(il)) THEN
273            asupmaxmin(il) = asupmax(il, i)
274            pmin(il) = p(il, i)
275          END IF
276        END IF
277      END IF
278    END DO
279  END DO
280
281  DO il = 1, ncum
282    ! IM
283    IF (prt_level>=20) THEN
284      PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
285        asupmaxmin(il), pzero(il), pmin(il)
286    END IF
287    IF (asupmax0(il)<asupmaxmin(il)) THEN
288      asupmaxmin(il) = asupmax0(il)
289      pmin(il) = pzero(il)
290    END IF
291  END DO
292  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 5.'
293
294
295  ! Compute Supmax at Pzero
296
297  DO i = 1, nl
298    DO il = 1, ncum
299      IF (i>icb(il) .AND. i<=inb(il)) THEN
300        IF (p(il,i)<=pzero(il)) THEN
301          supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, &
302            i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
303          GO TO 425
304        END IF ! end IF (P(i) ... )
305      END IF ! end IF (icb+1 le i le inb)
306    END DO
307  END DO
308
309425 CONTINUE
310  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.'
311
312  ! c 6. Calculate ptop2
313
314  DO il = 1, ncum
315    IF (asupmaxmin(il)<supcrit1) THEN
316      ptop2(il) = pmin(il)
317    END IF
318
319    IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
320      ptop2(il) = ptop2old(il)
321    END IF
322
323    IF (asupmaxmin(il)>supcrit2) THEN
324      ptop2(il) = ph(il, inb(il))
325    END IF
326  END DO
327
328  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.'
329
330  ! c 7. Compute multiplying factor for adiabatic updraught mass flux
331
332
333  IF (ok_inhib) THEN
334
335    DO i = 1, nl
336      DO il = 1, ncum
337        IF (i<=nl) THEN
338          coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( &
339            il,i))
340          coefmix(il, i) = min(coefmix(il,i), 1.)
341        END IF
342      END DO
343    END DO
344
345
346  ELSE ! when inhibition is not taken into account, coefmix=1
347
348
349
350    DO i = 1, nl
351      DO il = 1, ncum
352        IF (i<=nl) THEN
353          coefmix(il, i) = 1.
354        END IF
355      END DO
356    END DO
357
358  END IF ! ok_inhib
359  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.'
360  ! -------------------------------------------------------------------
361  ! -------------------------------------------------------------------
362
363
364  ! jyg2
365
366  ! ==========================================================================
367
368
369  ! -------------------------------------------------------------
370  ! -- Calculate convective inhibition (CIN)
371  ! -------------------------------------------------------------
372
373  ! do i=1,nloc
374  ! print*,'avant cine p',pbase(i),plcl(i)
375  ! enddo
376  ! do j=1,nd
377  ! do i=1,nloc
378  ! print*,'avant cine t',tv(i),tvp(i)
379  ! enddo
380  ! enddo
381  CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
382    cinb, plfc)
383
384  DO il = 1, ncum
385    cin(il) = cina(il) + cinb(il)
386  END DO
387  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'
388  ! -------------------------------------------------------------
389  ! --Update buoyancies to account for Ale
390  ! -------------------------------------------------------------
391
392  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
393    tvp, buoy)
394  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy'
395
396  ! -------------------------------------------------------------
397  ! -- Calculate convective available potential energy (cape),
398  ! -- vertical velocity (w), fractional area covered by
399  ! -- undilute updraft (sig), and updraft mass flux (m)
400  ! -------------------------------------------------------------
401
402  DO il = 1, ncum
403    cape(il) = 0.0
404  END DO
405
406  ! compute dtmin (minimum buoyancy between ICB and given level k):
407
408  DO k = 1, nl
409    DO il = 1, ncum
410      dtmin(il, k) = 100.0
411    END DO
412  END DO
413
414  DO k = 1, nl
415    DO j = minorig, nl
416      DO il = 1, ncum
417        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= &
418            (k-1))) THEN
419          dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
420        END IF
421      END DO
422    END DO
423  END DO
424
425  ! the interval on which cape is computed starts at pbase :
426
427  DO k = 1, nl
428    DO il = 1, ncum
429
430      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
431
432        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
433        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
434        cape(il) = amax1(0.0, cape(il))
435        sigold(il, k) = sig(il, k)
436
437
438        ! jyg       Coefficient coefmix limits convection to levels where a
439        ! sufficient
440        ! fraction of mixed draughts are ascending.
441        siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
442        siglim(il, k) = amax1(siglim(il,k), 0.0)
443        siglim(il, k) = amin1(siglim(il,k), 0.01)
444        ! c         fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
445        fac = 1.
446        wlim(il, k) = fac*sqrt(cape(il))
447        amu = siglim(il, k)*wlim(il, k)
448        rhodp = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
449        mlim(il, k) = amu*rhodp
450        ! print*, 'siglim ', k,siglim(1,k)
451      END IF
452
453    END DO
454  END DO
455  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600'
456
457  DO il = 1, ncum
458    ! IM beg
459    IF (prt_level>=20) THEN
460      PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
461        mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
462        ph(il, icb(il)+2)
463    END IF
464
465    IF (icb(il)+1<=inb(il)) THEN
466      ! IM end
467      mlim(il, icb(il)) = 0.5*mlim(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb( &
468        il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2))
469      ! IM beg
470    END IF !(icb(il.le.inb(il))) then
471    ! IM end
472  END DO
473  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700'
474
475  ! jyg1
476  ! ------------------------------------------------------------------------
477  ! c     Correct mass fluxes so that power used to overcome CIN does not
478  ! c     exceed Power Available for Lifting (PAL).
479  ! ------------------------------------------------------------------------
480
481  DO il = 1, ncum
482    cbmflim(il) = 0.
483    cbmf(il) = 0.
484  END DO
485
486  ! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim)
487
488  DO k = 1, nl
489    DO il = 1, ncum
490      ! old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
491      ! IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
492      IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
493          .AND. icb(il)+1<=inb(il)) THEN !cor jyg
494        cbmflim(il) = cbmflim(il) + mlim(il, k)
495      END IF
496    END DO
497  END DO
498  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim'
499
500  ! c 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
501  ! c     allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
502  ! c     Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud)
503  ! is
504  ! --    exceedingly small.
505
506  DO il = 1, ncum
507    wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
508  END DO
509
510  DO il = 1, ncum
511    IF (plfc(il)<100.) THEN
512      ! This is an irealistic value for plfc => no calculation of wbeff
513      wbeff(il) = 100.1
514    ELSE
515      ! Calculate wbeff
516      IF (flag_wb==0) THEN
517        wbeff(il) = wbmax
518      ELSE IF (flag_wb==1) THEN
519        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
520      ELSE IF (flag_wb==2) THEN
521        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
522      END IF
523    END IF
524  END DO
525
526!CR:Compute k at plfc
527  DO il=1,ncum
528           klfc(il)=nl
529  ENDDO
530  DO k=1,nl
531     DO il=1,ncum
532        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
533           klfc(il)=k
534        endif
535     ENDDO
536  ENDDO
537!RC
538
539  DO il = 1, ncum
540    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
541    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
542    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
543!CR: Add large-scale component to the mass-flux
544!encore connu sous le nom "Experience du tube de dentifrice"
545    if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
546       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
547    endif
548!RC
549    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
550      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
551        , il, alp2(il), alp(il), cin(il)
552      abort_message = ''
553      CALL abort_gcm(modname, abort_message, 1)
554    END IF
555    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
556  END DO
557
558  DO il = 1, ncum
559    IF (cbmflim(il)>1.E-6) THEN
560      ! ATTENTION TEST CR
561      ! if (cbmfmax(il).lt.1.e-12) then
562      cbmf(il) = min(cbmf1(il), cbmfmax(il))
563      ! else
564      ! cbmf(il) = cbmf1(il)
565      ! endif
566      ! print*,'cbmf',cbmf1(il),cbmfmax(il)
567    END IF
568  END DO
569  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
570
571  ! c 2. Compute coefficient and apply correction
572
573  DO il = 1, ncum
574    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
575  END DO
576  IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
577
578  DO k = 1, nl
579    DO il = 1, ncum
580      IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
581        amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)* &
582          wlim(il, k)
583        w0(il, k) = wlim(il, k)
584        w0(il, k) = max(w0(il,k), 1.E-10)
585        sig(il, k) = amu/w0(il, k)
586        sig(il, k) = min(sig(il,k), 1.)
587        ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
588        m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
589      END IF
590    END DO
591  END DO
592  ! jyg2
593  DO il = 1, ncum
594    w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
595    m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
596      (ph(il,icb(il)+1)-ph(il,icb(il)+2))
597    sig(il, icb(il)) = sig(il, icb(il)+1)
598    sig(il, icb(il)-1) = sig(il, icb(il))
599  END DO
600  IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
601
602  ! c 3. Compute final cloud base mass flux and set iflag to 3 if
603  ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if
604  ! c    the final mass flux (cbmflast) is greater than the target mass flux
605  ! c    (cbmf)).
606
607  DO il = 1, ncum
608    cbmflast(il) = 0.
609  END DO
610
611  DO k = 1, nl
612    DO il = 1, ncum
613      IF (k>=icb(il) .AND. k<=inb(il)) THEN
614          !IMpropo??      IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
615        cbmflast(il) = cbmflast(il) + m(il, k)
616      END IF
617    END DO
618  END DO
619
620  DO il = 1, ncum
621    IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
622      iflag(il) = 3
623    END IF
624  END DO
625
626  DO k = 1, nl
627    DO il = 1, ncum
628      IF (iflag(il)>=3) THEN
629        m(il, k) = 0.
630        sig(il, k) = 0.
631        w0(il, k) = 0.
632      END IF
633    END DO
634  END DO
635  IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'
636
637  ! c 4. Introduce a correcting factor for coef, in order to obtain an
638  ! effective
639  ! c    sigdz larger in the present case (using cv3p1_closure) than in the
640  ! old
641  ! c    closure (using cv3_closure).
642  IF (1==0) THEN
643    DO il = 1, ncum
644      ! c      coef(il) = 2.*coef(il)
645      coef(il) = 5.*coef(il)
646    END DO
647    ! version CVS du ..2008
648  ELSE
649    IF (iflag_cvl_sigd==0) THEN
650      ! test pour verifier qu on fait la meme chose qu avant: sid constant
651      coef(1:ncum) = 1.
652    ELSE
653      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
654      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
655    END IF
656  END IF
657
658  IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
659  RETURN
660END SUBROUTINE cv3p1_closure
661
662
Note: See TracBrowser for help on using the repository browser.