LMDZ
cv3_cine.F90
Go to the documentation of this file.
1 
2 ! $Id: cv3_cine.F90 1992 2014-03-05 13:19:12Z lguez $
3 
4 SUBROUTINE 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
445 END SUBROUTINE cv3_cine
cv3_cine
subroutine cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, cinb, plfc)
Definition: cv3_cine.F90:6
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
nl
!$Id!Parameters for nl
Definition: cv30param.h:5
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139