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

Last change on this file since 5715 was 5692, checked in by yann meurdesoif, 10 days ago

Convection GPU porting : set convection subroutines into module

Files will be renamed later to *_mod.f90

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