Context Navigation

cv3_cine.F @ 5429

Last change on this file since 5429 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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 1403 2010-07-01 09:02:53Z fhourdin $
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