source: LMDZ5/branches/LF-private/libf/phylmd/cv3_cine.F @ 5456

Last change on this file since 5456 was 1861, checked in by lguez, 11 years ago

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