Context Navigation

diagedyn.F @ 4785

Last change on this file since 4785 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.7 KB

Line
1	!
2	! $Id: diagedyn.F 1279 2009-12-10 09:02:56Z idelkadi $
3	!
4
5	C======================================================================
6	SUBROUTINE diagedyn(tit,iprt,idiag,idiag2,dtime
7	e , ucov , vcov , ps, p ,pk , teta , q, ql)
8	C======================================================================
9	C
10	C Purpose:
11	C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
12	C et calcul le flux de chaleur et le flux d'eau necessaire a ces
13	C changements. Ces valeurs sont moyennees sur la surface de tout
14	C le globe et sont exprime en W/2 et kg/s/m2
15	C Outil pour diagnostiquer la conservation de l'energie
16	C et de la masse dans la dynamique.
17	C
18	C
19	c======================================================================
20	C Arguments:
21	C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
22	C iprt----input-I- PRINT level ( <=1 : no PRINT)
23	C idiag---input-I- indice dans lequel sera range les nouveaux
24	C bilans d' entalpie et de masse
25	C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
26	C sont compare au bilan de d'enthalpie de masse de
27	C l'indice numero idiag2
28	C Cas parriculier : si idiag2=0, pas de comparaison, on
29	c sort directement les bilans d'enthalpie et de masse
30	C dtime----input-R- time step (s)
31	C uconv, vconv-input-R- vents covariants (m/s)
32	C ps-------input-R- Surface pressure (Pa)
33	C p--------input-R- pressure at the interfaces
34	C pk-------input-R- pk= (p/Pref)**kappa
35	c teta-----input-R- potential temperature (K)
36	c q--------input-R- vapeur d'eau (kg/kg)
37	c ql-------input-R- liquid watter (kg/kg)
38	c aire-----input-R- mesh surafce (m2)
39	c
40	C the following total value are computed by UNIT of earth surface
41	C
42	C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
43	c change (J/m2) during one time step (dtime) for the whole
44	C atmosphere (air, watter vapour, liquid and solid)
45	C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
46	C total watter (kg/m2) change during one time step (dtime),
47	C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
48	C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
49	C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
50	C
51	C
52	C J.L. Dufresne, July 2002
53	c======================================================================
54
55	IMPLICIT NONE
56	C
57	#include "dimensions.h"
58	#include "paramet.h"
59	#include "comgeom.h"
60	#include "iniprint.h"
61
62	#ifdef CPP_EARTH
63	#include "../phylmd/YOMCST.h"
64	#include "../phylmd/YOETHF.h"
65	#endif
66	C
67	INTEGER imjmp1
68	PARAMETER( imjmp1=iim*jjp1)
69	c Input variables
70	CHARACTER*15 tit
71	INTEGER iprt,idiag, idiag2
72	REAL dtime
73	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
74	REAL ps(ip1jmp1) ! pression au sol
75	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
76	REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa
77	REAL teta(ip1jmp1,llm) ! temperature potentielle
78	REAL q(ip1jmp1,llm) ! champs eau vapeur
79	REAL ql(ip1jmp1,llm) ! champs eau liquide
80
81
82	c Output variables
83	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
84	C
85	C Local variables
86	c
87	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
88	. , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
89	c h_vcol_tot-- total enthalpy of vertical air column
90	C (air with watter vapour, liquid and solid) (J/m2)
91	c h_dair_tot-- total enthalpy of dry air (J/m2)
92	c h_qw_tot---- total enthalpy of watter vapour (J/m2)
93	c h_ql_tot---- total enthalpy of liquid watter (J/m2)
94	c h_qs_tot---- total enthalpy of solid watter (J/m2)
95	c qw_tot------ total mass of watter vapour (kg/m2)
96	c ql_tot------ total mass of liquid watter (kg/m2)
97	c qs_tot------ total mass of solid watter (kg/m2)
98	c ec_tot------ total cinetic energy (kg/m2)
99	C
100	REAL masse(ip1jmp1,llm) ! masse d'air
101	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
102	REAL ecin(ip1jmp1,llm)
103
104	REAL zaire(imjmp1)
105	REAL zps(imjmp1)
106	REAL zairm(imjmp1,llm)
107	REAL zecin(imjmp1,llm)
108	REAL zpaprs(imjmp1,llm)
109	REAL zpk(imjmp1,llm)
110	REAL zt(imjmp1,llm)
111	REAL zh(imjmp1,llm)
112	REAL zqw(imjmp1,llm)
113	REAL zql(imjmp1,llm)
114	REAL zqs(imjmp1,llm)
115
116	REAL zqw_col(imjmp1)
117	REAL zql_col(imjmp1)
118	REAL zqs_col(imjmp1)
119	REAL zec_col(imjmp1)
120	REAL zh_dair_col(imjmp1)
121	REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
122	C
123	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
124	C
125	REAL airetot, zcpvap, zcwat, zcice
126	C
127	INTEGER i, k, jj, ij , l ,ip1jjm1
128	C
129	INTEGER ndiag ! max number of diagnostic in parallel
130	PARAMETER (ndiag=10)
131	integer pas(ndiag)
132	save pas
133	data pas/ndiag*0/
134	C
135	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
136	$ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
137	$ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
138	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
139	$ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
140
141
142	#ifdef CPP_EARTH
143	c======================================================================
144	C Compute Kinetic enrgy
145	CALL covcont ( llm , ucov , vcov , ucont, vcont )
146	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
147	CALL massdair( p, masse )
148	c======================================================================
149	C
150	C
151	print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
152	return
153	C On ne garde les donnees que dans les colonnes i=1,iim
154	DO jj = 1,jjp1
155	ip1jjm1=iip1*(jj-1)
156	DO ij = 1,iim
157	i=iim*(jj-1)+ij
158	zaire(i)=aire(ij+ip1jjm1)
159	zps(i)=ps(ij+ip1jjm1)
160	ENDDO
161	ENDDO
162	C 3D arrays
163	DO l = 1, llm
164	DO jj = 1,jjp1
165	ip1jjm1=iip1*(jj-1)
166	DO ij = 1,iim
167	i=iim*(jj-1)+ij
168	zairm(i,l) = masse(ij+ip1jjm1,l)
169	zecin(i,l) = ecin(ij+ip1jjm1,l)
170	zpaprs(i,l) = p(ij+ip1jjm1,l)
171	zpk(i,l) = pk(ij+ip1jjm1,l)
172	zh(i,l) = teta(ij+ip1jjm1,l)
173	zqw(i,l) = q(ij+ip1jjm1,l)
174	zql(i,l) = ql(ij+ip1jjm1,l)
175	zqs(i,l) = 0.
176	ENDDO
177	ENDDO
178	ENDDO
179	C
180	C Reset variables
181	DO i = 1, imjmp1
182	zqw_col(i)=0.
183	zql_col(i)=0.
184	zqs_col(i)=0.
185	zec_col(i) = 0.
186	zh_dair_col(i) = 0.
187	zh_qw_col(i) = 0.
188	zh_ql_col(i) = 0.
189	zh_qs_col(i) = 0.
190	ENDDO
191	C
192	zcpvap=RCPV
193	zcwat=RCW
194	zcice=RCS
195	C
196	C Compute vertical sum for each atmospheric column
197	C ================================================
198	DO k = 1, llm
199	DO i = 1, imjmp1
200	C Watter mass
201	zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
202	zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
203	zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
204	C Cinetic Energy
205	zec_col(i) = zec_col(i)
206	$ +zecin(i,k)*zairm(i,k)
207	C Air enthalpy
208	zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
209	zh_dair_col(i) = zh_dair_col(i)
210	$ + RCPD(1.-zqw(i,k)-zql(i,k)-zqs(i,k))zairm(i,k)*zt(i,k)
211	zh_qw_col(i) = zh_qw_col(i)
212	$ + zcpvapzqw(i,k)zairm(i,k)*zt(i,k)
213	zh_ql_col(i) = zh_ql_col(i)
214	$ + zcwatzql(i,k)zairm(i,k)*zt(i,k)
215	$ - RLVTTzql(i,k)zairm(i,k)
216	zh_qs_col(i) = zh_qs_col(i)
217	$ + zcicezqs(i,k)zairm(i,k)*zt(i,k)
218	$ - RLSTTzqs(i,k)zairm(i,k)
219
220	END DO
221	ENDDO
222	C
223	C Mean over the planete surface
224	C =============================
225	qw_tot = 0.
226	ql_tot = 0.
227	qs_tot = 0.
228	ec_tot = 0.
229	h_vcol_tot = 0.
230	h_dair_tot = 0.
231	h_qw_tot = 0.
232	h_ql_tot = 0.
233	h_qs_tot = 0.
234	airetot=0.
235	C
236	do i=1,imjmp1
237	qw_tot = qw_tot + zqw_col(i)
238	ql_tot = ql_tot + zql_col(i)
239	qs_tot = qs_tot + zqs_col(i)
240	ec_tot = ec_tot + zec_col(i)
241	h_dair_tot = h_dair_tot + zh_dair_col(i)
242	h_qw_tot = h_qw_tot + zh_qw_col(i)
243	h_ql_tot = h_ql_tot + zh_ql_col(i)
244	h_qs_tot = h_qs_tot + zh_qs_col(i)
245	airetot=airetot+zaire(i)
246	END DO
247	C
248	qw_tot = qw_tot/airetot
249	ql_tot = ql_tot/airetot
250	qs_tot = qs_tot/airetot
251	ec_tot = ec_tot/airetot
252	h_dair_tot = h_dair_tot/airetot
253	h_qw_tot = h_qw_tot/airetot
254	h_ql_tot = h_ql_tot/airetot
255	h_qs_tot = h_qs_tot/airetot
256	C
257	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
258	C
259	C Compute the change of the atmospheric state compare to the one
260	C stored in "idiag2", and convert it in flux. THis computation
261	C is performed IF idiag2 /= 0 and IF it is not the first CALL
262	c for "idiag"
263	C ===================================
264	C
265	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
266	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
267	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
268	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
269	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
270	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
271	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
272	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
273	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
274	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
275	d_qt = d_qw + d_ql + d_qs
276	ELSE
277	d_h_vcol = 0.
278	d_h_dair = 0.
279	d_h_qw = 0.
280	d_h_ql = 0.
281	d_h_qs = 0.
282	d_qw = 0.
283	d_ql = 0.
284	d_qs = 0.
285	d_ec = 0.
286	d_qt = 0.
287	ENDIF
288	C
289	IF (iprt.ge.2) THEN
290	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
291	9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
292	$ ,1i6,10(1pE14.6))
293	WRITE(6,9001) tit,pas(idiag), d_h_vcol
294	9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
295	WRITE(6,9002) tit,pas(idiag), d_ec
296	9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
297	C WRITE(6,9003) tit,pas(idiag), ec_tot
298	9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
299	WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
300	9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
301	END IF
302	C
303	C Store the new atmospheric state in "idiag"
304	C
305	pas(idiag)=pas(idiag)+1
306	h_vcol_pre(idiag) = h_vcol_tot
307	h_dair_pre(idiag) = h_dair_tot
308	h_qw_pre(idiag) = h_qw_tot
309	h_ql_pre(idiag) = h_ql_tot
310	h_qs_pre(idiag) = h_qs_tot
311	qw_pre(idiag) = qw_tot
312	ql_pre(idiag) = ql_tot
313	qs_pre(idiag) = qs_tot
314	ec_pre (idiag) = ec_tot
315	C
316	#else
317	write(lunout,*)'diagedyn: Needs Earth physics to function'
318	#endif
319	! #endif of #ifdef CPP_EARTH
320	RETURN
321	END

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

Download in other formats:

Original Format