Context Navigation

diagedyn.F @ 517

Last change on this file since 517 was 431, checked in by lmdzadmin, 22 years ago

Convergence avec la version de Ionela Musat

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

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

Download in other formats:

Original Format