source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/cv3p1_closure.F90

Last change on this file was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.
  • 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: 17.2 KB
Line 
1
2! $Id: cv3p1_closure.F90 1992 2014-03-05 13:19:12Z evignon $
3
4SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
5    tvp, buoy, supmax, ok_inhib, ale, alp, 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
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)
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
526  DO il = 1, ncum
527    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
528    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
529    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
530    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
531      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
532        , il, alp2(il), alp(il), cin(il)
533      abort_message = ''
534      CALL abort_gcm(modname, abort_message, 1)
535    END IF
536    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
537  END DO
538
539  DO il = 1, ncum
540    IF (cbmflim(il)>1.E-6) THEN
541      ! ATTENTION TEST CR
542      ! if (cbmfmax(il).lt.1.e-12) then
543      cbmf(il) = min(cbmf1(il), cbmfmax(il))
544      ! else
545      ! cbmf(il) = cbmf1(il)
546      ! endif
547      ! print*,'cbmf',cbmf1(il),cbmfmax(il)
548    END IF
549  END DO
550  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
551
552  ! c 2. Compute coefficient and apply correction
553
554  DO il = 1, ncum
555    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
556  END DO
557  IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
558
559  DO k = 1, nl
560    DO il = 1, ncum
561      IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
562        amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)* &
563          wlim(il, k)
564        w0(il, k) = wlim(il, k)
565        w0(il, k) = max(w0(il,k), 1.E-10)
566        sig(il, k) = amu/w0(il, k)
567        sig(il, k) = min(sig(il,k), 1.)
568        ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
569        m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
570      END IF
571    END DO
572  END DO
573  ! jyg2
574  DO il = 1, ncum
575    w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
576    m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
577      (ph(il,icb(il)+1)-ph(il,icb(il)+2))
578    sig(il, icb(il)) = sig(il, icb(il)+1)
579    sig(il, icb(il)-1) = sig(il, icb(il))
580  END DO
581  IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
582
583  ! c 3. Compute final cloud base mass flux and set iflag to 3 if
584  ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if
585  ! c    the final mass flux (cbmflast) is greater than the target mass flux
586  ! c    (cbmf)).
587
588  DO il = 1, ncum
589    cbmflast(il) = 0.
590  END DO
591
592  DO k = 1, nl
593    DO il = 1, ncum
594      IF (k>=icb(il) .AND. k<=inb(il)) THEN
595          !IMpropo??      IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
596        cbmflast(il) = cbmflast(il) + m(il, k)
597      END IF
598    END DO
599  END DO
600
601  DO il = 1, ncum
602    IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
603      iflag(il) = 3
604    END IF
605  END DO
606
607  DO k = 1, nl
608    DO il = 1, ncum
609      IF (iflag(il)>=3) THEN
610        m(il, k) = 0.
611        sig(il, k) = 0.
612        w0(il, k) = 0.
613      END IF
614    END DO
615  END DO
616  IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'
617
618  ! c 4. Introduce a correcting factor for coef, in order to obtain an
619  ! effective
620  ! c    sigdz larger in the present case (using cv3p1_closure) than in the
621  ! old
622  ! c    closure (using cv3_closure).
623  IF (1==0) THEN
624    DO il = 1, ncum
625      ! c      coef(il) = 2.*coef(il)
626      coef(il) = 5.*coef(il)
627    END DO
628    ! version CVS du ..2008
629  ELSE
630    IF (iflag_cvl_sigd==0) THEN
631      ! test pour verifier qu on fait la meme chose qu avant: sid constant
632      coef(1:ncum) = 1.
633    ELSE
634      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
635      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
636    END IF
637  END IF
638
639  IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
640  RETURN
641END SUBROUTINE cv3p1_closure
642
643
Note: See TracBrowser for help on using the repository browser.