source: trunk/libf/phylmd/cv3p1_closure.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

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