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

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