source: LMDZ4/trunk/libf/phylmd/cv3p1_closure.F @ 900

Last change on this file since 900 was 879, checked in by Laurent Fairhead, 17 years ago

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