source: LMDZ4/trunk/libf/phylmd/cv3_cine.F @ 1372

Last change on this file since 1372 was 1146, checked in by Laurent Fairhead, 16 years ago

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF
  • 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)
[1146]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)
[1146]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--------------------------------------------------------------
[1146]54      DO k = 1,nd
[879]55        DO il = 1,ncum
[1146]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---------------------------------------------------------------
[1146]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
[1146]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
[1146]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 ----------------------------------------------------
[1146]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
[1146]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.