source: lmdz_wrf/trunk/WRFV3/lmdz/cv3_cine.F90 @ 1554

Last change on this file since 1554 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 14.3 KB
Line 
1!
2! $Id: cv3_cine.F 1403 2010-07-01 09:02:53Z fairhead $
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
19!c
20#include "YOMCST.h"
21#include "cvthermo.h"
22#include "cv3param.h"
23!c 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)
29!c
30!c output
31      real cina(nloc),cinb(nloc),plfc(nloc)
32!c
33!c 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)
44!c
45!c-------------------------------------------------------------
46!c     Initialization
47!c-------------------------------------------------------------
48      do il = 1,ncum
49       cina(il) = 0.
50       cinb(il) = 0.
51      enddo
52!c
53!c--------------------------------------------------------------
54!c      Recompute buoyancies
55!c--------------------------------------------------------------
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
62!c---------------------------------------------------------------
63!c
64!c   calcul de la flottabilite a LCL (Buoylcl)
65!c     ifst = first P-level above lcl
66!c     isublcl = highest P-level below lcl.
67!c---------------------------------------------------------------
68!c
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
72!c
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
82!c
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
87!c
88      do il = 1,ncum
89        BUOYlcl(il) = TVPlcl(il)-TVlcl(il)
90      enddo
91!c
92!c---------------------------------------------------------------
93!c premiere couche contenant un  niveau de flotabilite positive
94!c et premiere couche contenant un  niveau de flotabilite negative
95!c  au dessus du niveau de condensation
96!c---------------------------------------------------------------
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
114!c
115!c---------------------------------------------------------------
116!c When there is no positive buoyancy level, set Plfc, Cina and Cinb
117!c to arbitrary extreme values.
118!c---------------------------------------------------------------
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
126!c
127!c
128!c---------------------------------------------------------------
129!c -- Two cases : BUOYlcl >= 0 and BUOYlcl < 0.
130!c---------------------------------------------------------------
131!C
132!C--------------------
133!C -- 1.0 BUOYlcl >=0.
134!C--------------------
135!c
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
141!c
142!c 1.1 No inhibition case
143!c ----------------------
144!C   If buoyancy is positive at LCL and stays positive over a large enough
145!C pressure interval (=DPMAX), inhibition is set to zero,
146!C
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
156!c
157!c 1.2 Upper inhibition only case
158!c ------------------------------
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
163!c
164      DO il = 1,ncum
165      IF (lswitch(il)) THEN
166          Cinb(il) = 0.
167!c
168!c 1.2.1  Calcul de la pression du niveau de flot. nulle juste au-dessus de LCL
169!c ---------------------------------------------------------------------------
170         IF (ineg(il) .GT. isublcl(il)+1) THEN
171!C In order to get P0, one may interpolate linearly buoyancies
172!C  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
177!C 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
183!c
184!c 1.2.2 Recompute itop (=1st layer with positive buoyancy above ineg)
185!c -------------------------------------------------------------------
186      do il = 1,ncum
187      IF (lswitch(il)) THEN
188        itop(il) =nl-1
189      ENDIF
190      enddo
191!c
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
201!c
202!c 1.2.3 Computation of PLFC
203!c -------------------------
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
211!c
212!c 1.2.4 Computation of CINA
213!c -------------------------
214!c
215!C   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
224!c
225!C   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
237!c
238!C   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
246!c
247!C
248!C ------------------
249!C -- 2.0 BUOYlcl <0.
250!C ------------------
251!C
252      DO il = 1,ncum
253        lswitch1(il)=BUOYlcl(il) .LT. 0. .AND. exist_lfc(il)
254        lswitch(il) = lswitch1(il)
255      ENDDO
256!c
257!c 2.0.1 Premiere  couche ou la flotabilite est negative au dessus du sol
258!c ----------------------------------------------------
259!c    au cas ou elle existe  sinon ilow=1 (nk apres)
260!c      on suppose que la parcelle part de la premiere couche
261!c
262      DO il = 1,ncum
263      IF (lswitch(il)) THEN
264       ilow(il)=1
265      ENDIF
266      ENDDO
267!c
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
278!c 2.0.2  Calcul de la pression du niveau de flot. nulle sous le nuage
279!c ----------------------------------------------------
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
293!c
294!C 2.1. Computation of CINB
295!C -----------------------
296!c
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
302!cc      IF (    (isublcl .EQ. 1 .AND. ilow .EQ. 1)
303!cc     $    .OR.(isublcl .EQ. ilow-1)) THEN
304!c
305!c 2.1.1 First case : Plcl just above P0
306!c -------------------------------------
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
314!c
315      DO il = 1,ncum
316        lswitch(il) = lswitch1(il) .AND. .NOT. lswitch2(il)
317      ENDDO
318!cc      ELSE
319!c
320!c 2.1.2 Second case : there is at least one P-level between P0 and Plcl
321!c ---------------------------------------------------------------------
322!c
323!C   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
332!c
333!c
334!C  Middle part of CINB : integral from P(ilow) to P(isublcl)
335!cc      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
346!c
347!C  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
356!C
357!c
358!cc      ENDIF
359!c
360!C 2.2 Computation of CINA
361!c ---------------------
362!c
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
367!c
368!c 2.2.1 FIrst case : Plcl > P(itop-1)
369!C ---------------------------------
370!C In order to get Plfc, one may interpolate linearly buoyancies
371!C  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
379!c
380!C   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
389!c
390!C   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
401!c
402!C   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
412!c
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
425!c
426      DO il = 1,ncum
427        lswitch(il) = lswitch1(il) .AND. .NOT. lswitch2(il)
428      ENDDO
429!cc      ELSE
430!c
431!c 2.2.2 Second case : Plcl lies between P(itop-1) and P(itop);
432!C ----------------------------------------------------------
433!C 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
441!c
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
449!cc      ENDIF
450!c
451
452
453      RETURN
454      END SUBROUTINE cv3_cine
Note: See TracBrowser for help on using the repository browser.