Context Navigation

diagedyn.F @ 7

Last change on this file since 7 was 5, checked in by slebonnois, 14 years ago

Ajout des fichiers .def Venus et Titan (tels qu'ils sont utilisés

actuellement) dans les deftanks.

Ajout d'une doc sur Cp(T).
Modifications dans dyn3d concernant Cp(T), cf le log (v5) dans chantiers
Premières modifs de l'appel à la physique dans dyn3d/calfis, cf log (v5)
Elimination des cpdet.* dans phytitan et phyvenus (remplacés par cpdet.F

dans dyn3d).

File size: 10.8 KB

Line
1	!
2	! $Id: diagedyn.F 1279 2009-12-10 09:02:56Z fairhead $
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	! ADAPTATION GCM POUR CP(T)
199	call tpot2t(imjmp1*llm,zh,zt,zpk)
200	DO k = 1, llm
201	DO i = 1, imjmp1
202	C Watter mass
203	zqw_col(i) = zqw_col(i) + zqw(i,k)*zairm(i,k)
204	zql_col(i) = zql_col(i) + zql(i,k)*zairm(i,k)
205	zqs_col(i) = zqs_col(i) + zqs(i,k)*zairm(i,k)
206	C Cinetic Energy
207	zec_col(i) = zec_col(i)
208	$ +zecin(i,k)*zairm(i,k)
209	C Air enthalpy
210	zh_dair_col(i) = zh_dair_col(i)
211	< ! ADAPTATION GCM POUR CP(T)
212	$ + cpdet(zt(i,k))*(1.-zqw(i,k)-zql(i,k)-zqs(i,k))
213	$ zairm(i,k)zt(i,k)
214	zh_qw_col(i) = zh_qw_col(i)
215	$ + zcpvapzqw(i,k)zairm(i,k)*zt(i,k)
216	zh_ql_col(i) = zh_ql_col(i)
217	$ + zcwatzql(i,k)zairm(i,k)*zt(i,k)
218	$ - RLVTTzql(i,k)zairm(i,k)
219	zh_qs_col(i) = zh_qs_col(i)
220	$ + zcicezqs(i,k)zairm(i,k)*zt(i,k)
221	$ - RLSTTzqs(i,k)zairm(i,k)
222
223	END DO
224	ENDDO
225	C
226	C Mean over the planete surface
227	C =============================
228	qw_tot = 0.
229	ql_tot = 0.
230	qs_tot = 0.
231	ec_tot = 0.
232	h_vcol_tot = 0.
233	h_dair_tot = 0.
234	h_qw_tot = 0.
235	h_ql_tot = 0.
236	h_qs_tot = 0.
237	airetot=0.
238	C
239	do i=1,imjmp1
240	qw_tot = qw_tot + zqw_col(i)
241	ql_tot = ql_tot + zql_col(i)
242	qs_tot = qs_tot + zqs_col(i)
243	ec_tot = ec_tot + zec_col(i)
244	h_dair_tot = h_dair_tot + zh_dair_col(i)
245	h_qw_tot = h_qw_tot + zh_qw_col(i)
246	h_ql_tot = h_ql_tot + zh_ql_col(i)
247	h_qs_tot = h_qs_tot + zh_qs_col(i)
248	airetot=airetot+zaire(i)
249	END DO
250	C
251	qw_tot = qw_tot/airetot
252	ql_tot = ql_tot/airetot
253	qs_tot = qs_tot/airetot
254	ec_tot = ec_tot/airetot
255	h_dair_tot = h_dair_tot/airetot
256	h_qw_tot = h_qw_tot/airetot
257	h_ql_tot = h_ql_tot/airetot
258	h_qs_tot = h_qs_tot/airetot
259	C
260	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
261	C
262	C Compute the change of the atmospheric state compare to the one
263	C stored in "idiag2", and convert it in flux. THis computation
264	C is performed IF idiag2 /= 0 and IF it is not the first CALL
265	c for "idiag"
266	C ===================================
267	C
268	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
269	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
270	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
271	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
272	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
273	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
274	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
275	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
276	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
277	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
278	d_qt = d_qw + d_ql + d_qs
279	ELSE
280	d_h_vcol = 0.
281	d_h_dair = 0.
282	d_h_qw = 0.
283	d_h_ql = 0.
284	d_h_qs = 0.
285	d_qw = 0.
286	d_ql = 0.
287	d_qs = 0.
288	d_ec = 0.
289	d_qt = 0.
290	ENDIF
291	C
292	IF (iprt.ge.2) THEN
293	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
294	9000 format('Dyn3d. Watter Mass Budget (kg/m2/s)',A15
295	$ ,1i6,10(1pE14.6))
296	WRITE(6,9001) tit,pas(idiag), d_h_vcol
297	9001 format('Dyn3d. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
298	WRITE(6,9002) tit,pas(idiag), d_ec
299	9002 format('Dyn3d. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
300	C WRITE(6,9003) tit,pas(idiag), ec_tot
301	9003 format('Dyn3d. Cinetic Energy (W/m2) ',A15,1i6,10(E15.6))
302	WRITE(6,9004) tit,pas(idiag), d_h_vcol+d_ec
303	9004 format('Dyn3d. Total Energy Budget (W/m2) ',A15,1i6,10(F8.2))
304	END IF
305	C
306	C Store the new atmospheric state in "idiag"
307	C
308	pas(idiag)=pas(idiag)+1
309	h_vcol_pre(idiag) = h_vcol_tot
310	h_dair_pre(idiag) = h_dair_tot
311	h_qw_pre(idiag) = h_qw_tot
312	h_ql_pre(idiag) = h_ql_tot
313	h_qs_pre(idiag) = h_qs_tot
314	qw_pre(idiag) = qw_tot
315	ql_pre(idiag) = ql_tot
316	qs_pre(idiag) = qs_tot
317	ec_pre (idiag) = ec_tot
318	C
319	#else
320	write(lunout,*)'diagedyn: Needs Earth physics to function'
321	#endif
322	! #endif of #ifdef CPP_EARTH
323	RETURN
324	END

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

Download in other formats:

Original Format