source: LMDZ5/trunk/libf/phylmd/cv3p1_closure.F @ 1990

Last change on this file since 1990 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

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