Context Navigation

cv3_cine.F @ 1105

Last change on this file since 1105 was 879, checked in by Laurent Fairhead, 16 years ago
Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ... LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format