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

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

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