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

Rev	Line
[1]	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
[142]	81
	82	save pas
	83	data pas/0/
	84	!$OMP THREADPRIVATE(pas)
	85	!C
	86
[31]	87	! L. Fita, LMD July 2014
[142]	88	CHARACTER(LEN=50) :: errmsg, fname, varname
	89	LOGICAL :: found
	90	REAL :: largest
[31]	91
	92	fname = 'diagphy'
	93	errmsg = 'ERROR -- error -- ERROR -- error'
[142]	94	largest = 10.e4
[31]	95
[1]	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
[31]	159	! L. Fita, LMD July 2014. Checking for consistency
[142]	160	IF (fs_bound .NE. fs_bound .OR. ABS(fs_bound) > largest) THEN
[31]	161	PRINT *,TRIM(errmsg)
[142]	162	PRINT *,' ' // TRIM(fname) // ': Wrong fs_bound= ',fs_bound,' !!!'
[31]	163	PRINT *,' fs_bound: Heat flux at atm. boundaries'
	164	PRINT *,' fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront'
[142]	165	PRINT *,' airetot= ',airetot
[240]	166	IF (airetot .NE. airetot .OR. ABS(airetot) > largest*10.e16) THEN
[142]	167	varname = 'airephy'
[240]	168	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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
[31]	180	PRINT *,' ssols= ',ssols
[142]	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
[31]	185	PRINT *,' ssoll= ',ssoll
[142]	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
[31]	190	PRINT *,' ssens= ',ssens
[142]	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
[31]	195	PRINT *,' sfront= ',sfront
[142]	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
[31]	206	STOP
	207	END IF
	208
[142]	209	IF (d_etp_tot .NE. d_etp_tot .OR. ABS(d_etp_tot) > largest) THEN
[31]	210	PRINT *,TRIM(errmsg)
[142]	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,' !!!'
[31]	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
[142]	224	PRINT *,' airetot= ',airetot
[240]	225	IF (airetot .NE. airetot .OR. ABS(airetot) > largest*10.e15) THEN
[142]	226	varname = 'airephy'
[240]	227	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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
[31]	233	PRINT *,' rain_fall_tot= ',rain_fall_tot
[142]	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
[31]	238	PRINT *,' snow_fall_tot= ',snow_fall_tot
[142]	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
[31]	243	STOP
	244	END IF
	245
[142]	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
[1]	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
[142]	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
[1]	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)
[142]	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
[1]	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
[142]	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
[240]	675	IF (airetot .NE. airetot .OR. ABS(airetot) > largest*10.e15) THEN
[142]	676	varname = 'airephy'
[240]	677	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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'
[240]	720	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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'
[240]	740	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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'
[240]	756	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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'
[240]	772	CALL check_var(fname, varname, airephy, klon, largest*10.e15, .FALSE.)
[142]	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
[1]	784	RETURN
[142]	785
[1]	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