Context Navigation

diagphy.F @ 1062

Last change on this file since 1062 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.1 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE diagphy(airephy,tit,iprt
5	$ , tops, topl, sols, soll, sens
6	$ , evap, rain_fall, snow_fall, ts
7	$ , d_etp_tot, d_qt_tot, d_ec_tot
8	$ , fs_bound, fq_bound)
9	C======================================================================
10	C
11	C Purpose:
12	C Compute the thermal flux and the watter mass flux at the atmosphere
13	c boundaries. Print them and also the atmospheric enthalpy change and
14	C the atmospheric mass change.
15	C
16	C Arguments:
17	C airephy-------input-R- grid area
18	C tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
19	C iprt--------input-I- PRINT level ( <=0 : no PRINT)
20	C tops(klon)--input-R- SW rad. at TOA (W/m2), positive up.
21	C topl(klon)--input-R- LW rad. at TOA (W/m2), positive down
22	C sols(klon)--input-R- Net SW flux above surface (W/m2), positive up
23	C (i.e. -1 * flux absorbed by the surface)
24	C soll(klon)--input-R- Net LW flux above surface (W/m2), positive up
25	C (i.e. flux emited - flux absorbed by the surface)
26	C sens(klon)--input-R- Sensible Flux at surface (W/m2), positive down
27	C evap(klon)--input-R- Evaporation + sublimation watter vapour mass flux
28	C (kg/m2/s), positive up
29	C rain_fall(klon)
30	C --input-R- Liquid watter mass flux (kg/m2/s), positive down
31	C snow_fall(klon)
32	C --input-R- Solid watter mass flux (kg/m2/s), positive down
33	C ts(klon)----input-R- Surface temperature (K)
34	C d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
35	C change (W/m2)
36	C d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
37	C change (kg/m2/s)
38	C d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
39	C change (W/m2)
40	C
41	C fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
42	C fq_bound---output-R- Watter mass flux at the atmosphere boundaries (kg/m2/s)
43	C
44	C J.L. Dufresne, July 2002
45	C======================================================================
46	C
47	implicit none
48
49	#include "dimensions.h"
50	#include "dimphy.h"
51	#include "YOMCST.h"
52	#include "YOETHF.h"
53	C
54	C Input variables
55	real airephy(klon)
56	CHARACTER*15 tit
57	INTEGER iprt
58	real tops(klon),topl(klon),sols(klon),soll(klon)
59	real sens(klon),evap(klon),rain_fall(klon),snow_fall(klon)
60	REAL ts(klon)
61	REAL d_etp_tot, d_qt_tot, d_ec_tot
62	c Output variables
63	REAL fs_bound, fq_bound
64	C
65	C Local variables
66	real stops,stopl,ssols,ssoll
67	real ssens,sfront,slat
68	real airetot, zcpvap, zcwat, zcice
69	REAL rain_fall_tot, snow_fall_tot, evap_tot
70	C
71	integer i
72	C
73	integer pas
74	save pas
75	data pas/0/
76	c$OMP THREADPRIVATE(pas)
77	C
78	pas=pas+1
79	stops=0.
80	stopl=0.
81	ssols=0.
82	ssoll=0.
83	ssens=0.
84	sfront = 0.
85	evap_tot = 0.
86	rain_fall_tot = 0.
87	snow_fall_tot = 0.
88	airetot=0.
89	C
90	C Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
91	C la glace, on travaille par difference a la chaleur specifique de l'
92	c air sec. En effet, comme on travaille a niveau de pression donne,
93	C toute variation de la masse d'un constituant est totalement
94	c compense par une variation de masse d'air.
95	C
96	zcpvap=RCPV-RCPD
97	zcwat=RCW-RCPD
98	zcice=RCS-RCPD
99	C
100	do i=1,klon
101	stops=stops+tops(i)*airephy(i)
102	stopl=stopl+topl(i)*airephy(i)
103	ssols=ssols+sols(i)*airephy(i)
104	ssoll=ssoll+soll(i)*airephy(i)
105	ssens=ssens+sens(i)*airephy(i)
106	sfront = sfront
107	$ + ( evap(i)zcpvap-rain_fall(i)zcwat-snow_fall(i)*zcice
108	$ ) ts(i) airephy(i)
109	evap_tot = evap_tot + evap(i)*airephy(i)
110	rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
111	snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
112	airetot=airetot+airephy(i)
113	enddo
114	stops=stops/airetot
115	stopl=stopl/airetot
116	ssols=ssols/airetot
117	ssoll=ssoll/airetot
118	ssens=ssens/airetot
119	sfront = sfront/airetot
120	evap_tot = evap_tot /airetot
121	rain_fall_tot = rain_fall_tot/airetot
122	snow_fall_tot = snow_fall_tot/airetot
123	C
124	slat = RLVTT * rain_fall_tot + RLSTT * snow_fall_tot
125	C Heat flux at atm. boundaries
126	fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront
127	$ + slat
128	C Watter flux at atm. boundaries
129	fq_bound = evap_tot - rain_fall_tot -snow_fall_tot
130	C
131	IF (iprt.ge.1) write(6,6666)
132	$ tit, pas, fs_bound, d_etp_tot, fq_bound, d_qt_tot
133	C
134	IF (iprt.ge.1) write(6,6668)
135	$ tit, pas, d_etp_tot+d_ec_tot-fs_bound, d_qt_tot-fq_bound
136	C
137	IF (iprt.ge.2) write(6,6667)
138	$ tit, pas, stops,stopl,ssols,ssoll,ssens,slat,evap_tot
139	$ ,rain_fall_tot+snow_fall_tot
140
141	return
142
143	6666 format('Phys. Flux Budget ',a15,1i6,2f8.2,2(1pE13.5))
144	6667 format('Phys. Boundary Flux ',a15,1i6,6f8.2,2(1pE13.5))
145	6668 format('Phys. Total Budget ',a15,1i6,f8.2,2(1pE13.5))
146
147	end
148
149	C======================================================================
150	SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime
151	e ,t,q,ql,qs,u,v,paprs,pplay
152	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
153	C======================================================================
154	C
155	C Purpose:
156	C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
157	C et calcul le flux de chaleur et le flux d'eau necessaire a ces
158	C changements. Ces valeurs sont moyennees sur la surface de tout
159	C le globe et sont exprime en W/2 et kg/s/m2
160	C Outil pour diagnostiquer la conservation de l'energie
161	C et de la masse dans la physique. Suppose que les niveau de
162	c pression entre couche ne varie pas entre 2 appels.
163	C
164	C Plusieurs de ces diagnostics peuvent etre fait en parallele: les
165	c bilans sont sauvegardes dans des tableaux indices. On parlera
166	C "d'indice de diagnostic"
167	c
168	C
169	c======================================================================
170	C Arguments:
171	C airephy-------input-R- grid area
172	C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
173	C iprt----input-I- PRINT level ( <=1 : no PRINT)
174	C idiag---input-I- indice dans lequel sera range les nouveaux
175	C bilans d' entalpie et de masse
176	C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
177	C sont compare au bilan de d'enthalpie de masse de
178	C l'indice numero idiag2
179	C Cas parriculier : si idiag2=0, pas de comparaison, on
180	c sort directement les bilans d'enthalpie et de masse
181	C dtime----input-R- time step (s)
182	c t--------input-R- temperature (K)
183	c q--------input-R- vapeur d'eau (kg/kg)
184	c ql-------input-R- liquid watter (kg/kg)
185	c qs-------input-R- solid watter (kg/kg)
186	c u--------input-R- vitesse u
187	c v--------input-R- vitesse v
188	c paprs----input-R- pression a intercouche (Pa)
189	c pplay----input-R- pression au milieu de couche (Pa)
190	c
191	C the following total value are computed by UNIT of earth surface
192	C
193	C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
194	c change (J/m2) during one time step (dtime) for the whole
195	C atmosphere (air, watter vapour, liquid and solid)
196	C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
197	C total watter (kg/m2) change during one time step (dtime),
198	C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
199	C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
200	C d_qs------output-R- same, for the solid watter only (kg/m2/s)
201	C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
202	C
203	C other (COMMON...)
204	C RCPD, RCPV, ....
205	C
206	C J.L. Dufresne, July 2002
207	c======================================================================
208
209	IMPLICIT NONE
210	C
211	#include "dimensions.h"
212	#include "dimphy.h"
213	#include "YOMCST.h"
214	#include "YOETHF.h"
215	C
216	c Input variables
217	real airephy(klon)
218	CHARACTER*15 tit
219	INTEGER iprt,idiag, idiag2
220	REAL dtime
221	REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
222	REAL u(klon,klev), v(klon,klev)
223	REAL paprs(klon,klev+1), pplay(klon,klev)
224	c Output variables
225	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
226	C
227	C Local variables
228	c
229	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
230	. , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
231	c h_vcol_tot-- total enthalpy of vertical air column
232	C (air with watter vapour, liquid and solid) (J/m2)
233	c h_dair_tot-- total enthalpy of dry air (J/m2)
234	c h_qw_tot---- total enthalpy of watter vapour (J/m2)
235	c h_ql_tot---- total enthalpy of liquid watter (J/m2)
236	c h_qs_tot---- total enthalpy of solid watter (J/m2)
237	c qw_tot------ total mass of watter vapour (kg/m2)
238	c ql_tot------ total mass of liquid watter (kg/m2)
239	c qs_tot------ total mass of solid watter (kg/m2)
240	c ec_tot------ total cinetic energy (kg/m2)
241	C
242	REAL zairm(klon,klev) ! layer air mass (kg/m2)
243	REAL zqw_col(klon)
244	REAL zql_col(klon)
245	REAL zqs_col(klon)
246	REAL zec_col(klon)
247	REAL zh_dair_col(klon)
248	REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
249	C
250	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
251	C
252	REAL airetot, zcpvap, zcwat, zcice
253	C
254	INTEGER i, k
255	C
256	INTEGER ndiag ! max number of diagnostic in parallel
257	PARAMETER (ndiag=10)
258	integer pas(ndiag)
259	save pas
260	data pas/ndiag*0/
261	c$OMP THREADPRIVATE(pas)
262	C
263	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
264	$ , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
265	$ , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
266	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
267	$ , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
268	c$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
269	c$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
270	c======================================================================
271	C
272	DO k = 1, klev
273	DO i = 1, klon
274	C layer air mass
275	zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
276	ENDDO
277	END DO
278	C
279	C Reset variables
280	DO i = 1, klon
281	zqw_col(i)=0.
282	zql_col(i)=0.
283	zqs_col(i)=0.
284	zec_col(i) = 0.
285	zh_dair_col(i) = 0.
286	zh_qw_col(i) = 0.
287	zh_ql_col(i) = 0.
288	zh_qs_col(i) = 0.
289	ENDDO
290	C
291	zcpvap=RCPV
292	zcwat=RCW
293	zcice=RCS
294	C
295	C Compute vertical sum for each atmospheric column
296	C ================================================
297	DO k = 1, klev
298	DO i = 1, klon
299	C Watter mass
300	zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
301	zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
302	zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
303	C Cinetic Energy
304	zec_col(i) = zec_col(i)
305	$ +0.5(u(i,k)2+v(i,k)2)zairm(i,k)
306	C Air enthalpy
307	zh_dair_col(i) = zh_dair_col(i)
308	$ + RCPD(1.-q(i,k)-ql(i,k)-qs(i,k))zairm(i,k)*t(i,k)
309	zh_qw_col(i) = zh_qw_col(i)
310	$ + zcpvapq(i,k)zairm(i,k)*t(i,k)
311	zh_ql_col(i) = zh_ql_col(i)
312	$ + zcwatql(i,k)zairm(i,k)*t(i,k)
313	$ - RLVTTql(i,k)zairm(i,k)
314	zh_qs_col(i) = zh_qs_col(i)
315	$ + zciceqs(i,k)zairm(i,k)*t(i,k)
316	$ - RLSTTqs(i,k)zairm(i,k)
317
318	END DO
319	ENDDO
320	C
321	C Mean over the planete surface
322	C =============================
323	qw_tot = 0.
324	ql_tot = 0.
325	qs_tot = 0.
326	ec_tot = 0.
327	h_vcol_tot = 0.
328	h_dair_tot = 0.
329	h_qw_tot = 0.
330	h_ql_tot = 0.
331	h_qs_tot = 0.
332	airetot=0.
333	C
334	do i=1,klon
335	qw_tot = qw_tot + zqw_col(i)*airephy(i)
336	ql_tot = ql_tot + zql_col(i)*airephy(i)
337	qs_tot = qs_tot + zqs_col(i)*airephy(i)
338	ec_tot = ec_tot + zec_col(i)*airephy(i)
339	h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
340	h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
341	h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
342	h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
343	airetot=airetot+airephy(i)
344	END DO
345	C
346	qw_tot = qw_tot/airetot
347	ql_tot = ql_tot/airetot
348	qs_tot = qs_tot/airetot
349	ec_tot = ec_tot/airetot
350	h_dair_tot = h_dair_tot/airetot
351	h_qw_tot = h_qw_tot/airetot
352	h_ql_tot = h_ql_tot/airetot
353	h_qs_tot = h_qs_tot/airetot
354	C
355	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
356	C
357	C Compute the change of the atmospheric state compare to the one
358	C stored in "idiag2", and convert it in flux. THis computation
359	C is performed IF idiag2 /= 0 and IF it is not the first CALL
360	c for "idiag"
361	C ===================================
362	C
363	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
364	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
365	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
366	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
367	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
368	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
369	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
370	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
371	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
372	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
373	d_qt = d_qw + d_ql + d_qs
374	ELSE
375	d_h_vcol = 0.
376	d_h_dair = 0.
377	d_h_qw = 0.
378	d_h_ql = 0.
379	d_h_qs = 0.
380	d_qw = 0.
381	d_ql = 0.
382	d_qs = 0.
383	d_ec = 0.
384	d_qt = 0.
385	ENDIF
386	C
387	IF (iprt.ge.2) THEN
388	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
389	9000 format('Phys. Watter Mass Budget (kg/m2/s)',A15
390	$ ,1i6,10(1pE14.6))
391	WRITE(6,9001) tit,pas(idiag), d_h_vcol
392	9001 format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
393	WRITE(6,9002) tit,pas(idiag), d_ec
394	9002 format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
395	END IF
396	C
397	C Store the new atmospheric state in "idiag"
398	C
399	pas(idiag)=pas(idiag)+1
400	h_vcol_pre(idiag) = h_vcol_tot
401	h_dair_pre(idiag) = h_dair_tot
402	h_qw_pre(idiag) = h_qw_tot
403	h_ql_pre(idiag) = h_ql_tot
404	h_qs_pre(idiag) = h_qs_tot
405	qw_pre(idiag) = qw_tot
406	ql_pre(idiag) = ql_tot
407	qs_pre(idiag) = qs_tot
408	ec_pre (idiag) = ec_tot
409	C
410	RETURN
411	END

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

Download in other formats:

Original Format