source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_cine.F @ 1276

Last change on this file since 1276 was 1127, checked in by idelkadi, 16 years ago

Corrections sur les wakes et la convection pour surmonter le probleme de l'eau negative
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
RevLine 
[879]1        SUBROUTINE cv3_cine(nloc,ncum,nd,icb,inb
2     :                      ,pbase,plcl,p,ph,tv,tvp
3     :                      ,cina,cinb)
4
5***************************************************************
6*                                                             *
7* CV3_CINE                                                    *
8*                                                             *
9*                                                             *
10* written by   :   Frederique Cheruy                          *
11* vectorization:   Jean-Yves Grandpeix, 19/06/2003, 11.54.43  *
12* modified by :                                               *
13***************************************************************
14*
15      implicit none
16c
17#include "YOMCST.h"
18#include "cvthermo.h"
19#include "cv3param.h"
20c input:
21      integer ncum, nd, nloc
22      integer icb(nloc), inb(nloc)
23      real pbase(nloc),plcl(nloc)
24      real p(nloc,nd), ph(nloc,nd+1)
25      real tv(nloc,nd),tvp(nloc,nd)
26c
27c output
28      real cina(nloc),cinb(nloc)
29c
30c local variables
31      integer il,i,j,k
32      integer itop(nloc),ineg(nloc),ilow(nloc)
33      integer ifst(nloc),isublcl(nloc)
34      logical lswitch(nloc),lswitch1(nloc),lswitch2(nloc)
[1127]35      logical exist_lfc(nloc)
36      real plfc(nloc)
[879]37      real dpmax
38      real deltap,dcin
39      real buoylcl(nloc),tvplcl(nloc),tvlcl(nloc)
[1127]40      real p0(nloc)
[879]41      real buoyz(nloc), buoy(nloc,nd)
42c
43c-------------------------------------------------------------
44c     Initialization
45c-------------------------------------------------------------
46      do il = 1,ncum
47       cina(il) = 0.
48       cinb(il) = 0.
49      enddo
50c
51c--------------------------------------------------------------
52c      Recompute buoyancies
53c--------------------------------------------------------------
[1127]54      DO k = 1,nd
[879]55        DO il = 1,ncum
[1127]56!      print*,'tvp tv=',tvp(il,k),tv(il,k)
[879]57          buoy(il,k) = tvp(il,k) - tv(il,k)
58        ENDDO
59      ENDDO
60c---------------------------------------------------------------
61c
62c   calcul de la flottabilite a LCL (Buoylcl)
63c     ifst = first P-level above lcl
64c     isublcl = highest P-level below lcl.
65c---------------------------------------------------------------
66c
67      do il = 1,ncum
68       TVPlcl(il) = TVP(il,1)*(Plcl(il)/P(il,1))**(2./7.) !For dry air, R/Cp=2/7
69      enddo
70c
71      do il = 1,ncum
72       IF (Plcl(il) .GT. P(il,icb(il))) THEN
73        ifst(il) = icb(il)
74        isublcl(il) = icb(il)-1
75       ELSE
76        ifst(il) = icb(il)+1
77        isublcl(il) = icb(il)
78       ENDIF
79      enddo
80c
81      do il = 1,ncum
82       TVlcl(il)=TV(il,ifst(il)-1)+(TV(il,ifst(il))-TV(il,ifst(il)-1))
83     $   *(Plcl(il)-P(il,ifst(il)-1))/(P(il,ifst(il))-P(il,ifst(il)-1))
84      enddo
85c
86      do il = 1,ncum
87        BUOYlcl(il) = TVPlcl(il)-TVlcl(il)
88      enddo
89c
90c---------------------------------------------------------------
[1127]91c premiere couche contenant un  niveau de flotabilite positive
92c et premiere couche contenant un  niveau de flotabilite negative
93c  au dessus du niveau de condensation
94c---------------------------------------------------------------
95      do il = 1,ncum
96        itop(il) =nl-1
97        ineg(il) = nl-1
98        exist_lfc(il) = .FALSE.
99      enddo
100      do 100 k=nl-1,1,-1
101       do 110 il=1,ncum
102        if (k .ge. ifst(il)) then
103         if (buoy(il,k) .gt. 0.) then
104          itop(il)=k
105          exist_lfc(il) = .TRUE.
106         else
107          ineg(il)=k
108         endif
109        endif
110110    continue
111100   continue
112c
113c---------------------------------------------------------------
114c When there is no positive buoyancy level, set Plfc, Cina and Cinb
115c to arbitrary extreme values.
116c---------------------------------------------------------------
117      DO il = 1,ncum
118       IF (.NOT.exist_lfc(il)) THEN
119         Plfc(il) = 1.111
120         Cinb(il) = -1111.
121         Cina(il) = -1112.
122       ENDIF
123      ENDDO
124c
125c
126c---------------------------------------------------------------
[879]127c -- Two cases : BUOYlcl >= 0 and BUOYlcl < 0.
128c---------------------------------------------------------------
129C
130C--------------------
131C -- 1.0 BUOYlcl >=0.
132C--------------------
133c
134      DPMAX = 50.
135      DO il = 1,ncum
[1127]136        lswitch1(il)=BUOYlcl(il) .GE. 0. .AND. exist_lfc(il)
[879]137        lswitch(il) = lswitch1(il)
138      ENDDO
139c
140c 1.1 No inhibition case
141c ----------------------
142C   If buoyancy is positive at LCL and stays positive over a large enough
143C pressure interval (=DPMAX), inhibition is set to zero,
144C
145      DO il = 1,ncum
146      IF (lswitch(il)) THEN
147        IF (P(il,ineg(il)) .LT. P(il,icb(il))-DPmax) THEN
148          PLFC(il) = Plcl(il)
149          Cina(il) = 0.
150          Cinb(il) = 0.
151        ENDIF
152      ENDIF
153      ENDDO
154c
155c 1.2 Upper inhibition only case
156c ------------------------------
157      DO il = 1,ncum
158        lswitch2(il)= P(il,ineg(il)) .GE. P(il,icb(il))-DPmax
159        lswitch(il) = lswitch1(il) .AND. lswitch2(il)
160      ENDDO
161c
162      DO il = 1,ncum
163      IF (lswitch(il)) THEN
164          Cinb(il) = 0.
165c
166c 1.2.1  Calcul de la pression du niveau de flot. nulle juste au-dessus de LCL
167c ---------------------------------------------------------------------------
168         IF (ineg(il) .GT. isublcl(il)+1) THEN
169C In order to get P0, one may interpolate linearly buoyancies
170C  between P(ineg) and P(ineg-1).
171        P0(il)=(buoy(il,ineg(il))*P(il,ineg(il)-1)
172     $         -buoy(il,ineg(il)-1)*P(il,ineg(il)))
173     :           / (buoy(il,ineg(il))-buoy(il,ineg(il)-1))
174         ELSE
175C In order to get P0, one has to interpolate between P(ineg) and Plcl.
176        P0(il) = (BUOY(il,ineg(il))*Plcl(il)-BUOYlcl(il)*P(il,ineg(il)))
177     $          /(BUOY(il,ineg(il))     -BUOYlcl(il))
178         ENDIF
179      ENDIF
180      ENDDO
181c
182c 1.2.2 Recompute itop (=1st layer with positive buoyancy above ineg)
183c -------------------------------------------------------------------
184      do il = 1,ncum
185      IF (lswitch(il)) THEN
186        itop(il) =nl-1
187      ENDIF
188      enddo
189c
190      do  k=nl,1,-1
191       do  il=1,ncum
192       IF (lswitch(il)) THEN
193        if (k .ge. ineg(il) .and. buoy(il,k) .gt. 0) then
194         itop(il)=k
195        endif
196       ENDIF
197       enddo
198      enddo
199c
200c 1.2.3 Computation of PLFC
201c -------------------------
202      DO il = 1,ncum
203      IF (lswitch(il)) THEN
204        PLFC(il)=(buoy(il,itop(il))*P(il,itop(il)-1)
205     $           -buoy(il,itop(il)-1)*P(il,itop(il)))
206     $           / (buoy(il,itop(il))-buoy(il,itop(il)-1))
207      ENDIF
208      ENDDO
209c
210c 1.2.4 Computation of CINA
211c -------------------------
212c
213C   Upper part of CINA : integral from P(itop-1) to Plfc
214      DO il = 1,ncum
215      IF (lswitch(il)) THEN
216        deltap = P(il,itop(il)-1)-Plfc(il)
217        dcin = RD*BUOY(il,itop(il)-1)*deltap
218     $        / (P(il,itop(il)-1)+Plfc(il))
219        CINA(il) = min(0.,dcin)
220      ENDIF
221      ENDDO
222c
223C   Middle part of CINA : integral from P(ineg) to P(itop-1)
224      DO k = 1,nl
225        DO il = 1,ncum
226        IF (lswitch(il)) THEN
227          IF (k .GE. ineg(il) .AND. k .LE. itop(il)-2) THEN
228           deltap = P(il,k)-P(il,k+1)
229           dcin = 0.5*RD*(BUOY(il,k)+BUOY(il,k+1))*deltap/PH(il,k+1)
230           CINA(il) = CINA(il) + min(0.,dcin)
231          ENDIF
232        ENDIF
233        ENDDO
234      ENDDO
235c
236C   Lower part of CINA : integral from P0 to P(ineg)
237      DO il = 1,ncum
238      IF (lswitch(il)) THEN
239        deltap = P0(il)-P(il,ineg(il))
240        dcin = RD*BUOY(il,ineg(il))*deltap/(P(il,ineg(il))+P0(il))
241        CINA(il) = CINA(il) + min(0.,dcin)
242      ENDIF
243      ENDDO
244c
245C
246C ------------------
247C -- 2.0 BUOYlcl <0.
248C ------------------
249C
250      DO il = 1,ncum
[1127]251        lswitch1(il)=BUOYlcl(il) .LT. 0. .AND. exist_lfc(il)
[879]252        lswitch(il) = lswitch1(il)
253      ENDDO
254c
255c 2.0.1 Premiere  couche ou la flotabilite est negative au dessus du sol
256c ----------------------------------------------------
[1127]257c    au cas ou elle existe  sinon ilow=1 (nk apres)
[879]258c      on suppose que la parcelle part de la premiere couche
259c
260      DO il = 1,ncum
261      IF (lswitch(il)) THEN
262       ilow(il)=1
263      ENDIF
264      ENDDO
265c
[1127]266      do 200 k=nl,1,-1
[879]267        DO il = 1,ncum
268        IF (lswitch(il) .AND. k .LE.icb(il)-1) THEN
269         if(buoy(il,k).lt. 0.) then
270           ilow(il) = k
271          endif
272        ENDIF
273        ENDDO
274 200  continue
275
276c 2.0.2  Calcul de la pression du niveau de flot. nulle sous le nuage
277c ----------------------------------------------------
278      DO il = 1,ncum
279      IF (lswitch(il)) THEN
280       if(ilow(il).gt. 1) then
281         P0(il)=(buoy(il,ilow(il))*P(il,ilow(il)-1)
282     $          -buoy(il,ilow(il)-1)*P(il,ilow(il)))
283     :            / (buoy(il,ilow(il))-buoy(il,ilow(il)-1))
284         BUOYz(il) = 0.
285       else
286         P0(il) = P(il,1)
287         BUOYz(il) = BUOY(il,1)
288       endif
289      ENDIF
290      ENDDO
291c
292C 2.1. Computation of CINB
293C -----------------------
294c
295      DO il = 1,ncum
296        lswitch2(il)= (isublcl(il) .EQ. 1 .AND. ilow(il) .EQ. 1)
297     $                  .OR.(isublcl(il) .EQ. ilow(il)-1)
298        lswitch(il) = lswitch1(il) .AND. lswitch2(il)
299      ENDDO
300cc      IF (    (isublcl .EQ. 1 .AND. ilow .EQ. 1)
301cc     $    .OR.(isublcl .EQ. ilow-1)) THEN
302c
303c 2.1.1 First case : Plcl just above P0
304c -------------------------------------
305      DO il = 1,ncum
306      IF (lswitch(il)) THEN
307        deltap = P0(il)-Plcl(il)
308        dcin = RD*(BUOYz(il)+BUOYlcl(il))*deltap/(P0(il)+Plcl(il))
309        CINB(il) = min(0.,dcin)
310      ENDIF
311      ENDDO
312c
313      DO il = 1,ncum
314        lswitch(il) = lswitch1(il) .AND. .NOT. lswitch2(il)
315      ENDDO
316cc      ELSE
317c
318c 2.1.2 Second case : there is at least one P-level between P0 and Plcl
319c ---------------------------------------------------------------------
320c
321C   Lower part of CINB : integral from P0 to P(ilow)
322      DO il = 1,ncum
323      IF (lswitch(il)) THEN
324        deltap = P0(il)-P(il,ilow(il))
325        dcin = RD*(BUOYz(il)+BUOY(il,ilow(il)))*deltap
326     $         /(P0(il)+P(il,ilow(il)))
327        CINB(il) = min(0.,dcin)
328      ENDIF
329      ENDDO
330c
331c
332C  Middle part of CINB : integral from P(ilow) to P(isublcl)
333cc      DO k = ilow,isublcl-1
334      DO k = 1,nl
335        DO il = 1,ncum
336        IF (lswitch(il)
337     $   .AND. k .GE. ilow(il) .AND. k .LE. isublcl(il)-1) THEN
338          deltap = P(il,k)-P(il,k+1)
339          dcin = 0.5*RD*(BUOY(il,k)+BUOY(il,k+1))*deltap/PH(il,k+1)
340          CINB(il) = CINB(il) + min(0.,dcin)
341        ENDIF
342        ENDDO
343      ENDDO
344c
345C  Upper part of CINB : integral from P(isublcl) to Plcl
346      DO il = 1,ncum
347      IF (lswitch(il)) THEN
348        deltap = P(il,isublcl(il)) - Plcl(il)
349        dcin = RD*(BUOY(il,isublcl(il))+BUOYlcl(il))*deltap
350     $         /(P(il,isublcl(il))+Plcl(il))
351        CINB(il) = CINB(il)+min(0.,dcin)
352      ENDIF
353      ENDDO
354C
355c
356cc      ENDIF
357c
358C 2.2 Computation of CINA
359c ---------------------
360c
361      DO il = 1,ncum
362        lswitch2(il)= Plcl(il) .GT. P(il,itop(il)-1)
363        lswitch(il) = lswitch1(il) .AND. lswitch2(il)
364      ENDDO
365c
366c 2.2.1 FIrst case : Plcl > P(itop-1)
367C ---------------------------------
368C In order to get Plfc, one may interpolate linearly buoyancies
369C  between P(itop) and P(itop-1).
370      DO il = 1,ncum
371      IF (lswitch(il)) THEN
372        PLFC(il)=(buoy(il,itop(il))*P(il,itop(il)-1)
373     $           -buoy(il,itop(il)-1)*P(il,itop(il)))
374     $           / (buoy(il,itop(il))-buoy(il,itop(il)-1))
375      ENDIF
376      ENDDO
377c
378C   Upper part of CINA : integral from P(itop-1) to Plfc
379      DO il = 1,ncum
380      IF (lswitch(il)) THEN
381        deltap = P(il,itop(il)-1)-Plfc(il)
382        dcin = RD*BUOY(il,itop(il)-1)*deltap
383     $         /(P(il,itop(il)-1)+Plfc(il))
384        CINA(il) = min(0.,dcin)
385      ENDIF
386      ENDDO
387c
388C   Middle part of CINA : integral from P(icb+1) to P(itop-1)
389      DO k = 1,nl
390        DO il = 1,ncum
391        IF (lswitch(il)
392     $     .AND. k .GE. icb(il)+1 .AND. k .LE. itop(il)-2) THEN
393          deltap = P(il,k)-P(il,k+1)
394          dcin = 0.5*RD*(BUOY(il,k)+BUOY(il,k+1))*deltap/PH(il,k+1)
395          CINA(il) = CINA(il) + min(0.,dcin)
396        ENDIF
397        ENDDO
398      ENDDO
399c
400C   Lower part of CINA : integral from Plcl to P(icb+1)
401      DO il = 1,ncum
402      IF (lswitch(il)) THEN
403        IF (Plcl(il) .GT. P(il,icb(il))) THEN
404          IF (icb(il) .LT. itop(il)-1) THEN
405            deltap = P(il,icb(il))-P(il,icb(il)+1)
406            dcin = 0.5*RD*(BUOY(il,icb(il))+BUOY(il,icb(il)+1))
407     $                   *deltap/PH(il,icb(il)+1)
408            CINA(il) = CINA(il)+min(0.,dcin)
409          ENDIF
410c
411          deltap = Plcl(il)-P(il,icb(il))
412          dcin = RD*(BUOYlcl(il)+BUOY(il,icb(il)))
413     $              *deltap/(Plcl(il)+P(il,icb(il)))
414          CINA(il) = CINA(il)+min(0.,dcin)
415        ELSE
416          deltap = Plcl(il)-P(il,icb(il)+1)
417          dcin = RD*(BUOYlcl(il)+BUOY(il,icb(il)+1))
418     $             *deltap/(Plcl(il)+P(il,icb(il)+1))
419          CINA(il) = CINA(il)+min(0.,dcin)
420        ENDIF
421      ENDIF
422      ENDDO
423c
424      DO il = 1,ncum
425        lswitch(il) = lswitch1(il) .AND. .NOT. lswitch2(il)
426      ENDDO
427cc      ELSE
428c
429c 2.2.2 Second case : Plcl lies between P(itop-1) and P(itop);
430C ----------------------------------------------------------
431C In order to get Plfc, one has to interpolate between P(itop) and Plcl.
432      DO il = 1,ncum
433      IF (lswitch(il)) THEN
434        PLFC(il) =
435     $    (BUOY(il,itop(il))*Plcl(il)-BUOYlcl(il)*P(il,itop(il)))
436     $          /(BUOY(il,itop(il))     -BUOYlcl(il))
437      ENDIF
438      ENDDO
439c
440      DO il = 1,ncum
441      IF (lswitch(il)) THEN
442        deltap = Plcl(il)-Plfc(il)
443        dcin = RD*BUOYlcl(il)*deltap/(Plcl(il)+Plfc(il))
444        CINA(il) = min(0.,dcin)
445      ENDIF
446      ENDDO
447cc      ENDIF
448c
449
450
451      RETURN
452      END
Note: See TracBrowser for help on using the repository browser.