source: LMDZ5/trunk/libf/phylmd/cv3_cine.F90 @ 2237

Last change on this file since 2237 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

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