Context Navigation

diagphy_mod.F90 @ 38

Last change on this file since 38 was 38, checked in by lfita, 10 years ago
Fixing string concatenation... This is Fortran not python!!!!
File size: 21.1 KB

Line
1	!
2	! $Header$
3	!
4	MODULE diagphy_mod
5
6	CONTAINS
7
8	SUBROUTINE diagphy(airephy,tit,iprt &
9	& , tops, topl, sols, soll, sens &
10	& , evap, rain_fall, snow_fall, ts &
11	& , d_etp_tot, d_qt_tot, d_ec_tot &
12	& , fs_bound, fq_bound)
13	!C======================================================================
14	!C
15	!C Purpose:
16	!C Compute the thermal flux and the watter mass flux at the atmosphere
17	!c boundaries. Print them and also the atmospheric enthalpy change and
18	!C the atmospheric mass change.
19	!C
20	!C Arguments:
21	!C airephy-------input-R- grid area
22	!C tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
23	!C iprt--------input-I- PRINT level ( <=0 : no PRINT)
24	!C tops(klon)--input-R- SW rad. at TOA (W/m2), positive up.
25	!C topl(klon)--input-R- LW rad. at TOA (W/m2), positive down
26	!C sols(klon)--input-R- Net SW flux above surface (W/m2), positive up
27	!C (i.e. -1 * flux absorbed by the surface)
28	!C soll(klon)--input-R- Net LW flux above surface (W/m2), positive up
29	!C (i.e. flux emited - flux absorbed by the surface)
30	!C sens(klon)--input-R- Sensible Flux at surface (W/m2), positive down
31	!C evap(klon)--input-R- Evaporation + sublimation watter vapour mass flux
32	!C (kg/m2/s), positive up
33	!C rain_fall(klon)
34	!C --input-R- Liquid watter mass flux (kg/m2/s), positive down
35	!C snow_fall(klon)
36	!C --input-R- Solid watter mass flux (kg/m2/s), positive down
37	!C ts(klon)----input-R- Surface temperature (K)
38	!C d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
39	!C change (W/m2)
40	!C d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
41	!C change (kg/m2/s)
42	!C d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
43	!C change (W/m2)
44	!C
45	!C fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
46	!C fq_bound---output-R- Watter mass flux at the atmosphere boundaries (kg/m2/s)
47	!C
48	!C J.L. Dufresne, July 2002
49	!C Version prise sur ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd
50	!C le 25 Novembre 2002.
51	!C======================================================================
52	!C
53	use dimphy
54	implicit none
55
56	#include "dimensions.h"
57	!ccccc#include "dimphy.h"
58	#include "YOMCST.h"
59	#include "YOETHF.h"
60	!C
61	!C Input variables
62	real airephy(klon)
63	CHARACTER*15 tit
64	INTEGER iprt
65	real tops(klon),topl(klon),sols(klon),soll(klon)
66	real sens(klon),evap(klon),rain_fall(klon),snow_fall(klon)
67	REAL ts(klon)
68	REAL d_etp_tot, d_qt_tot, d_ec_tot
69	!c Output variables
70	REAL fs_bound, fq_bound
71	!C
72	!C Local variables
73	real stops,stopl,ssols,ssoll
74	real ssens,sfront,slat
75	real airetot, zcpvap, zcwat, zcice
76	REAL rain_fall_tot, snow_fall_tot, evap_tot
77	!C
78	integer i
79	!C
80	integer pas
81
82	save pas
83	data pas/0/
84	!$OMP THREADPRIVATE(pas)
85	!C
86
87	! L. Fita, LMD July 2014
88	CHARACTER(LEN=50) :: errmsg, fname, varname
89	LOGICAL :: found
90	REAL :: largest
91
92	fname = 'diagphy'
93	errmsg = 'ERROR -- error -- ERROR -- error'
94	largest = 10000.
95
96	pas=pas+1
97	stops=0.
98	stopl=0.
99	ssols=0.
100	ssoll=0.
101	ssens=0.
102	sfront = 0.
103	evap_tot = 0.
104	rain_fall_tot = 0.
105	snow_fall_tot = 0.
106	airetot=0.
107	!C
108	!C Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
109	!C la glace, on travaille par difference a la chaleur specifique de l'
110	!c air sec. En effet, comme on travaille a niveau de pression donne,
111	!C toute variation de la masse d'un constituant est totalement
112	!c compense par une variation de masse d'air.
113	!C
114	zcpvap=RCPV-RCPD
115	zcwat=RCW-RCPD
116	zcice=RCS-RCPD
117	!C
118	do i=1,klon
119	stops=stops+tops(i)*airephy(i)
120	stopl=stopl+topl(i)*airephy(i)
121	ssols=ssols+sols(i)*airephy(i)
122	ssoll=ssoll+soll(i)*airephy(i)
123	ssens=ssens+sens(i)*airephy(i)
124	sfront = sfront &
125	& + ( evap(i)zcpvap-rain_fall(i)zcwat-snow_fall(i)*zcice &
126	& ) ts(i) airephy(i)
127	evap_tot = evap_tot + evap(i)*airephy(i)
128	rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
129	snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
130	airetot=airetot+airephy(i)
131	enddo
132	stops=stops/airetot
133	stopl=stopl/airetot
134	ssols=ssols/airetot
135	ssoll=ssoll/airetot
136	ssens=ssens/airetot
137	sfront = sfront/airetot
138	evap_tot = evap_tot /airetot
139	rain_fall_tot = rain_fall_tot/airetot
140	snow_fall_tot = snow_fall_tot/airetot
141	!C
142	slat = RLVTT * rain_fall_tot + RLSTT * snow_fall_tot
143	!C Heat flux at atm. boundaries
144	fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront &
145	& + slat
146	!C Watter flux at atm. boundaries
147	fq_bound = evap_tot - rain_fall_tot -snow_fall_tot
148	!C
149	IF (iprt.ge.1) write(6,6666) &
150	& tit, pas, fs_bound, d_etp_tot, fq_bound, d_qt_tot
151	!C
152	IF (iprt.ge.1) write(6,6668) &
153	& tit, pas, d_etp_tot+d_ec_tot-fs_bound, d_qt_tot-fq_bound
154	!C
155	IF (iprt.ge.2) write(6,6667) &
156	& tit, pas, stops,stopl,ssols,ssoll,ssens,slat,evap_tot &
157	& ,rain_fall_tot+snow_fall_tot
158
159	! L. Fita, LMD July 2014. Checking for consistency
160	IF (fs_bound .NE. fs_bound .OR. ABS(fs_bound) > largest) THEN
161	PRINT *,TRIM(errmsg)
162	PRINT *,' ' // TRIM(fname) // ': Wrong fs_bound= ',fs_bound,' !!!'
163	PRINT *,' fs_bound: Heat flux at atm. boundaries'
164	PRINT *,' fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront'
165	PRINT *,' airetot= ',airetot
166	IF (airetot .NE. airetot .OR. ABS(airetot) > largest) THEN
167	varname = 'airephy'
168	CALL check_var(fname, varname, airephy, klon, largest, .FALSE.)
169	END IF
170	PRINT *,' stops= ',stops
171	IF (stops .NE. stops .OR. ABS(stops) > largest) THEN
172	varname = 'tops'
173	CALL check_var(fname, varname, tops, klon, largest, .FALSE.)
174	END IF
175	PRINT *,' stopl= ',stopl
176	IF (stopl .NE. stopl .OR. ABS(stopl) > largest) THEN
177	varname = 'topl'
178	CALL check_var(fname, varname, topl, klon, largest, .FALSE.)
179	END IF
180	PRINT *,' ssols= ',ssols
181	IF (ssols .NE. ssols .OR. ABS(ssols) > largest) THEN
182	varname = 'sols'
183	CALL check_var(fname, varname, sols, klon, largest, .FALSE.)
184	END IF
185	PRINT *,' ssoll= ',ssoll
186	IF (ssoll .NE. ssoll .OR. ABS(ssoll) > largest) THEN
187	varname = 'soll'
188	CALL check_var(fname, varname, soll, klon, largest, .FALSE.)
189	END IF
190	PRINT *,' ssens= ',ssens
191	IF (ssens .NE. ssens .OR. ABS(ssens) > largest) THEN
192	varname = 'sens'
193	CALL check_var(fname, varname, sens, klon, largest, .FALSE.)
194	END IF
195	PRINT *,' sfront= ',sfront
196	IF (sfront .NE. sfront .OR. ABS(sfront) > largest) THEN
197	varname = 'evap'
198	CALL check_var(fname, varname, evap, klon, largest, .FALSE.)
199	varname = 'rain_fall'
200	CALL check_var(fname, varname, rain_fall, klon, largest, .FALSE.)
201	varname = 'snow_fall'
202	CALL check_var(fname, varname, snow_fall, klon, largest, .FALSE.)
203	varname = 'ts'
204	CALL check_var(fname, varname, ts, klon, largest, .FALSE.)
205	END IF
206	STOP
207	END IF
208
209	IF (fq_bound .NE. fq_bound .OR. ABS(fq_bound) > 100000.) THEN
210	PRINT *,TRIM(errmsg)
211	PRINT *,' ' // TRIM(fname) // ': Wrong fq_bound= ',fs_bound,' !!!'
212	PRINT *,' fq_bound: Watter flux at atm. boundaries'
213	PRINT *,' fq_bound = evap_tot - rain_fall_tot -snow_fall_tot'
214	PRINT *,' evap_tot= ',evap_tot
215	PRINT *,' airetot= ',airetot
216	IF (airetot .NE. airetot .OR. ABS(airetot) > largest) THEN
217	varname = 'airephy'
218	CALL check_var(fname, varname, airephy, klon, largest, .FALSE.)
219	END IF
220	IF (evap_tot .NE. evap_tot .OR. ABS(evap_tot) > largest) THEN
221	varname = 'evap'
222	CALL check_var(fname, varname, evap, klon, largest, .FALSE.)
223	END IF
224	PRINT *,' rain_fall_tot= ',rain_fall_tot
225	IF (rain_fall_tot .NE. rain_fall_tot .OR. ABS(rain_fall_tot) > largest) THEN
226	varname = 'rain_fall'
227	CALL check_var(fname, varname, rain_fall, klon, largest, .FALSE.)
228	END IF
229	PRINT *,' snow_fall_tot= ',snow_fall_tot
230	IF (snow_fall_tot .NE. snow_fall_tot .OR. ABS(snow_fall_tot) > largest) THEN
231	varname = 'snow_fall'
232	CALL check_var(fname, varname, snow_fall, klon, largest, .FALSE.)
233	END IF
234	STOP
235	END IF
236
237	return
238
239	6666 format('Phys. Flux Budget ',a15,1i6,2f8.2,2(1pE13.5))
240	6667 format('Phys. Boundary Flux ',a15,1i6,6f8.2,2(1pE13.5))
241	6668 format('Phys. Total Budget ',a15,1i6,f8.2,2(1pE13.5))
242
243	end SUBROUTINE diagphy
244
245	!L. Fita, LMD 2004. Check variable
246	SUBROUTINE check_var(funcn, varn, var, sizev, bigvalue, stoprun)
247	! Subroutine to check the consistency of a variable
248	! * NaN value: by definition is variable /= variable
249	! * bigvalue: allowd threshold for variable
250
251	IMPLICIT NONE
252
253	INTEGER, INTENT(IN) :: sizev
254	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
255	REAL, DIMENSION(sizev), INTENT(IN) :: var
256	REAL, INTENT(IN) :: bigvalue
257	LOGICAL, INTENT(IN) :: stoprun
258
259	! Local
260	INTEGER :: i, wrongi
261	CHARACTER(LEN=50) :: errmsg
262	LOGICAL :: found
263	REAL, DIMENSION(sizev) :: wrongvalues
264	INTEGER, DIMENSION(sizev) :: wronggridpt
265
266	!!!!!!! Variables
267	! funcn: at which functino of part of the program variable is checked
268	! varn: name of the variable
269	! var: variable to check
270	! sizev: sixe of the variable
271	! bigvalue: biggest attenaible value for the variable
272	! stoprun: Should the run stop if it founds a problem?
273
274	errmsg = 'ERROR -- error -- ERROR -- error'
275
276	found = .FALSE.
277	wrongi = 0
278	DO i=1,sizev
279	IF (var(i) /= var(i) .OR. ABS(var(i)) > bigvalue ) THEN
280	IF (wrongi == 0) found = .TRUE.
281	wrongi = wrongi + 1
282	wrongvalues(wrongi) = var(i)
283	wronggridpt(wrongi) = i
284	END IF
285	END DO
286
287	IF (found) THEN
288	PRINT *,TRIM(errmsg)
289	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
290	"' is wrong at i value___"
291	DO i=1,wrongi
292	PRINT *,wronggridpt(i), wrongvalues(i)
293	END DO
294	IF (stoprun) THEN
295	STOP
296	END IF
297	END IF
298
299	RETURN
300
301	END SUBROUTINE check_var
302
303	!C======================================================================
304	SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime &
305	& ,t,q,ql,qs,u,v,paprs,pplay &
306	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
307	!C======================================================================
308	!C
309	!C Purpose:
310	!C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
311	!C et calcul le flux de chaleur et le flux d'eau necessaire a ces
312	!C changements. Ces valeurs sont moyennees sur la surface de tout
313	!C le globe et sont exprime en W/2 et kg/s/m2
314	!C Outil pour diagnostiquer la conservation de l'energie
315	!C et de la masse dans la physique. Suppose que les niveau de
316	!c pression entre couche ne varie pas entre 2 appels.
317	!C
318	!C Plusieurs de ces diagnostics peuvent etre fait en parallele: les
319	!c bilans sont sauvegardes dans des tableaux indices. On parlera
320	!C "d'indice de diagnostic"
321	!c
322	!C
323	!c======================================================================
324	!C Arguments:
325	!C airephy-------input-R- grid area
326	!C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
327	!C iprt----input-I- PRINT level ( <=1 : no PRINT)
328	!C idiag---input-I- indice dans lequel sera range les nouveaux
329	!C bilans d' entalpie et de masse
330	!C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
331	!C sont compare au bilan de d'enthalpie de masse de
332	!C l'indice numero idiag2
333	!C Cas parriculier : si idiag2=0, pas de comparaison, on
334	!c sort directement les bilans d'enthalpie et de masse
335	!C dtime----input-R- time step (s)
336	!c t--------input-R- temperature (K)
337	!c q--------input-R- vapeur d'eau (kg/kg)
338	!c ql-------input-R- liquid watter (kg/kg)
339	!c qs-------input-R- solid watter (kg/kg)
340	!c u--------input-R- vitesse u
341	!c v--------input-R- vitesse v
342	!c paprs----input-R- pression a intercouche (Pa)
343	!c pplay----input-R- pression au milieu de couche (Pa)
344	!c
345	!C the following total value are computed by UNIT of earth surface
346	!C
347	!C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
348	!c change (J/m2) during one time step (dtime) for the whole
349	!C atmosphere (air, watter vapour, liquid and solid)
350	!C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
351	!C total watter (kg/m2) change during one time step (dtime),
352	!C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
353	!C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
354	!C d_qs------output-R- same, for the solid watter only (kg/m2/s)
355	!C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
356	!C
357	!C other (COMMON...)
358	!C RCPD, RCPV, ....
359	!C
360	!C J.L. Dufresne, July 2002
361	!c======================================================================
362
363	USE dimphy
364	IMPLICIT NONE
365	!C
366	#include "dimensions.h"
367	!cccccc#include "dimphy.h"
368	#include "YOMCST.h"
369	#include "YOETHF.h"
370	!C
371	!c Input variables
372	real airephy(klon)
373	CHARACTER*15 tit
374	INTEGER iprt,idiag, idiag2
375	REAL dtime
376	REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
377	REAL u(klon,klev), v(klon,klev)
378	REAL paprs(klon,klev+1), pplay(klon,klev)
379	!c Output variables
380	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
381	!C
382	!C Local variables
383	!c
384	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
385	& , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
386	!c h_vcol_tot-- total enthalpy of vertical air column
387	!C (air with watter vapour, liquid and solid) (J/m2)
388	!c h_dair_tot-- total enthalpy of dry air (J/m2)
389	!c h_qw_tot---- total enthalpy of watter vapour (J/m2)
390	!c h_ql_tot---- total enthalpy of liquid watter (J/m2)
391	!c h_qs_tot---- total enthalpy of solid watter (J/m2)
392	!c qw_tot------ total mass of watter vapour (kg/m2)
393	!c ql_tot------ total mass of liquid watter (kg/m2)
394	!c qs_tot------ total mass of solid watter (kg/m2)
395	!c ec_tot------ total cinetic energy (kg/m2)
396	!C
397	REAL zairm(klon,klev) ! layer air mass (kg/m2)
398	REAL zqw_col(klon)
399	REAL zql_col(klon)
400	REAL zqs_col(klon)
401	REAL zec_col(klon)
402	REAL zh_dair_col(klon)
403	REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
404	!C
405	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
406	!C
407	REAL airetot, zcpvap, zcwat, zcice
408	!C
409	INTEGER i, k
410	!C
411	INTEGER ndiag ! max number of diagnostic in parallel
412	PARAMETER (ndiag=10)
413	integer pas(ndiag)
414	save pas
415	data pas/ndiag*0/
416	!$OMP THREADPRIVATE(pas)
417	!C
418	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag) &
419	& , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag) &
420	& , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
421	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre &
422	& , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
423	!$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
424	!$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
425	!c======================================================================
426	!C
427	DO k = 1, klev
428	DO i = 1, klon
429	!C layer air mass
430	zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
431	ENDDO
432	END DO
433	!C
434	!C Reset variables
435	DO i = 1, klon
436	zqw_col(i)=0.
437	zql_col(i)=0.
438	zqs_col(i)=0.
439	zec_col(i) = 0.
440	zh_dair_col(i) = 0.
441	zh_qw_col(i) = 0.
442	zh_ql_col(i) = 0.
443	zh_qs_col(i) = 0.
444	ENDDO
445	!C
446	zcpvap=RCPV
447	zcwat=RCW
448	zcice=RCS
449	!C
450	!C Compute vertical sum for each atmospheric column
451	!C ================================================
452	DO k = 1, klev
453	DO i = 1, klon
454	!C Watter mass
455	zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
456	zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
457	zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
458	!C Cinetic Energy
459	zec_col(i) = zec_col(i) &
460	& +0.5(u(i,k)2+v(i,k)2)zairm(i,k)
461	!C Air enthalpy
462	zh_dair_col(i) = zh_dair_col(i) &
463	& + RCPD(1.-q(i,k)-ql(i,k)-qs(i,k))zairm(i,k)*t(i,k)
464	zh_qw_col(i) = zh_qw_col(i) &
465	& + zcpvapq(i,k)zairm(i,k)*t(i,k)
466	zh_ql_col(i) = zh_ql_col(i) &
467	& + zcwatql(i,k)zairm(i,k)*t(i,k) &
468	& - RLVTTql(i,k)zairm(i,k)
469	zh_qs_col(i) = zh_qs_col(i) &
470	& + zciceqs(i,k)zairm(i,k)*t(i,k) &
471	& - RLSTTqs(i,k)zairm(i,k)
472
473	END DO
474	ENDDO
475
476	!C
477	!C Mean over the planete surface
478	!C =============================
479	qw_tot = 0.
480	ql_tot = 0.
481	qs_tot = 0.
482	ec_tot = 0.
483	h_vcol_tot = 0.
484	h_dair_tot = 0.
485	h_qw_tot = 0.
486	h_ql_tot = 0.
487	h_qs_tot = 0.
488	airetot=0.
489	!C
490	do i=1,klon
491	qw_tot = qw_tot + zqw_col(i)*airephy(i)
492	ql_tot = ql_tot + zql_col(i)*airephy(i)
493	qs_tot = qs_tot + zqs_col(i)*airephy(i)
494	ec_tot = ec_tot + zec_col(i)*airephy(i)
495	h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
496	h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
497	h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
498	h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
499	airetot=airetot+airephy(i)
500	END DO
501	!C
502	qw_tot = qw_tot/airetot
503	ql_tot = ql_tot/airetot
504	qs_tot = qs_tot/airetot
505	ec_tot = ec_tot/airetot
506	h_dair_tot = h_dair_tot/airetot
507	h_qw_tot = h_qw_tot/airetot
508	h_ql_tot = h_ql_tot/airetot
509	h_qs_tot = h_qs_tot/airetot
510	!C
511	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
512	!C
513	!C Compute the change of the atmospheric state compare to the one
514	!C stored in "idiag2", and convert it in flux. THis computation
515	!C is performed IF idiag2 /= 0 and IF it is not the first CALL
516	!c for "idiag"
517	!C ===================================
518	!C
519	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
520	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
521	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
522	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
523	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
524	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
525	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
526	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
527	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
528	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
529	d_qt = d_qw + d_ql + d_qs
530
531	ELSE
532	d_h_vcol = 0.
533	d_h_dair = 0.
534	d_h_qw = 0.
535	d_h_ql = 0.
536	d_h_qs = 0.
537	d_qw = 0.
538	d_ql = 0.
539	d_qs = 0.
540	d_ec = 0.
541	d_qt = 0.
542	ENDIF
543	!C
544	IF (iprt.ge.2) THEN
545	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
546	9000 format('Phys. Watter Mass Budget (kg/m2/s)',A15 &
547	& ,1i6,10(1pE14.6))
548	WRITE(6,9001) tit,pas(idiag), d_h_vcol
549	9001 format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
550	WRITE(6,9002) tit,pas(idiag), d_ec
551	9002 format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
552	END IF
553	!C
554	!C Store the new atmospheric state in "idiag"
555	!C
556	pas(idiag)=pas(idiag)+1
557	h_vcol_pre(idiag) = h_vcol_tot
558	h_dair_pre(idiag) = h_dair_tot
559	h_qw_pre(idiag) = h_qw_tot
560	h_ql_pre(idiag) = h_ql_tot
561	h_qs_pre(idiag) = h_qs_tot
562	qw_pre(idiag) = qw_tot
563	ql_pre(idiag) = ql_tot
564	qs_pre(idiag) = qs_tot
565	ec_pre (idiag) = ec_tot
566	!C
567	RETURN
568	END SUBROUTINE diagetpq
569
570	END MODULE diagphy_mod
571

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

Download in other formats:

Original Format