source: LMDZ6/branches/Amaury_dev/libf/phylmd/cv3_cine.F90 @ 5218

Last change on this file since 5218 was 5144, checked in by abarral, 16 months ago

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