source: LMDZ6/trunk/libf/phylmd/cv3_cine.f90 @ 5489

Last change on this file since 5489 was 5346, checked in by fhourdin, 12 months ago

Debut de replaysation de la convection profonde.

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