Context Navigation

cv3_cine.F @ 2440

Last change on this file since 2440 was 1322, checked in by Laurent Fairhead, 15 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 13.3 KB

Line
1	!
2	! $Id: cv3_cine.F 1322 2010-03-12 10:54:11Z lguez $
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 = RDBUOY(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.5RD(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 = RDBUOY(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.5RD(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 = RDBUOY(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.5RD(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.5RD(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 = RDBUOYlcl(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

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

Download in other formats:

Original Format