My Project
 All Classes Files Functions Variables Macros
cv3_cine.F
Go to the documentation of this file.
1 !
2 ! $Id: cv3_cine.F 1403 2010-07-01 09:02:53Z fairhead $
3 !
4  SUBROUTINE cv3_cine(nloc,ncum,nd,icb,inb
5  : ,pbase,plcl,p,ph,tv,tvp
6  : ,cina,cinb,plfc)
7 
8 ***************************************************************
9 * *
10 * cv3_cine *
11 * *
12 * *
13 * written by : frederique cheruy *
14 * vectorization: jean-yves grandpeix, 19/06/2003, 11.54.43 *
15 * modified by : *
16 ***************************************************************
17 *
18  implicit none
19 c
20 #include "YOMCST.h"
21 #include "cvthermo.h"
22 #include "cv3param.h"
23 c 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 c
30 c output
31  real cina(nloc),cinb(nloc),plfc(nloc)
32 c
33 c 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)
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 c
45 c-------------------------------------------------------------
46 c initialization
47 c-------------------------------------------------------------
48  do il = 1,ncum
49  cina(il) = 0.
50  cinb(il) = 0.
51  enddo
52 c
53 c--------------------------------------------------------------
54 c recompute buoyancies
55 c--------------------------------------------------------------
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  ENDDO
61  ENDDO
62 c---------------------------------------------------------------
63 c
64 c calcul de la flottabilite a lcl(buoylcl)
65 c ifst = first p-level above lcl
66 c isublcl = highest p-level below lcl.
67 c---------------------------------------------------------------
68 c
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  enddo
72 c
73  do il = 1,ncum
74  IF (plcl(il) .GT. 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  ENDIF
81  enddo
82 c
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  enddo
87 c
88  do il = 1,ncum
89  buoylcl(il) = tvplcl(il)-tvlcl(il)
90  enddo
91 c
92 c---------------------------------------------------------------
93 c premiere couche contenant un niveau de flotabilite positive
94 c et premiere couche contenant un niveau de flotabilite negative
95 c au dessus du niveau de condensation
96 c---------------------------------------------------------------
97  do il = 1,ncum
98  itop(il) =nl-1
99  ineg(il) = nl-1
100  exist_lfc(il) = .false.
101  enddo
102  do 100 k=nl-1,1,-1
103  do 110 il=1,ncum
104  if (k .ge. ifst(il)) then
105  if (buoy(il,k) .gt. 0.) then
106  itop(il)=k
107  exist_lfc(il) = .true.
108  else
109  ineg(il)=k
110  endif
111  endif
112 110 continue
113 100 continue
114 c
115 c---------------------------------------------------------------
116 c when there is no positive buoyancy level, set plfc, cina and cinb
117 c to arbitrary extreme values.
118 c---------------------------------------------------------------
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  ENDIF
125  ENDDO
126 c
127 c
128 c---------------------------------------------------------------
129 c -- two cases : buoylcl >= 0 and buoylcl < 0.
130 c---------------------------------------------------------------
131 c
132 c--------------------
133 c -- 1.0 buoylcl >=0.
134 c--------------------
135 c
136  dpmax = 50.
137  DO il = 1,ncum
138  lswitch1(il)=buoylcl(il) .GE. 0. .AND. exist_lfc(il)
139  lswitch(il) = lswitch1(il)
140  ENDDO
141 c
142 c 1.1 no inhibition case
143 c ----------------------
144 c If buoyancy is positive at lcl and stays positive over a large enough
145 c pressure interval(=dpmax), inhibition is set to zero,
146 c
147  DO il = 1,ncum
148  IF (lswitch(il)) THEN
149  IF (p(il,ineg(il)) .LT. p(il,icb(il))-dpmax) THEN
150  plfc(il) = plcl(il)
151  cina(il) = 0.
152  cinb(il) = 0.
153  ENDIF
154  ENDIF
155  ENDDO
156 c
157 c 1.2 upper inhibition only case
158 c ------------------------------
159  DO il = 1,ncum
160  lswitch2(il)= p(il,ineg(il)) .GE. p(il,icb(il))-dpmax
161  lswitch(il) = lswitch1(il) .AND. lswitch2(il)
162  ENDDO
163 c
164  DO il = 1,ncum
165  IF (lswitch(il)) THEN
166  cinb(il) = 0.
167 c
168 c 1.2.1 calcul de la pression du niveau de flot. nulle juste au-dessus de lcl
169 c ---------------------------------------------------------------------------
170  IF (ineg(il) .GT. isublcl(il)+1) THEN
171 c in order to get p0, one may interpolate linearly buoyancies
172 c between p(ineg) and p(ineg-1).
173  p0(il)=(buoy(il,ineg(il))*p(il,ineg(il)-1)
174  $ -buoy(il,ineg(il)-1)*p(il,ineg(il)))
175  : / (buoy(il,ineg(il))-buoy(il,ineg(il)-1))
176  ELSE
177 c in order to get p0, one has to interpolate between p(ineg) and plcl.
178  p0(il) = (buoy(il,ineg(il))*plcl(il)-buoylcl(il)*p(il,ineg(il)))
179  $ /(buoy(il,ineg(il)) -buoylcl(il))
180  ENDIF
181  ENDIF
182  ENDDO
183 c
184 c 1.2.2 recompute itop(=1st layer with positive buoyancy above ineg)
185 c -------------------------------------------------------------------
186  do il = 1,ncum
187  IF (lswitch(il)) THEN
188  itop(il) =nl-1
189  ENDIF
190  enddo
191 c
192  do k=nl,1,-1
193  do il=1,ncum
194  IF (lswitch(il)) THEN
195  if (k .ge. ineg(il) .and. buoy(il,k) .gt. 0) then
196  itop(il)=k
197  endif
198  ENDIF
199  enddo
200  enddo
201 c
202 c 1.2.3 computation of plfc
203 c -------------------------
204  DO il = 1,ncum
205  IF (lswitch(il)) THEN
206  plfc(il)=(buoy(il,itop(il))*p(il,itop(il)-1)
207  $ -buoy(il,itop(il)-1)*p(il,itop(il)))
208  $ / (buoy(il,itop(il))-buoy(il,itop(il)-1))
209  ENDIF
210  ENDDO
211 c
212 c 1.2.4 computation of cina
213 c -------------------------
214 c
215 c upper part of cina : integral from p(itop-1) to plfc
216  DO il = 1,ncum
217  IF (lswitch(il)) THEN
218  deltap = p(il,itop(il)-1)-plfc(il)
219  dcin = rd*buoy(il,itop(il)-1)*deltap
220  $ / (p(il,itop(il)-1)+plfc(il))
221  cina(il) = min(0.,dcin)
222  ENDIF
223  ENDDO
224 c
225 c middle part of cina : integral from p(ineg) to p(itop-1)
226  DO k = 1,nl
227  DO il = 1,ncum
228  IF (lswitch(il)) THEN
229  IF (k .GE. ineg(il) .AND. k .LE. itop(il)-2) THEN
230  deltap = p(il,k)-p(il,k+1)
231  dcin = 0.5*rd*(buoy(il,k)+buoy(il,k+1))*deltap/ph(il,k+1)
232  cina(il) = cina(il) + min(0.,dcin)
233  ENDIF
234  ENDIF
235  ENDDO
236  ENDDO
237 c
238 c lower part of cina : integral from p0 to p(ineg)
239  DO il = 1,ncum
240  IF (lswitch(il)) THEN
241  deltap = p0(il)-p(il,ineg(il))
242  dcin = rd*buoy(il,ineg(il))*deltap/(p(il,ineg(il))+p0(il))
243  cina(il) = cina(il) + min(0.,dcin)
244  ENDIF
245  ENDDO
246 c
247 c
248 c ------------------
249 c -- 2.0 buoylcl <0.
250 c ------------------
251 c
252  DO il = 1,ncum
253  lswitch1(il)=buoylcl(il) .LT. 0. .AND. exist_lfc(il)
254  lswitch(il) = lswitch1(il)
255  ENDDO
256 c
257 c 2.0.1 premiere couche ou la flotabilite est negative au dessus du sol
258 c ----------------------------------------------------
259 c au cas ou elle existe sinon ilow=1 (nk apres)
260 c on suppose que la parcelle part de la premiere couche
261 c
262  DO il = 1,ncum
263  IF (lswitch(il)) THEN
264  ilow(il)=1
265  ENDIF
266  ENDDO
267 c
268  do 200 k=nl,1,-1
269  DO il = 1,ncum
270  IF (lswitch(il) .AND. k .LE.icb(il)-1) THEN
271  if(buoy(il,k).lt. 0.) then
272  ilow(il) = k
273  endif
274  ENDIF
275  ENDDO
276  200 continue
277 
278 c 2.0.2 calcul de la pression du niveau de flot. nulle sous le nuage
279 c ----------------------------------------------------
280  DO il = 1,ncum
281  IF (lswitch(il)) THEN
282  if(ilow(il).gt. 1) then
283  p0(il)=(buoy(il,ilow(il))*p(il,ilow(il)-1)
284  $ -buoy(il,ilow(il)-1)*p(il,ilow(il)))
285  : / (buoy(il,ilow(il))-buoy(il,ilow(il)-1))
286  buoyz(il) = 0.
287  else
288  p0(il) = p(il,1)
289  buoyz(il) = buoy(il,1)
290  endif
291  ENDIF
292  ENDDO
293 c
294 c 2.1. computation of cinb
295 c -----------------------
296 c
297  DO il = 1,ncum
298  lswitch2(il)= (isublcl(il) .EQ. 1 .AND. ilow(il) .EQ. 1)
299  $ .OR.(isublcl(il) .EQ. ilow(il)-1)
300  lswitch(il) = lswitch1(il) .AND. lswitch2(il)
301  ENDDO
302 cc IF ( (isublcl .EQ. 1 .AND. ilow .EQ. 1)
303 cc $ .OR.(isublcl .EQ. ilow-1)) THEN
304 c
305 c 2.1.1 first case : plcl just above p0
306 c -------------------------------------
307  DO il = 1,ncum
308  IF (lswitch(il)) THEN
309  deltap = p0(il)-plcl(il)
310  dcin = rd*(buoyz(il)+buoylcl(il))*deltap/(p0(il)+plcl(il))
311  cinb(il) = min(0.,dcin)
312  ENDIF
313  ENDDO
314 c
315  DO il = 1,ncum
316  lswitch(il) = lswitch1(il) .AND. .NOT. lswitch2(il)
317  ENDDO
318 cc ELSE
319 c
320 c 2.1.2 second case : there is at least one p-level between p0 and plcl
321 c ---------------------------------------------------------------------
322 c
323 c lower part of cinb : integral from p0 to p(ilow)
324  DO il = 1,ncum
325  IF (lswitch(il)) THEN
326  deltap = p0(il)-p(il,ilow(il))
327  dcin = rd*(buoyz(il)+buoy(il,ilow(il)))*deltap
328  $ /(p0(il)+p(il,ilow(il)))
329  cinb(il) = min(0.,dcin)
330  ENDIF
331  ENDDO
332 c
333 c
334 c middle part of cinb : integral from p(ilow) to p(isublcl)
335 cc DO k = ilow,isublcl-1
336  DO k = 1,nl
337  DO il = 1,ncum
338  IF (lswitch(il)
339  $ .AND. k .GE. ilow(il) .AND. k .LE. isublcl(il)-1) THEN
340  deltap = p(il,k)-p(il,k+1)
341  dcin = 0.5*rd*(buoy(il,k)+buoy(il,k+1))*deltap/ph(il,k+1)
342  cinb(il) = cinb(il) + min(0.,dcin)
343  ENDIF
344  ENDDO
345  ENDDO
346 c
347 c upper part of cinb : integral from p(isublcl) to plcl
348  DO il = 1,ncum
349  IF (lswitch(il)) THEN
350  deltap = p(il,isublcl(il)) - plcl(il)
351  dcin = rd*(buoy(il,isublcl(il))+buoylcl(il))*deltap
352  $ /(p(il,isublcl(il))+plcl(il))
353  cinb(il) = cinb(il)+min(0.,dcin)
354  ENDIF
355  ENDDO
356 c
357 c
358 cc ENDIF
359 c
360 c 2.2 computation of cina
361 c ---------------------
362 c
363  DO il = 1,ncum
364  lswitch2(il)= plcl(il) .GT. p(il,itop(il)-1)
365  lswitch(il) = lswitch1(il) .AND. lswitch2(il)
366  ENDDO
367 c
368 c 2.2.1 first case : plcl > p(itop-1)
369 c ---------------------------------
370 c in order to get plfc, one may interpolate linearly buoyancies
371 c between p(itop) and p(itop-1).
372  DO il = 1,ncum
373  IF (lswitch(il)) THEN
374  plfc(il)=(buoy(il,itop(il))*p(il,itop(il)-1)
375  $ -buoy(il,itop(il)-1)*p(il,itop(il)))
376  $ / (buoy(il,itop(il))-buoy(il,itop(il)-1))
377  ENDIF
378  ENDDO
379 c
380 c upper part of cina : integral from p(itop-1) to plfc
381  DO il = 1,ncum
382  IF (lswitch(il)) THEN
383  deltap = p(il,itop(il)-1)-plfc(il)
384  dcin = rd*buoy(il,itop(il)-1)*deltap
385  $ /(p(il,itop(il)-1)+plfc(il))
386  cina(il) = min(0.,dcin)
387  ENDIF
388  ENDDO
389 c
390 c middle part of cina : integral from p(icb+1) to p(itop-1)
391  DO k = 1,nl
392  DO il = 1,ncum
393  IF (lswitch(il)
394  $ .AND. k .GE. icb(il)+1 .AND. k .LE. itop(il)-2) THEN
395  deltap = p(il,k)-p(il,k+1)
396  dcin = 0.5*rd*(buoy(il,k)+buoy(il,k+1))*deltap/ph(il,k+1)
397  cina(il) = cina(il) + min(0.,dcin)
398  ENDIF
399  ENDDO
400  ENDDO
401 c
402 c lower part of cina : integral from plcl to p(icb+1)
403  DO il = 1,ncum
404  IF (lswitch(il)) THEN
405  IF (plcl(il) .GT. p(il,icb(il))) THEN
406  IF (icb(il) .LT. itop(il)-1) THEN
407  deltap = p(il,icb(il))-p(il,icb(il)+1)
408  dcin = 0.5*rd*(buoy(il,icb(il))+buoy(il,icb(il)+1))
409  $ *deltap/ph(il,icb(il)+1)
410  cina(il) = cina(il)+min(0.,dcin)
411  ENDIF
412 c
413  deltap = plcl(il)-p(il,icb(il))
414  dcin = rd*(buoylcl(il)+buoy(il,icb(il)))
415  $ *deltap/(plcl(il)+p(il,icb(il)))
416  cina(il) = cina(il)+min(0.,dcin)
417  ELSE
418  deltap = plcl(il)-p(il,icb(il)+1)
419  dcin = rd*(buoylcl(il)+buoy(il,icb(il)+1))
420  $ *deltap/(plcl(il)+p(il,icb(il)+1))
421  cina(il) = cina(il)+min(0.,dcin)
422  ENDIF
423  ENDIF
424  ENDDO
425 c
426  DO il = 1,ncum
427  lswitch(il) = lswitch1(il) .AND. .NOT. lswitch2(il)
428  ENDDO
429 cc ELSE
430 c
431 c 2.2.2 second case : plcl lies between p(itop-1) and p(itop);
432 c ----------------------------------------------------------
433 c in order to get plfc, one has to interpolate between p(itop) and plcl.
434  DO il = 1,ncum
435  IF (lswitch(il)) THEN
436  plfc(il) =
437  $ (buoy(il,itop(il))*plcl(il)-buoylcl(il)*p(il,itop(il)))
438  $ /(buoy(il,itop(il)) -buoylcl(il))
439  ENDIF
440  ENDDO
441 c
442  DO il = 1,ncum
443  IF (lswitch(il)) THEN
444  deltap = plcl(il)-plfc(il)
445  dcin = rd*buoylcl(il)*deltap/(plcl(il)+plfc(il))
446  cina(il) = min(0.,dcin)
447  ENDIF
448  ENDDO
449 cc ENDIF
450 c
451 
452 
453  RETURN
454  END