Context Navigation

diagedyn.F @ 1178

Last change on this file since 1178 was 545, checked in by lmdzadmin, 20 years ago
Suite des modifs pour adaptation a Prism MED LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.5 KB

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

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

Download in other formats:

Original Format