Context Navigation

diagedyn.F @ 839

Last change on this file since 839 was 774, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.6 KB

Line
1	!
2	! $Header$
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
61	#ifdef CPP_PHYS
62	#include "../phylmd/YOMCST.h"
63	#include "../phylmd/YOETHF.h"
64	#endif
65	C
66	INTEGER imjmp1
67	PARAMETER( imjmp1=iim*jjp1)
68	c Input variables
69	CHARACTER*15 tit
70	INTEGER iprt,idiag, idiag2
71	REAL dtime
72	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
73	REAL ps(ip1jmp1) ! pression au sol
74	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
75	REAL pk (ip1jmp1,llm ) ! = (p/Pref)**kappa
76	REAL teta(ip1jmp1,llm) ! temperature potentielle
77	REAL q(ip1jmp1,llm) ! champs eau vapeur
78	REAL ql(ip1jmp1,llm) ! champs eau liquide
79
80
81	c Output variables
82	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
83	C
84	C Local variables
85	c
86	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
87	. , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
88	c h_vcol_tot-- total enthalpy of vertical air column
89	C (air with watter vapour, liquid and solid) (J/m2)
90	c h_dair_tot-- total enthalpy of dry air (J/m2)
91	c h_qw_tot---- total enthalpy of watter vapour (J/m2)
92	c h_ql_tot---- total enthalpy of liquid watter (J/m2)
93	c h_qs_tot---- total enthalpy of solid watter (J/m2)
94	c qw_tot------ total mass of watter vapour (kg/m2)
95	c ql_tot------ total mass of liquid watter (kg/m2)
96	c qs_tot------ total mass of solid watter (kg/m2)
97	c ec_tot------ total cinetic energy (kg/m2)
98	C
99	REAL masse(ip1jmp1,llm) ! masse d'air
100	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
101	REAL ecin(ip1jmp1,llm)
102
103	REAL zaire(imjmp1)
104	REAL zps(imjmp1)
105	REAL zairm(imjmp1,llm)
106	REAL zecin(imjmp1,llm)
107	REAL zpaprs(imjmp1,llm)
108	REAL zpk(imjmp1,llm)
109	REAL zt(imjmp1,llm)
110	REAL zh(imjmp1,llm)
111	REAL zqw(imjmp1,llm)
112	REAL zql(imjmp1,llm)
113	REAL zqs(imjmp1,llm)
114
115	REAL zqw_col(imjmp1)
116	REAL zql_col(imjmp1)
117	REAL zqs_col(imjmp1)
118	REAL zec_col(imjmp1)
119	REAL zh_dair_col(imjmp1)
120	REAL zh_qw_col(imjmp1), zh_ql_col(imjmp1), zh_qs_col(imjmp1)
121	C
122	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
123	C
124	REAL airetot, zcpvap, zcwat, zcice
125	C
126	INTEGER i, k, jj, ij , l ,ip1jjm1
127	C
128	INTEGER ndiag ! max number of diagnostic in parallel
129	PARAMETER (ndiag=10)
130	integer pas(ndiag)
131	save pas
132	data pas/ndiag*0/
133	C
134	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
135	$ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
136	$ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
137	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
138	$ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
139
140
141	#ifdef CPP_PHYS
142	c======================================================================
143	C Compute Kinetic enrgy
144	CALL covcont ( llm , ucov , vcov , ucont, vcont )
145	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
146	CALL massdair( p, masse )
147	c======================================================================
148	C
149	C
150	print*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
151	return
152	C On ne garde les donnees que dans les colonnes i=1,iim
153	DO jj = 1,jjp1
154	ip1jjm1=iip1*(jj-1)
155	DO ij = 1,iim
156	i=iim*(jj-1)+ij
157	zaire(i)=aire(ij+ip1jjm1)
158	zps(i)=ps(ij+ip1jjm1)
159	ENDDO
160	ENDDO
161	C 3D arrays
162	DO l = 1, llm
163	DO jj = 1,jjp1
164	ip1jjm1=iip1*(jj-1)
165	DO ij = 1,iim
166	i=iim*(jj-1)+ij
167	zairm(i,l) = masse(ij+ip1jjm1,l)
168	zecin(i,l) = ecin(ij+ip1jjm1,l)
169	zpaprs(i,l) = p(ij+ip1jjm1,l)
170	zpk(i,l) = pk(ij+ip1jjm1,l)
171	zh(i,l) = teta(ij+ip1jjm1,l)
172	zqw(i,l) = q(ij+ip1jjm1,l)
173	zql(i,l) = ql(ij+ip1jjm1,l)
174	zqs(i,l) = 0.
175	ENDDO
176	ENDDO
177	ENDDO
178	C
179	C Reset variables
180	DO i = 1, imjmp1
181	zqw_col(i)=0.
182	zql_col(i)=0.
183	zqs_col(i)=0.
184	zec_col(i) = 0.
185	zh_dair_col(i) = 0.
186	zh_qw_col(i) = 0.
187	zh_ql_col(i) = 0.
188	zh_qs_col(i) = 0.
189	ENDDO
190	C
191	zcpvap=RCPV
192	zcwat=RCW
193	zcice=RCS
194	C
195	C Compute vertical sum for each atmospheric column
196	C ================================================
197	DO k = 1, llm
198	DO i = 1, imjmp1
199	C Watter mass
200	zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
201	zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
202	zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
203	C Cinetic Energy
204	zec_col(i) = zec_col(i)
205	$ +zecin(i,k)*zairm(i,k)
206	C Air enthalpy
207	zt(i,k)= zh(i,k) * zpk(i,k) / RCPD
208	zh_dair_col(i) = zh_dair_col(i)
209	$ + RCPD(1.-zqw(i,k)-zql(i,k)-zqs(i,k))zairm(i,k)*zt(i,k)
210	zh_qw_col(i) = zh_qw_col(i)
211	$ + zcpvapzqw(i,k)zairm(i,k)*zt(i,k)
212	zh_ql_col(i) = zh_ql_col(i)
213	$ + zcwatzql(i,k)zairm(i,k)*zt(i,k)
214	$ - RLVTTzql(i,k)zairm(i,k)
215	zh_qs_col(i) = zh_qs_col(i)
216	$ + zcicezqs(i,k)zairm(i,k)*zt(i,k)
217	$ - RLSTTzqs(i,k)zairm(i,k)
218
219	END DO
220	ENDDO
221	C
222	C Mean over the planete surface
223	C =============================
224	qw_tot = 0.
225	ql_tot = 0.
226	qs_tot = 0.
227	ec_tot = 0.
228	h_vcol_tot = 0.
229	h_dair_tot = 0.
230	h_qw_tot = 0.
231	h_ql_tot = 0.
232	h_qs_tot = 0.
233	airetot=0.
234	C
235	do i=1,imjmp1
236	qw_tot = qw_tot + zqw_col(i)
237	ql_tot = ql_tot + zql_col(i)
238	qs_tot = qs_tot + zqs_col(i)
239	ec_tot = ec_tot + zec_col(i)
240	h_dair_tot = h_dair_tot + zh_dair_col(i)
241	h_qw_tot = h_qw_tot + zh_qw_col(i)
242	h_ql_tot = h_ql_tot + zh_ql_col(i)
243	h_qs_tot = h_qs_tot + zh_qs_col(i)
244	airetot=airetot+zaire(i)
245	END DO
246	C
247	qw_tot = qw_tot/airetot
248	ql_tot = ql_tot/airetot
249	qs_tot = qs_tot/airetot
250	ec_tot = ec_tot/airetot
251	h_dair_tot = h_dair_tot/airetot
252	h_qw_tot = h_qw_tot/airetot
253	h_ql_tot = h_ql_tot/airetot
254	h_qs_tot = h_qs_tot/airetot
255	C
256	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
257	C
258	C Compute the change of the atmospheric state compare to the one
259	C stored in "idiag2", and convert it in flux. THis computation
260	C is performed IF idiag2 /= 0 and IF it is not the first CALL
261	c for "idiag"
262	C ===================================
263	C
264	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
265	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
266	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
267	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
268	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
269	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
270	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
271	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
272	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
273	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
274	d_qt = d_qw + d_ql + d_qs
275	ELSE
276	d_h_vcol = 0.
277	d_h_dair = 0.
278	d_h_qw = 0.
279	d_h_ql = 0.
280	d_h_qs = 0.
281	d_qw = 0.
282	d_ql = 0.
283	d_qs = 0.
284	d_ec = 0.
285	d_qt = 0.
286	ENDIF
287	C
288	IF (iprt.ge.2) THEN
289	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
290	9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
291	$ ,1i6,10(1pE14.6))
292	WRITE(6,9001) tit,pas(idiag), d_h_vcol
293	9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
294	WRITE(6,9002) tit,pas(idiag), d_ec
295	9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
296	C WRITE(6,9003) tit,pas(idiag), ec_tot
297	9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
298	WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
299	9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
300	END IF
301	C
302	C Store the new atmospheric state in "idiag"
303	C
304	pas(idiag)=pas(idiag)+1
305	h_vcol_pre(idiag) = h_vcol_tot
306	h_dair_pre(idiag) = h_dair_tot
307	h_qw_pre(idiag) = h_qw_tot
308	h_ql_pre(idiag) = h_ql_tot
309	h_qs_pre(idiag) = h_qs_tot
310	qw_pre(idiag) = qw_tot
311	ql_pre(idiag) = ql_tot
312	qs_pre(idiag) = qs_tot
313	ec_pre (idiag) = ec_tot
314	C
315	#else
316	print*,'Pour l instant diagedyn a besoin de la physique'
317	#endif
318	RETURN
319	END

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

Download in other formats:

Original Format