source: LMDZ5/branches/testing/libf/phylmd/cv3p1_closure.F90 @ 2435

Last change on this file since 2435 was 2435, checked in by Laurent Fairhead, 9 years ago

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