Context Navigation

diagphy_mod.F90 @ 354

Last change on this file since 354 was 240, checked in by lfita, 10 years ago
Adding map factor on tendencies
File size: 30.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 = 10.e4
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*10.e16) THEN
167	varname = 'airephy'
168	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .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 (d_etp_tot .NE. d_etp_tot .OR. ABS(d_etp_tot) > largest) THEN
210	PRINT *,TRIM(errmsg)
211	PRINT *,' ' // TRIM(fname) // ': Wrong d_etp_tot= ',d_etp_tot,' !!!'
212	PRINT *,' d_etp_tot: Heat flux equivalent to atmospheric enthalpy' // &
213	' change (W/m2)'
214	PRINT *,' d_etp_tot = input value !'
215	STOP
216	END IF
217
218	IF (fq_bound .NE. fq_bound .OR. ABS(fq_bound) > largest) THEN
219	PRINT *,TRIM(errmsg)
220	PRINT *,' ' // TRIM(fname) // ': Wrong fq_bound= ',fs_bound,' !!!'
221	PRINT *,' fq_bound: Watter flux at atm. boundaries'
222	PRINT *,' fq_bound = evap_tot - rain_fall_tot -snow_fall_tot'
223	PRINT *,' evap_tot= ',evap_tot
224	PRINT *,' airetot= ',airetot
225	IF (airetot .NE. airetot .OR. ABS(airetot) > largest*10.e15) THEN
226	varname = 'airephy'
227	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
228	END IF
229	IF (evap_tot .NE. evap_tot .OR. ABS(evap_tot) > largest) THEN
230	varname = 'evap'
231	CALL check_var(fname, varname, evap, klon, largest, .FALSE.)
232	END IF
233	PRINT *,' rain_fall_tot= ',rain_fall_tot
234	IF (rain_fall_tot .NE. rain_fall_tot .OR. ABS(rain_fall_tot) > largest) THEN
235	varname = 'rain_fall'
236	CALL check_var(fname, varname, rain_fall, klon, largest, .FALSE.)
237	END IF
238	PRINT *,' snow_fall_tot= ',snow_fall_tot
239	IF (snow_fall_tot .NE. snow_fall_tot .OR. ABS(snow_fall_tot) > largest) THEN
240	varname = 'snow_fall'
241	CALL check_var(fname, varname, snow_fall, klon, largest, .FALSE.)
242	END IF
243	STOP
244	END IF
245
246	IF (d_qt_tot .NE. d_qt_tot .OR. ABS(d_qt_tot) > largest) THEN
247	PRINT *,TRIM(errmsg)
248	PRINT *,' ' // TRIM(fname) // ': Wrong d_qt_tot= ',d_qt_tot,' !!!'
249	PRINT *,' d_qt_tot: Mass flux equivalent to atmospheric watter mass' // &
250	' change (kg/m2/s)'
251	PRINT *,' d_qt_tot = input value !'
252	STOP
253	END IF
254
255	return
256
257	6666 format('Phys. Flux Budget ',a15,1i6,2f8.2,2(1pE13.5))
258	6667 format('Phys. Boundary Flux ',a15,1i6,6f8.2,2(1pE13.5))
259	6668 format('Phys. Total Budget ',a15,1i6,f8.2,2(1pE13.5))
260
261	end SUBROUTINE diagphy
262
263	!L. Fita, LMD 2004. Check variable
264	SUBROUTINE check_var(funcn, varn, var, sizev, bigvalue, stoprun)
265	! Subroutine to check the consistency of a variable
266	! * NaN value: by definition is variable /= variable
267	! * bigvalue: threshold for the variable
268
269	IMPLICIT NONE
270
271	#include "dimensions.h"
272
273	INTEGER, INTENT(IN) :: sizev
274	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
275	REAL, DIMENSION(sizev), INTENT(IN) :: var
276	REAL, INTENT(IN) :: bigvalue
277	LOGICAL, INTENT(IN) :: stoprun
278
279	! Local
280	INTEGER :: i, wrongi, xpt, ypt
281	CHARACTER(LEN=50) :: errmsg
282	LOGICAL :: found
283	REAL, DIMENSION(sizev) :: wrongvalues
284	INTEGER, DIMENSION(sizev) :: wronggridpt
285
286	!!!!!!! Variables
287	! funcn: at which functino of part of the program variable is checked
288	! varn: name of the variable
289	! var: variable to check
290	! sizev: size of the variable
291	! bigvalue: biggest attenaible value for the variable
292	! stoprun: Should the run stop if it founds a problem?
293
294	errmsg = 'ERROR -- error -- ERROR -- error'
295
296	found = .FALSE.
297	wrongi = 0
298	DO i=1,sizev
299	IF (var(i) /= var(i) .OR. ABS(var(i)) > bigvalue ) THEN
300	IF (wrongi == 0) found = .TRUE.
301	wrongi = wrongi + 1
302	wrongvalues(wrongi) = var(i)
303	wronggridpt(wrongi) = i
304	END IF
305	END DO
306
307	IF (found) THEN
308	PRINT *,TRIM(errmsg)
309	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
310	"' is wrong in Nvalues= ",wrongi,' at i (x, y) value___'
311	DO i=1,wrongi
312	ypt = INT(wronggridpt(i)/wiim) + 1
313	xpt = wronggridpt(i) - (ypt-1)*wiim
314	PRINT *,wronggridpt(i), '(',xpt,', ',ypt,')', wrongvalues(i)
315	END DO
316	IF (stoprun) THEN
317	STOP
318	END IF
319	END IF
320
321	RETURN
322
323	END SUBROUTINE check_var
324
325	SUBROUTINE check_var3D(funcn, varn, var, sizev, zsize, bigvalue, stoprun)
326	! Subroutine to check the consistency of a 3D LMDSZ - variable (klon, klev) !
327	! * NaN value: by definition is variable /= variable
328	! * bigvalue: threshold for the variable
329
330	IMPLICIT NONE
331
332	#include "dimensions.h"
333
334	INTEGER, INTENT(IN) :: sizev, zsize
335	CHARACTER(LEN=50), INTENT(IN) :: funcn, varn
336	REAL, DIMENSION(sizev,zsize), INTENT(IN) :: var
337	REAL, INTENT(IN) :: bigvalue
338	LOGICAL, INTENT(IN) :: stoprun
339
340	! Local
341	INTEGER :: i, k, wrongi, xpt, ypt
342	CHARACTER(LEN=50) :: errmsg
343	LOGICAL :: found
344	REAL, DIMENSION(sizev*zsize) :: wrongvalues
345	INTEGER, DIMENSION(sizev*zsize,2) :: wronggridpt
346
347	!!!!!!! Variables
348	! funcn: at which functino of part of the program variable is checked
349	! varn: name of the variable
350	! var: variable to check
351	! sizev: size of the variable
352	! zsize: vertical size of the variable
353	! bigvalue: biggest attenaible value for the variable
354	! stoprun: Should the run stop if it founds a problem?
355
356	errmsg = 'ERROR -- error -- ERROR -- error'
357
358	found = .FALSE.
359	wrongi = 0
360	DO i=1,sizev
361	DO k=1,zsize
362	IF (var(i,k) /= var(i,k) .OR. ABS(var(i,k)) > bigvalue ) THEN
363	IF (wrongi == 0) found = .TRUE.
364	wrongi = wrongi + 1
365	wrongvalues(wrongi) = var(i,k)
366	wronggridpt(wrongi,1) = i
367	wronggridpt(wrongi,2) = k
368	END IF
369	END DO
370	END DO
371
372	IF (found) THEN
373	PRINT *,TRIM(errmsg)
374	PRINT *," at '" // TRIM(funcn) // "' variable '" //TRIM(varn)// &
375	"' is wrong in Nvalues= ",wrongi,' at i (x,y) k value___'
376	DO i=1,wrongi
377	ypt = INT(wronggridpt(i,1)/wiim) + 1
378	xpt = wronggridpt(i,1) - (ypt-1)*wiim
379	PRINT *,wronggridpt(i,1), '(',xpt,', ',ypt,')', wronggridpt(i,2), wrongvalues(i)
380	END DO
381	IF (stoprun) THEN
382	STOP
383	END IF
384	END IF
385
386	RETURN
387
388	END SUBROUTINE check_var3D
389
390
391	!C======================================================================
392	SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime &
393	& ,t,q,ql,qs,u,v,paprs,pplay &
394	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
395	!C======================================================================
396	!C
397	!C Purpose:
398	!C Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
399	!C et calcul le flux de chaleur et le flux d'eau necessaire a ces
400	!C changements. Ces valeurs sont moyennees sur la surface de tout
401	!C le globe et sont exprime en W/2 et kg/s/m2
402	!C Outil pour diagnostiquer la conservation de l'energie
403	!C et de la masse dans la physique. Suppose que les niveau de
404	!c pression entre couche ne varie pas entre 2 appels.
405	!C
406	!C Plusieurs de ces diagnostics peuvent etre fait en parallele: les
407	!c bilans sont sauvegardes dans des tableaux indices. On parlera
408	!C "d'indice de diagnostic"
409	!c
410	!C
411	!c======================================================================
412	!C Arguments:
413	!C airephy-------input-R- grid area
414	!C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
415	!C iprt----input-I- PRINT level ( <=1 : no PRINT)
416	!C idiag---input-I- indice dans lequel sera range les nouveaux
417	!C bilans d' entalpie et de masse
418	!C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
419	!C sont compare au bilan de d'enthalpie de masse de
420	!C l'indice numero idiag2
421	!C Cas parriculier : si idiag2=0, pas de comparaison, on
422	!c sort directement les bilans d'enthalpie et de masse
423	!C dtime----input-R- time step (s)
424	!c t--------input-R- temperature (K)
425	!c q--------input-R- vapeur d'eau (kg/kg)
426	!c ql-------input-R- liquid watter (kg/kg)
427	!c qs-------input-R- solid watter (kg/kg)
428	!c u--------input-R- vitesse u
429	!c v--------input-R- vitesse v
430	!c paprs----input-R- pression a intercouche (Pa)
431	!c pplay----input-R- pression au milieu de couche (Pa)
432	!c
433	!C the following total value are computed by UNIT of earth surface
434	!C
435	!C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
436	!c change (J/m2) during one time step (dtime) for the whole
437	!C atmosphere (air, watter vapour, liquid and solid)
438	!C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
439	!C total watter (kg/m2) change during one time step (dtime),
440	!C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
441	!C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
442	!C d_qs------output-R- same, for the solid watter only (kg/m2/s)
443	!C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
444	!C
445	!C other (COMMON...)
446	!C RCPD, RCPV, ....
447	!C
448	!C J.L. Dufresne, July 2002
449	!c======================================================================
450
451	USE dimphy
452	IMPLICIT NONE
453	!C
454	#include "dimensions.h"
455	!cccccc#include "dimphy.h"
456	#include "YOMCST.h"
457	#include "YOETHF.h"
458	!C
459	!c Input variables
460	real airephy(klon)
461	CHARACTER*15 tit
462	INTEGER iprt,idiag, idiag2
463	REAL dtime
464	REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
465	REAL u(klon,klev), v(klon,klev)
466	REAL paprs(klon,klev+1), pplay(klon,klev)
467	!c Output variables
468	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
469	!C
470	!C Local variables
471	!c
472	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
473	& , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
474	!c h_vcol_tot-- total enthalpy of vertical air column
475	!C (air with watter vapour, liquid and solid) (J/m2)
476	!c h_dair_tot-- total enthalpy of dry air (J/m2)
477	!c h_qw_tot---- total enthalpy of watter vapour (J/m2)
478	!c h_ql_tot---- total enthalpy of liquid watter (J/m2)
479	!c h_qs_tot---- total enthalpy of solid watter (J/m2)
480	!c qw_tot------ total mass of watter vapour (kg/m2)
481	!c ql_tot------ total mass of liquid watter (kg/m2)
482	!c qs_tot------ total mass of solid watter (kg/m2)
483	!c ec_tot------ total cinetic energy (kg/m2)
484	!C
485	REAL zairm(klon,klev) ! layer air mass (kg/m2)
486	REAL zqw_col(klon)
487	REAL zql_col(klon)
488	REAL zqs_col(klon)
489	REAL zec_col(klon)
490	REAL zh_dair_col(klon)
491	REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
492	!C
493	REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
494	!C
495	REAL airetot, zcpvap, zcwat, zcice
496	!C
497	INTEGER i, k
498	!C
499	INTEGER ndiag ! max number of diagnostic in parallel
500	PARAMETER (ndiag=10)
501	integer pas(ndiag)
502	save pas
503	data pas/ndiag*0/
504	!$OMP THREADPRIVATE(pas)
505	!C
506	REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag) &
507	& , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag) &
508	& , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
509	SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre &
510	& , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
511	!$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
512	!$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
513
514	! L. Fita, LMD July 2014
515	CHARACTER(LEN=50) :: errmsg, fname, varname
516	LOGICAL :: found
517	REAL :: largest
518
519	fname = 'diagetpq'
520	errmsg = 'ERROR -- error -- ERROR -- error'
521	largest = 10.e4
522	!c======================================================================
523	!C
524	DO k = 1, klev
525	DO i = 1, klon
526	!C layer air mass
527	zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
528	ENDDO
529	END DO
530	!C
531	!C Reset variables
532	DO i = 1, klon
533	zqw_col(i)=0.
534	zql_col(i)=0.
535	zqs_col(i)=0.
536	zec_col(i) = 0.
537	zh_dair_col(i) = 0.
538	zh_qw_col(i) = 0.
539	zh_ql_col(i) = 0.
540	zh_qs_col(i) = 0.
541	ENDDO
542	!C
543	zcpvap=RCPV
544	zcwat=RCW
545	zcice=RCS
546	!C
547	!C Compute vertical sum for each atmospheric column
548	!C ================================================
549	DO k = 1, klev
550	DO i = 1, klon
551	!C Watter mass
552	zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
553	zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
554	zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
555	!C Cinetic Energy
556	zec_col(i) = zec_col(i) &
557	& +0.5(u(i,k)2+v(i,k)2)zairm(i,k)
558	!C Air enthalpy
559	zh_dair_col(i) = zh_dair_col(i) &
560	& + RCPD(1.-q(i,k)-ql(i,k)-qs(i,k))zairm(i,k)*t(i,k)
561	zh_qw_col(i) = zh_qw_col(i) &
562	& + zcpvapq(i,k)zairm(i,k)*t(i,k)
563	zh_ql_col(i) = zh_ql_col(i) &
564	& + zcwatql(i,k)zairm(i,k)*t(i,k) &
565	& - RLVTTql(i,k)zairm(i,k)
566	zh_qs_col(i) = zh_qs_col(i) &
567	& + zciceqs(i,k)zairm(i,k)*t(i,k) &
568	& - RLSTTqs(i,k)zairm(i,k)
569
570	END DO
571	ENDDO
572
573	!C
574	!C Mean over the planete surface
575	!C =============================
576	qw_tot = 0.
577	ql_tot = 0.
578	qs_tot = 0.
579	ec_tot = 0.
580	h_vcol_tot = 0.
581	h_dair_tot = 0.
582	h_qw_tot = 0.
583	h_ql_tot = 0.
584	h_qs_tot = 0.
585	airetot=0.
586	!C
587	do i=1,klon
588	qw_tot = qw_tot + zqw_col(i)*airephy(i)
589	ql_tot = ql_tot + zql_col(i)*airephy(i)
590	qs_tot = qs_tot + zqs_col(i)*airephy(i)
591	ec_tot = ec_tot + zec_col(i)*airephy(i)
592	h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
593	h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
594	h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
595	h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
596	airetot=airetot+airephy(i)
597	END DO
598	!C
599	qw_tot = qw_tot/airetot
600	ql_tot = ql_tot/airetot
601	qs_tot = qs_tot/airetot
602	ec_tot = ec_tot/airetot
603	h_dair_tot = h_dair_tot/airetot
604	h_qw_tot = h_qw_tot/airetot
605	h_ql_tot = h_ql_tot/airetot
606	h_qs_tot = h_qs_tot/airetot
607	!C
608	h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
609	!C
610	!C Compute the change of the atmospheric state compare to the one
611	!C stored in "idiag2", and convert it in flux. THis computation
612	!C is performed IF idiag2 /= 0 and IF it is not the first CALL
613	!c for "idiag"
614	!C ===================================
615	!C
616	IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
617	d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
618	d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
619	d_h_qw = (h_qw_tot - h_qw_pre(idiag2) )/dtime
620	d_h_ql = (h_ql_tot - h_ql_pre(idiag2) )/dtime
621	d_h_qs = (h_qs_tot - h_qs_pre(idiag2) )/dtime
622	d_qw = (qw_tot - qw_pre(idiag2) )/dtime
623	d_ql = (ql_tot - ql_pre(idiag2) )/dtime
624	d_qs = (qs_tot - qs_pre(idiag2) )/dtime
625	d_ec = (ec_tot - ec_pre(idiag2) )/dtime
626	d_qt = d_qw + d_ql + d_qs
627
628	ELSE
629	d_h_vcol = 0.
630	d_h_dair = 0.
631	d_h_qw = 0.
632	d_h_ql = 0.
633	d_h_qs = 0.
634	d_qw = 0.
635	d_ql = 0.
636	d_qs = 0.
637	d_ec = 0.
638	d_qt = 0.
639	ENDIF
640	!C
641	IF (iprt.ge.2) THEN
642	WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
643	9000 format('Phys. Watter Mass Budget (kg/m2/s)',A15 &
644	& ,1i6,10(1pE14.6))
645	WRITE(6,9001) tit,pas(idiag), d_h_vcol
646	9001 format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
647	WRITE(6,9002) tit,pas(idiag), d_ec
648	9002 format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
649	END IF
650	!C
651	!C Store the new atmospheric state in "idiag"
652	!C
653	pas(idiag)=pas(idiag)+1
654	h_vcol_pre(idiag) = h_vcol_tot
655	h_dair_pre(idiag) = h_dair_tot
656	h_qw_pre(idiag) = h_qw_tot
657	h_ql_pre(idiag) = h_ql_tot
658	h_qs_pre(idiag) = h_qs_tot
659	qw_pre(idiag) = qw_tot
660	ql_pre(idiag) = ql_tot
661	qs_pre(idiag) = qs_tot
662	ec_pre (idiag) = ec_tot
663	!C
664	IF (d_h_vcol .NE. d_h_vcol .OR. ABS(d_h_vcol) > largest) THEN
665	PRINT *,TRIM(errmsg)
666	PRINT *,' ' // TRIM(fname) // ': Wrong d_h_vcol= ',d_h_vcol,' !!!'
667	PRINT *,' d_h_vcol: Heat flux (W/m2) define as the Enthalpy change' // &
668	' (J/m2) during one time step (dtime) for the whole atmosphere (air,' // &
669	' watter vapour, liquid and solid)'
670	PRINT *,' d_h_vcol = (h_vcol_tot - h_vcol_pre(idiag2) )'
671	PRINT *,' h_vcol_tot= ',h_vcol_tot
672	IF (h_vcol_tot .NE. h_vcol_tot .OR. ABS(h_vcol_tot) > largest) THEN
673	PRINT *,' h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot'
674	PRINT *,' airetot= ',airetot
675	IF (airetot .NE. airetot .OR. ABS(airetot) > largest*10.e15) THEN
676	varname = 'airephy'
677	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
678	END IF
679	PRINT *,' h_dair_tot= ',h_dair_tot
680	IF (h_dair_tot .NE. h_dair_tot .OR. ABS(h_dair_tot) > largest*10.e6) THEN
681	varname = 'zh_dair_col'
682	CALL check_var(fname, varname, zh_dair_col, klon, largest*10.e6, .FALSE.)
683	PRINT *,' zh_dair_col = func(q,ql,qs,paprs,t)'
684	varname = 'q'
685	CALL check_var3D(fname, varname, q, klon, klev, largest, .FALSE.)
686	varname = 'ql'
687	CALL check_var3D(fname, varname, ql, klon, klev, largest, .FALSE.)
688	varname = 'qs'
689	CALL check_var3D(fname, varname, qs, klon, klev, largest, .FALSE.)
690	varname = 'paprs'
691	CALL check_var3D(fname, varname, paprs, klon, klev, largest*10.e6, .FALSE.)
692	varname = 't'
693	CALL check_var3D(fname, varname, t, klon, klev, largest, .FALSE.)
694	END IF
695	PRINT *,' h_qw_tot= ',h_qw_tot
696	IF (h_qw_tot .NE. h_qw_tot .OR. ABS(h_qw_tot) > largest*10.e4) THEN
697	varname = 'zh_qw_col'
698	CALL check_var(fname, varname, zh_qw_col, klon, largest*10.e4, .FALSE.)
699	END IF
700	PRINT *,' h_ql_tot= ',h_ql_tot
701	IF (h_ql_tot .NE. h_ql_tot .OR. ABS(h_ql_tot) > largest*10.e2) THEN
702	varname = 'zh_ql_col'
703	CALL check_var(fname, varname, zh_ql_col, klon, largest*10.e2, .FALSE.)
704	END IF
705	PRINT *,' h_qs_tot= ',h_qs_tot
706	IF (h_qs_tot .NE. h_qs_tot .OR. ABS(h_qs_tot) > largest) THEN
707	varname = 'zh_qs_col'
708	CALL check_var(fname, varname, zh_qs_col, klon, largest, .FALSE.)
709	END IF
710	END IF
711	STOP
712	END IF
713
714	IF (qw_tot .NE. qw_tot .OR. ABS(qw_tot) > largest) THEN
715	PRINT *,TRIM(errmsg)
716	PRINT *,' ' // TRIM(fname) // ': Wrong qw_tot= ',qw_tot,' !!!'
717	PRINT *,' qw_tot: total mass of watter vapour (kg/m2)'
718	PRINT *,' qw_tot = f(q,ql,qs,airephy,t,paprs)'
719	varname = 'airephy'
720	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
721	varname = 'q'
722	CALL check_var3D(fname, varname, q, klon, klev, largest, .FALSE.)
723	varname = 'ql'
724	CALL check_var3D(fname, varname, ql, klon, klev, largest, .FALSE.)
725	varname = 'qs'
726	CALL check_var3D(fname, varname, qs, klon, klev, largest, .FALSE.)
727	varname = 't'
728	CALL check_var3D(fname, varname, t, klon, klev, largest, .FALSE.)
729	varname = 'paprs'
730	CALL check_var3D(fname, varname, paprs, klon, klev, largest, .FALSE.)
731	STOP
732	END IF
733
734	IF (ql_tot .NE. ql_tot .OR. ABS(ql_tot) > largest) THEN
735	PRINT *,TRIM(errmsg)
736	PRINT *,' ' // TRIM(fname) // ': Wrong ql_tot= ',ql_tot,' !!!'
737	PRINT *,' ql_tot: total mass of liquid watter (kg/m2)'
738	PRINT *,' ql_tot = f(ql,airephy,t,paprs)'
739	varname = 'airephy'
740	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
741	varname = 'ql'
742	CALL check_var3D(fname, varname, ql, klon, klev, largest, .FALSE.)
743	varname = 't'
744	CALL check_var3D(fname, varname, t, klon, klev, largest, .FALSE.)
745	varname = 'paprs'
746	CALL check_var3D(fname, varname, paprs, klon, klev, largest, .FALSE.)
747	STOP
748	END IF
749
750	IF (qs_tot .NE. qs_tot .OR. ABS(qs_tot) > largest) THEN
751	PRINT *,TRIM(errmsg)
752	PRINT *,' ' // TRIM(fname) // ': Wrong qs_tot= ',qs_tot,' !!!'
753	PRINT *,' qs_tot: total mass of solid watter (kg/m2)'
754	PRINT *,' qs_tot = f(qs,airephy,t,paprs)'
755	varname = 'airephy'
756	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
757	varname = 'qs'
758	CALL check_var3D(fname, varname, qs, klon, klev, largest, .FALSE.)
759	varname = 't'
760	CALL check_var3D(fname, varname, t, klon, klev, largest, .FALSE.)
761	varname = 'paprs'
762	CALL check_var3D(fname, varname, paprs, klon, klev, largest, .FALSE.)
763	STOP
764	END IF
765
766	IF (ec_tot .NE. ec_tot .OR. ABS(ec_tot) > largest*10.e4) THEN
767	PRINT *,TRIM(errmsg)
768	PRINT *,' ' // TRIM(fname) // ': Wrong ec_tot= ',ec_tot,' !!!'
769	PRINT *,' ec_tot: total cinetic energy (kg/m2)'
770	PRINT *,' ec_tot = f(u,v,airephy,t,paprs)'
771	varname = 'airephy'
772	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
773	varname = 'u'
774	CALL check_var3D(fname, varname, u, klon, klev, largest, .FALSE.)
775	varname = 'v'
776	CALL check_var3D(fname, varname, v, klon, klev, largest, .FALSE.)
777	varname = 't'
778	CALL check_var3D(fname, varname, t, klon, klev, largest, .FALSE.)
779	varname = 'paprs'
780	CALL check_var3D(fname, varname, paprs, klon, klev, largest, .FALSE.)
781	STOP
782	END IF
783
784	RETURN
785
786	END SUBROUTINE diagetpq
787
788	END MODULE diagphy_mod
789

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

Download in other formats:

Original Format