source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3p1_closure.F @ 1085

Last change on this file since 1085 was 973, checked in by lmdzadmin, 17 years ago

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