Context Navigation

add_phys_tend_mod.F90 @ 4797

Last change on this file since 4797 was 4523, checked in by evignon, 20 months ago
merge de la branche blowing snow vers la trunk premiere tentative Etienne
Property svn:keywords set to `Id`
File size: 36.0 KB

Rev	Line
[3927]	1	!
	2	! $Id: add_phys_tend_mod.F90 4523 2023-05-03 16:21:08Z nfevrier $
	3	!
	4	!
	5	MODULE add_phys_tend_mod
	6
	7	IMPLICIT NONE
	8	! flag to compute diagnostics to check energy conservation. If fl_ebil==0, no check
	9	INTEGER, SAVE :: fl_ebil
	10	!$OMP THREADPRIVATE(fl_ebil)
	11	! flag to include modifcations to ensure energy conservation. If fl_cor_ebil==0, no corrections
	12	! Note that with time, all these modifications should be included by default
	13	INTEGER, SAVE :: fl_cor_ebil
	14	!$OMP THREADPRIVATE(fl_cor_ebil)
	15
	16	CONTAINS
	17
[4523]	18	SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap &
[3927]	19	#ifdef ISO
	20	,zdxt,zdxtl,zdxti &
	21	#endif
	22	)
	23	! ======================================================================
	24	! Ajoute les tendances de couche limite, soit determinees par la
	25	! parametrisation
	26	! physique, soit forcees, aux variables d etat de la dynamique t_seri,
	27	! q_seri ...
	28	! ======================================================================
	29
	30
	31	! ======================================================================
	32	! Declarations
	33	! ======================================================================
	34
	35	USE dimphy, ONLY: klon, klev
	36	! USE dimphy
	37	USE phys_local_var_mod
	38	USE phys_state_var_mod
	39	USE mod_grid_phy_lmdz, ONLY: nbp_lev
	40	#ifdef ISO
[4143]	41	USE infotrac_phy, ONLY: ntraciso=>ntiso
[3927]	42	USE isotopes_mod, ONLY: iso_eau
	43	#endif
	44	IMPLICIT NONE
	45	REAL,SAVE,ALLOCATABLE :: hthturb_gcssold(:)
	46	REAL,SAVE,ALLOCATABLE :: hqturb_gcssold(:)
	47	!$OMP THREADPRIVATE(hthturb_gcssold,hqturb_gcssold)
	48	REAL,SAVE :: dtime_frcg
	49	LOGICAL,SAVE :: turb_fcg_gcssold
	50	LOGICAL,SAVE :: firstcall=.true.
	51	!$OMP THREADPRIVATE(firstcall,turb_fcg_gcssold,dtime_frcg)
	52	INTEGER abortphy
	53	! COMMON /turb_forcing/dtime_frcg, hthturb_gcssold, hqturb_gcssold, &
	54	! turb_fcg_gcssold
	55	#ifdef ISO
	56	REAL,SAVE,ALLOCATABLE :: hxtturb_gcssold(:,:)
	57	!$OMP THREADPRIVATE(hxtturb_gcssold)
	58	! variable à définir dans le 1D
	59	#endif
	60
	61	! Arguments :
	62	! ------------
	63	REAL zdu(klon, klev), zdv(klon, klev)
[4523]	64	REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev), zdqbs(klon,klev)
[3927]	65	CHARACTER () text
	66	REAL paprs(klon,klev+1)
	67	INTEGER flag_inhib_tend ! if flag_inhib_tend != 0, tendencies are not added
	68	INTEGER itap ! time step number
	69	#ifdef ISO
	70	REAL zdxt(ntraciso,klon, klev), zdxtl(ntraciso,klon, klev), zdxti(ntraciso,klon, klev)
	71	#endif
	72
	73	! Local :
	74	! --------
	75	REAL zzdt(klon, klev), zzdq(klon, klev)
	76	INTEGER i, k
	77	#ifdef ISO
	78	REAL zzdxt(ntraciso,klon, klev)
	79	INTEGER ixt
	80	#endif
	81
	82	IF (firstcall) THEN
	83	ALLOCATE(hthturb_gcssold(nbp_lev))
	84	ALLOCATE(hqturb_gcssold(nbp_lev))
	85	#ifdef ISO
	86	ALLOCATE(hxtturb_gcssold(ntraciso,nbp_lev))
	87	#endif
	88	firstcall=.false.
	89	ENDIF
	90
	91	IF (turb_fcg_gcssold) THEN
	92	DO k = 1, klev
	93	DO i = 1, klon
	94	zzdt(i, k) = hthturb_gcssold(k)*dtime_frcg
	95	zzdq(i, k) = hqturb_gcssold(k)*dtime_frcg
	96	#ifdef ISO
	97	do ixt=1,ntraciso
	98	zzdxt(ixt,i, k) = hxtturb_gcssold(ixt,k)*dtime_frcg
	99	enddo
	100	#endif
	101	END DO
	102	END DO
	103	PRINT *, ' add_pbl_tend, dtime_frcg ', dtime_frcg
	104	PRINT *, ' add_pbl_tend, zzdt ', zzdt
	105	PRINT *, ' add_pbl_tend, zzdq ', zzdq
[4523]	106	CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
[3927]	107	#ifdef ISO
	108	,zzdxt,zdxtl,zdxti &
	109	#endif
	110	)
	111	ELSE
[4523]	112	CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
[3927]	113	#ifdef ISO
	114	,zdxt,zdxtl,zdxti &
	115	#endif
	116	)
	117	END IF
	118
	119
	120	RETURN
	121	END SUBROUTINE add_pbl_tend
	122	!
[4004]	123	! $Id: add_phys_tend_mod.F90 4523 2023-05-03 16:21:08Z nfevrier $
[3927]	124	!
[4523]	125	SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,zdqbs,paprs,text, &
[3927]	126	abortphy,flag_inhib_tend, itap, diag_mode &
	127	#ifdef ISO
	128	,zdxt,zdxtl,zdxti &
	129	#endif
	130	)
	131	!======================================================================
	132	! Ajoute les tendances des variables physiques aux variables
	133	! d'etat de la dynamique t_seri, q_seri ...
	134	! On en profite pour faire des tests sur les tendances en question.
	135	!======================================================================
	136
	137
	138	!======================================================================
	139	! Declarations
	140	!======================================================================
	141
	142	USE dimphy, ONLY: klon, klev
	143	USE phys_state_var_mod, ONLY : phys_tstep
[4523]	144	USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, qbs_seri, q_seri, t_seri
[3927]	145	#ifdef ISO
	146	USE phys_local_var_mod, ONLY: xtl_seri, xts_seri, xt_seri
	147	#endif
	148	USE phys_state_var_mod, ONLY: ftsol
	149	USE geometry_mod, ONLY: longitude_deg, latitude_deg
	150	USE print_control_mod, ONLY: prt_level
	151	USE cmp_seri_mod
[4523]	152	USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
	153	& , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
[3927]	154
	155	#ifdef ISO
[4143]	156	USE infotrac_phy, ONLY: ntraciso=>ntiso
[3927]	157	#ifdef ISOVERIF
	158	USE isotopes_mod, ONLY: iso_eau
	159	USE isotopes_verif_mod
	160	#endif
	161	#endif
	162	IMPLICIT none
	163	include "YOMCST.h"
	164	include "clesphys.h"
	165
	166	! Arguments :
	167	!------------
	168	REAL, DIMENSION(klon,klev), INTENT(IN) :: zdu, zdv
[4523]	169	REAL, DIMENSION(klon,klev), INTENT(IN) :: zdt, zdql, zdqi, zdqbs
[3927]	170	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
	171	CHARACTER(), INTENT(IN) :: text
	172	INTEGER, INTENT(IN) :: abortphy
	173	INTEGER, INTENT(IN) :: flag_inhib_tend ! if not 0, tendencies are not added
	174	INTEGER, INTENT(IN) :: itap ! time step number
	175	INTEGER, INTENT(IN) :: diag_mode ! 0 -> normal effective mode
	176	! 1 -> only conservation stats are computed
	177	!
	178	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zdq
	179	#ifdef ISO
	180	REAL zdxt(ntraciso,klon,klev),zdxtl(ntraciso,klon,klev) &
	181	,zdxti(ntraciso,klon,klev)
	182	integer ixt
	183	#endif
	184
	185	! Local :
	186	!--------
	187	REAL zt,zq
	188	#ifdef ISO
	189	REAL zxt(ntraciso)
	190	real qprec(klon,klev)
	191	REAL zxt_int(ntraciso), zxtp_int(ntraciso), zxtp(ntraciso,klev),zxt_new
	192	#endif
	193	REAL zq_int, zqp_int, zq_new
	194
	195	REAL zqp(klev)
	196
	197	! Save variables, used in diagnostic mode (diag_mode=1).
	198	REAL, DIMENSION(klon,klev) :: sav_u_seri, sav_v_seri
[4523]	199	REAL, DIMENSION(klon,klev) :: sav_ql_seri, sav_qs_seri, sav_qbs_seri, sav_q_seri
[3927]	200	REAL, DIMENSION(klon,klev) :: sav_t_seri
	201	REAL, DIMENSION(klon,klev) :: sav_zdq
	202	#ifdef ISO
	203	REAL, DIMENSION(ntraciso,klon,klev) :: sav_zdxt
	204	REAL, DIMENSION(ntraciso,klon,klev) :: sav_xtl_seri, sav_xts_seri, sav_xt_seri
	205	#endif
	206	!
	207	INTEGER i, k,j, n
	208	INTEGER jadrs(klon*klev), jbad
	209	INTEGER jqadrs(klon*klev), jqbad
	210	INTEGER kadrs(klon*klev)
	211	INTEGER kqadrs(klon*klev)
	212
	213	LOGICAL done(klon)
	214
	215	integer debug_level
	216	logical, save :: first=.true.
	217	!$OMP THREADPRIVATE(first)
	218	!
	219	!======================================================================
	220	! Variables for energy conservation tests
	221	!======================================================================
	222	!
	223
	224	! zh_col------- total enthalpy of vertical air column
	225	! (air with watter vapour, liquid and solid) (J/m2)
	226	! zh_dair_col--- total enthalpy of dry air (J/m2)
	227	! zh_qw_col---- total enthalpy of watter vapour (J/m2)
	228	! zh_ql_col---- total enthalpy of liquid watter (J/m2)
	229	! zh_qs_col---- total enthalpy of solid watter (J/m2)
[4523]	230	! zh_qbs_col---- total enthalpy of blowing snow (J/m2)
[3927]	231	! zqw_col------ total mass of watter vapour (kg/m2)
	232	! zql_col------ total mass of liquid watter (kg/m2)
[4523]	233	! zqs_col------ total mass of cloud ice (kg/m2)
	234	! zqbs_col------ total mass of blowing snow (kg/m2)
[3927]	235	! zek_col------ total kinetic energy (kg/m2)
	236	!
	237	REAL zairm(klon, klev) ! layer air mass (kg/m2)
	238	REAL zqw_col(klon,2)
	239	REAL zql_col(klon,2)
	240	REAL zqs_col(klon,2)
[4523]	241	REAL zqbs_col(klon,2)
[3927]	242	REAL zek_col(klon,2)
	243	REAL zh_dair_col(klon,2)
[4523]	244	REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2), zh_qbs_col(klon,2)
[3927]	245	REAL zh_col(klon,2)
	246
	247	REAL zcpvap, zcwat, zcice
	248
	249	!======================================================================
	250	! Initialisations
	251
	252	IF (prt_level >= 5) then
	253	write (,) "In add_phys_tend, after ",text
	254	call flush
	255	end if
	256
	257	! if flag_inhib_tend != 0, tendencies are not added
	258	IF (flag_inhib_tend /= 0) then
	259	! If requiered, diagnostics are shown
	260	IF (flag_inhib_tend > 0) then
	261	! print some diagnostics if xxx_seri have changed
	262	call cmp_seri(flag_inhib_tend,text)
	263	END IF
	264	RETURN ! on n ajoute pas les tendance
	265	END IF
	266
	267	IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
	268	! a deja plante.
	269
	270	debug_level=10
	271	if (first) then
	272	print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
	273	first=.false.
	274	endif
	275	!
	276	! print *,'add_phys_tend: paprs ',paprs
	277	! When in diagnostic mode, save initial values of out variables
	278	IF (diag_mode == 1) THEN
	279	sav_u_seri(:,:) = u_seri(:,:)
	280	sav_v_seri(:,:) = v_seri(:,:)
	281	sav_ql_seri(:,:) = ql_seri(:,:)
	282	sav_qs_seri(:,:) = qs_seri(:,:)
[4523]	283	sav_qbs_seri(:,:) = qbs_seri(:,:)
[3927]	284	sav_q_seri(:,:) = q_seri(:,:)
	285	sav_t_seri(:,:) = t_seri(:,:)
	286	sav_zdq(:,:) = zdq(:,:)
	287	#ifdef ISO
	288	sav_zdxt(:,:,:) = zdxt(:,:,:)
	289	#endif
	290	ENDIF ! (diag_mode == 1)
	291	!======================================================================
	292	! Diagnostics for energy conservation tests
	293	!======================================================================
	294	DO k = 1, klev
	295	! layer air mass
	296	zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
	297	END DO
	298
	299	if (fl_ebil .GT. 0) then
	300	! ------------------------------------------------
	301	! Compute vertical sum for each atmospheric column
	302	! ------------------------------------------------
	303	n=1 ! begining of time step
	304
	305	zcpvap = rcpv
	306	zcwat = rcw
	307	zcice = rcs
	308
	309	CALL integr_v(klon, klev, zcpvap, &
[4523]	310	t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, zairm, &
	311	zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
	312	zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
[3927]	313
	314	end if ! end if (fl_ebil .GT. 0)
	315
	316	!======================================================================
	317	! Ajout des tendances sur le vent et l'eau liquide
	318	!======================================================================
	319
[4491]	320	#ifdef ISO
	321	#ifdef ISOVERIF
	322	if (iso_eau.gt.0) then
	323	call iso_verif_egalite_vect2D( &
	324	xt_seri,q_seri, &
	325	'add_phys_tend 321a: '//text,ntraciso,klon,klev)
	326	call iso_verif_egalite_vect2D( &
	327	xtl_seri,ql_seri, &
	328	'add_phys_tend 321b: '//text,ntraciso,klon,klev)
	329	endif !if (iso_eau.gt.0) then
	330	#ifdef ISOTRAC
	331	call iso_verif_traceur_vect(xt_seri,klon,klev, &
	332	& 'add_phys_tend 328a: '//text)
	333	call iso_verif_traceur_vect(xtl_seri,klon,klev, &
	334	& 'add_phys_tend 328b: '//text)
	335	#endif
	336	#endif
	337	#endif
	338
[3927]	339	u_seri(:,:)=u_seri(:,:)+zdu(:,:)
	340	v_seri(:,:)=v_seri(:,:)+zdv(:,:)
	341	ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
	342	qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
[4523]	343	qbs_seri(:,:)=qbs_seri(:,:)+zdqbs(:,:)
[3927]	344	#ifdef ISO
	345	xtl_seri(:,:,:)=xtl_seri(:,:,:)+zdxtl(:,:,:)
	346	xts_seri(:,:,:)=xts_seri(:,:,:)+zdxti(:,:,:)
	347	#endif
	348
	349	!======================================================================
	350	! On ajoute les tendances de la temperature et de la vapeur d'eau
	351	! en verifiant que ca ne part pas dans les choux
	352	!======================================================================
	353
	354	jbad=0
	355	jqbad=0
	356	DO k = 1, klev
	357	DO i = 1, klon
	358	zt=t_seri(i,k)+zdt(i,k)
	359	zq=q_seri(i,k)+zdq(i,k)
	360	#ifdef ISO
	361	do ixt=1,ntraciso
	362	zxt(ixt)=xt_seri(ixt,i,k)+zdxt(ixt,i,k)
	363	enddo !do ixt=1,ntraciso
	364	qprec(i,k)=q_seri(i,k)
	365	#endif
	366	IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
	367	jbad = jbad + 1
	368	jadrs(jbad) = i
	369	kadrs(jbad) = k
	370	ENDIF
	371	IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
	372	jqbad = jqbad + 1
	373	jqadrs(jqbad) = i
	374	kqadrs(jqbad) = k
	375	ENDIF
	376	t_seri(i,k)=zt
	377	q_seri(i,k)=zq
	378	#ifdef ISO
	379	do ixt=1,ntraciso
	380	xt_seri(ixt,i,k)=zxt(ixt)
	381	enddo !do ixt=1,ntraciso
	382	#endif
	383	ENDDO
	384	ENDDO
	385
	386	#ifdef ISO
	387	#ifdef ISOVERIF
	388	if (iso_eau.gt.0) then
	389	call iso_verif_egalite_vect2D( &
	390	xt_seri,q_seri, &
	391	'add_phys_tend 122: '//text,ntraciso,klon,klev)
	392	call iso_verif_egalite_vect2D( &
	393	xtl_seri,ql_seri, &
	394	'add_phys_tend 138: '//text,ntraciso,klon,klev)
	395	endif !if (iso_eau.gt.0) then
[4491]	396	#ifdef ISOTRAC
	397	call iso_verif_traceur_vect(xt_seri,klon,klev, &
	398	& 'add_phys_tend 374a: '//text)
	399	call iso_verif_traceur_vect(xtl_seri,klon,klev, &
	400	& 'add_phys_tend 374b: '//text)
[3927]	401	#endif
	402	#endif
[4491]	403	#endif
[3927]	404
	405
	406	!=====================================================================================
	407	! Impression et stop en cas de probleme important
	408	!=====================================================================================
	409
	410	IF (jbad .GT. 0) THEN
	411	DO j = 1, jbad
	412	i=jadrs(j)
	413	if(prt_level.ge.debug_level) THEN
	414	print*,'PLANTAGE POUR LE POINT i lon lat =',&
	415	i,longitude_deg(i),latitude_deg(i),text
	416	print*,'l T dT Q dQ '
	417	DO k = 1, klev
	418	write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
	419	ENDDO
	420	call print_debug_phys(i,debug_level,text)
	421	endif
	422	ENDDO
	423	ENDIF
	424	!
	425	!=====================================================================================
	426	! Impression, warning et correction en cas de probleme moins important
	427	!=====================================================================================
	428	IF (jqbad .GT. 0) THEN
	429	done(:) = .false. !jyg
	430	DO j = 1, jqbad
	431	i=jqadrs(j)
	432	if(prt_level.ge.debug_level) THEN
	433	print*,'WARNING : EAU POUR LE POINT i lon lat =',&
	434	i,longitude_deg(i),latitude_deg(i),text
	435	print*,'l T dT Q dQ '
	436	DO k = 1, klev
	437	write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
	438	ENDDO
	439	endif
	440	IF (ok_conserv_q) THEN
	441	!jyg<20140228 Corrections pour conservation de l'eau
	442	IF (.NOT.done(i)) THEN !jyg
	443	DO k = 1, klev
	444	zqp(k) = max(q_seri(i,k),1.e-15)
	445	#ifdef ISO
	446	do ixt=1,ntraciso
	447	zxtp(ixt,k)=xt_seri(ixt,i,k)*max(1.0,1.e-15/qprec(i,k))
	448	enddo
	449	#endif
	450	ENDDO
	451	zq_int = 0.
	452	zqp_int = 0.
	453	DO k = 1, klev
	454	zq_int = zq_int + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
	455	zqp_int = zqp_int + zqp(k) *(paprs(i,k)-paprs(i,k+1))/Rg
	456	ENDDO
	457	#ifdef ISO
	458	do ixt=1,ntraciso
	459	zxt_int(ixt) = 0.
	460	zxtp_int(ixt) = 0.
	461	DO k = 1, klev
	462	zxt_int(ixt) = zxt_int(ixt) + xt_seri(ixt,i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
	463	zxtp_int(ixt) = zxtp_int(ixt) + zxtp(ixt,k) *(paprs(i,k)-paprs(i,k+1))/Rg
	464	ENDDO
	465	enddo
	466	#endif
	467	if(prt_level.ge.debug_level) THEN
	468	print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
	469	i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
	470	endif
	471	DO k = 1, klev
	472	zq_new = zqp(k)*zq_int/zqp_int
	473	zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
	474	q_seri(i,k) = zq_new
	475	#ifdef ISO
	476	do ixt=1,ntraciso
	477	zxt_new = zxtp(ixt,k)*zxt_int(ixt)/zxtp_int(ixt)
	478	zdxt(ixt,i,k) = zdxt(ixt,i,k) + zxt_new - xt_seri(ixt,i,k)
	479	xt_seri(ixt,i,k) = zxt_new
	480	enddo !ixt=1,ntraciso
	481	#endif
	482	ENDDO
	483	done(i) = .true.
	484	ENDIF !(.NOT.done(i))
	485	ELSE
	486	!jyg>
	487	DO k = 1, klev
	488	zq=q_seri(i,k)+zdq(i,k)
	489	if (zq.lt.1.e-15) then
	490	if (q_seri(i,k).lt.1.e-15) then
	491	if(prt_level.ge.debug_level) THEN
	492	print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
	493	endif
	494	q_seri(i,k)=1.e-15
	495	zdq(i,k)=(1.e-15-q_seri(i,k))
	496	#ifdef ISO
	497	do ixt=1,ntraciso
	498	xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
	499	! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
	500	! zdxt(ixt,i,k)=0.0
	501	zdxt(ixt,i,k)=xt_seri(ixt,i,k)*(1.e-15/qprec(i,k)-1)
	502	enddo !do ixt=1,ntraciso
	503	#endif
	504	endif
	505	endif
	506	! zq=q_seri(i,k)+zdq(i,k)
	507	! if (zq.lt.1.e-15) then
	508	! zdq(i,k)=(1.e-15-q_seri(i,k))
	509	! endif
	510	ENDDO
	511	!jyg<20140228
	512	ENDIF ! (ok_conserv_q)
	513	!jyg>
	514	ENDDO ! j = 1, jqbad
	515	ENDIF
	516
	517	#ifdef ISO
	518	#ifdef ISOVERIF
	519	if (iso_eau.gt.0) then
	520	call iso_verif_egalite_vect2D( &
	521	xt_seri,q_seri, &
	522	'add_phys_tend 173'//text,ntraciso,klon,klev)
	523	call iso_verif_egalite_vect2D( &
	524	zdxt,zdq, &
	525	'add_phys_tend 175'//text,ntraciso,klon,klev)
	526	endif !if (iso_eau.gt.0) then
[4491]	527	#ifdef ISOTRAC
	528	call iso_verif_traceur_vect(xt_seri,klon,klev, &
	529	& 'add_phys_tend 499'//text)
[3927]	530	#endif
	531	#endif
[4491]	532	#endif
[3927]	533	!
	534
	535	!IM ajout memes tests pour reverifier les jbad, jqbad beg
	536	jbad=0
	537	jqbad=0
	538	DO k = 1, klev
	539	DO i = 1, klon
	540	IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
	541	jbad = jbad + 1
	542	jadrs(jbad) = i
	543	! if(prt_level.ge.debug_level) THEN
	544	! print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
	545	! endif
	546	ENDIF
	547	IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
	548	jqbad = jqbad + 1
	549	jqadrs(jqbad) = i
	550	kqadrs(jqbad) = k
	551	! if(prt_level.ge.debug_level) THEN
	552	! print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
	553	! endif
	554	ENDIF
	555	ENDDO
	556	ENDDO
	557	IF (jbad .GT. 0) THEN
	558	DO j = 1, jbad
	559	i=jadrs(j)
	560	k=kadrs(j)
	561	if(prt_level.ge.debug_level) THEN
	562	print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
	563	i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
	564	& zdt(i,k),t_seri(i,k)-zdt(i,k)
	565	!!! if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
	566	print*,'l T dT Q dQ '
	567	DO k = 1, klev
	568	write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
	569	ENDDO
	570	call print_debug_phys(i,debug_level,text)
	571	endif
	572	ENDDO
	573	ENDIF
	574	!
	575	IF (jqbad .GT. 0) THEN
	576	DO j = 1, jqbad
	577	i=jqadrs(j)
	578	k=kqadrs(j)
	579	if(prt_level.ge.debug_level) THEN
	580	print*,'WARNING : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
	581	i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
	582	& zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
	583	!!! if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
	584	print*,'l T dT Q dQ '
	585	DO k = 1, klev
	586	write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
	587	ENDDO
	588	call print_debug_phys(i,debug_level,text)
	589	endif
	590	ENDDO
	591	ENDIF
	592
	593	!======================================================================
	594	! Contrôle des min/max pour arrêt du modèle
	595	! Si le modele est en mode abortphy, on retire les tendances qu'on
	596	! vient d'ajouter. Pas exactement parce qu'on ne tient pas compte des
	597	! seuils.
	598	!======================================================================
	599
	600	CALL hgardfou(t_seri,ftsol,text,abortphy)
	601	IF (abortphy==1) THEN
	602	Print*,'ERROR ABORT hgardfou dans ',text
	603	! JLD pourquoi on ne modifie pas de meme t_seri et q_seri ?
	604	u_seri(:,:)=u_seri(:,:)-zdu(:,:)
	605	v_seri(:,:)=v_seri(:,:)-zdv(:,:)
	606	ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
	607	qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
[4523]	608	qbs_seri(:,:)=qbs_seri(:,:)-zdqbs(:,:)
[3927]	609	ENDIF
	610
	611	!======================================================================
	612	! Diagnostics for energy conservation tests
	613	!======================================================================
	614
	615	if (fl_ebil .GT. 0) then
	616
	617	! ------------------------------------------------
	618	! Compute vertical sum for each atmospheric column
	619	! ------------------------------------------------
	620	n=2 ! end of time step
	621
	622	CALL integr_v(klon, klev, zcpvap, &
[4523]	623	t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, zairm, &
	624	zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
	625	zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
[3927]	626
	627	! ------------------------------------------------
	628	! Compute the changes by unit of time
	629	! ------------------------------------------------
	630
	631	d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/phys_tstep
	632	d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
	633	d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
[4523]	634	d_qbs_col(:) = (zqbs_col(:,2)-zqbs_col(:,1))/phys_tstep
	635	d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + d_qbs_col(:)
[3927]	636
	637	d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
	638
	639	d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep
	640	d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/phys_tstep
	641	d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
	642	d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
[4523]	643	d_h_qbs_col(:) = (zh_qbs_col(:,2)-zh_qbs_col(:,1))/phys_tstep
[3927]	644
	645	d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep
	646
	647	end if ! end if (fl_ebil .GT. 0)
	648	!
	649	! When in diagnostic mode, restore "out" variables to initial values.
	650	IF (diag_mode == 1) THEN
	651	u_seri(:,:) = sav_u_seri(:,:)
	652	v_seri(:,:) = sav_v_seri(:,:)
	653	ql_seri(:,:) = sav_ql_seri(:,:)
	654	qs_seri(:,:) = sav_qs_seri(:,:)
[4523]	655	qbs_seri(:,:) = sav_qbs_seri(:,:)
[3927]	656	q_seri(:,:) = sav_q_seri(:,:)
	657	t_seri(:,:) = sav_t_seri(:,:)
	658	zdq(:,:) = sav_zdq(:,:)
	659	#ifdef ISO
	660	xtl_seri(:,:,:) = sav_xtl_seri(:,:,:)
	661	xts_seri(:,:,:) = sav_xts_seri(:,:,:)
	662	xt_seri(:,:,:) = sav_xt_seri(:,:,:)
	663	zdxt(:,:,:) = sav_zdxt(:,:,:)
	664	#endif
	665	ENDIF ! (mode == 1)
	666
	667	RETURN
	668	END SUBROUTINE add_phys_tend
	669
[4523]	670	SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, qbs, &
	671	zdu,zdv,zdt,zdq,zdql,zdqs,zdqbs,paprs,text)
[3927]	672	!======================================================================
	673	! Ajoute les tendances des variables physiques aux variables
	674	! d'etat de la dynamique t_seri, q_seri ...
	675	! On en profite pour faire des tests sur les tendances en question.
	676	!======================================================================
	677	! C Risi: on ne met pas les isos dedans car elle n'est jamais apelée.
	678
	679	!======================================================================
	680	! Declarations
	681	!======================================================================
	682
	683	USE phys_state_var_mod, ONLY : phys_tstep, ftsol
	684	USE geometry_mod, ONLY: longitude_deg, latitude_deg
	685	USE print_control_mod, ONLY: prt_level
	686	USE cmp_seri_mod
[4523]	687	USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
	688	& , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
[3927]	689	IMPLICIT none
	690	include "YOMCST.h"
	691	include "clesphys.h"
	692
	693	! Arguments :
	694	!------------
	695	INTEGER, INTENT(IN) :: nlon, nlev
	696	REAL, DIMENSION(nlon,nlev), INTENT(IN) :: uu, vv
[4523]	697	REAL, DIMENSION(nlon,nlev), INTENT(IN) :: temp, qv, ql, qs, qbs
[3927]	698	REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zdu, zdv
[4523]	699	REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zdt, zdq, zdql, zdqs, zdqbs
[3927]	700	REAL, DIMENSION(nlon,nlev+1), INTENT(IN) :: paprs
	701	CHARACTER(), INTENT(IN) :: text
	702
	703	! Local :
	704	!--------
	705	REAL, DIMENSION(nlon,nlev) :: uu_n, vv_n
[4523]	706	REAL, DIMENSION(nlon,nlev) :: temp_n, qv_n, ql_n, qs_n, qbs_n
[3927]	707
	708
	709	!
	710	INTEGER k, n
	711
	712	integer debug_level
	713	logical, save :: first=.true.
	714	!$OMP THREADPRIVATE(first)
	715	!
	716	!======================================================================
	717	! Variables for energy conservation tests
	718	!======================================================================
	719	!
	720
	721	! zh_col------- total enthalpy of vertical air column
	722	! (air with watter vapour, liquid and solid) (J/m2)
	723	! zh_dair_col--- total enthalpy of dry air (J/m2)
	724	! zh_qw_col---- total enthalpy of watter vapour (J/m2)
	725	! zh_ql_col---- total enthalpy of liquid watter (J/m2)
	726	! zh_qs_col---- total enthalpy of solid watter (J/m2)
	727	! zqw_col------ total mass of watter vapour (kg/m2)
	728	! zql_col------ total mass of liquid watter (kg/m2)
[4523]	729	! zqs_col------ total mass of cloud ice (kg/m2)
	730	! zqbs_col------ total mass of blowing snow (kg/m2)
[3927]	731	! zek_col------ total kinetic energy (kg/m2)
	732	!
	733	REAL zairm(nlon, nlev) ! layer air mass (kg/m2)
	734	REAL zqw_col(nlon,2)
	735	REAL zql_col(nlon,2)
	736	REAL zqs_col(nlon,2)
[4523]	737	REAL zqbs_col(nlon,2)
[3927]	738	REAL zek_col(nlon,2)
	739	REAL zh_dair_col(nlon,2)
[4523]	740	REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2), zh_qbs_col(nlon,2)
[3927]	741	REAL zh_col(nlon,2)
	742
	743	!======================================================================
	744	! Initialisations
	745
	746	IF (prt_level >= 5) then
	747	write (,) "In diag_phys_tend, after ",text
	748	call flush
	749	end if
	750
	751	debug_level=10
	752	if (first) then
	753	print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
	754	first=.false.
	755	endif
	756	!
	757	! print *,'add_phys_tend: paprs ',paprs
	758	!======================================================================
	759	! Diagnostics for energy conservation tests
	760	!======================================================================
	761	DO k = 1, nlev
	762	! layer air mass
	763	zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
	764	END DO
	765
	766	if (fl_ebil .GT. 0) then
	767	! ------------------------------------------------
	768	! Compute vertical sum for each atmospheric column
	769	! ------------------------------------------------
	770	n=1 ! begining of time step
	771
	772	CALL integr_v(nlon, nlev, rcpv, &
[4523]	773	temp, qv, ql, qs, qbs, uu, vv, zairm, &
	774	zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
	775	zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
[3927]	776
	777	end if ! end if (fl_ebil .GT. 0)
	778
	779	!======================================================================
	780	! Ajout des tendances sur le vent, la temperature et les diverses phases de l'eau
	781	!======================================================================
	782
	783	uu_n(:,:)=uu(:,:)+zdu(:,:)
	784	vv_n(:,:)=vv(:,:)+zdv(:,:)
	785	qv_n(:,:)=qv(:,:)+zdq(:,:)
	786	ql_n(:,:)=ql(:,:)+zdql(:,:)
	787	qs_n(:,:)=qs(:,:)+zdqs(:,:)
[4523]	788	qbs_n(:,:)=qbs(:,:)+zdqbs(:,:)
[3927]	789	temp_n(:,:)=temp(:,:)+zdt(:,:)
	790
	791
	792
	793	!======================================================================
	794	! Diagnostics for energy conservation tests
	795	!======================================================================
	796
	797	if (fl_ebil .GT. 0) then
	798
	799	! ------------------------------------------------
	800	! Compute vertical sum for each atmospheric column
	801	! ------------------------------------------------
	802	n=2 ! end of time step
	803
	804	CALL integr_v(nlon, nlev, rcpv, &
[4523]	805	temp_n, qv_n, ql_n, qs_n, qbs_n, uu_n, vv_n, zairm, &
	806	zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
	807	zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
[3927]	808
	809	! ------------------------------------------------
	810	! Compute the changes by unit of time
	811	! ------------------------------------------------
	812
	813	d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/phys_tstep
	814	d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
	815	d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
[4523]	816	d_qbs_col(:) = (zqbs_col(:,2)-zqbs_col(:,1))/phys_tstep
	817	d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + d_qbs_col(:)
[3927]	818
	819	d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
	820
	821	print *,'zdu ', zdu
	822	print *,'zdv ', zdv
	823	print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1)
	824
	825	d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep
	826	d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/phys_tstep
	827	d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
	828	d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
[4523]	829	d_h_qbs_col(:) = (zh_qbs_col(:,2)-zh_qbs_col(:,1))/phys_tstep
[3927]	830
	831	d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep
	832
	833	end if ! end if (fl_ebil .GT. 0)
	834	!
	835
	836	RETURN
	837	END SUBROUTINE diag_phys_tend
	838
	839	SUBROUTINE integr_v(nlon, nlev, zcpvap, &
[4523]	840	temp, qv, ql, qs, qbs, uu, vv, zairm, &
	841	zqw_col, zql_col, zqs_col, zqbs_col, zek_col, zh_dair_col, &
	842	zh_qw_col, zh_ql_col, zh_qs_col, zh_qbs_col, zh_col)
[3927]	843
	844	IMPLICIT none
	845	include "YOMCST.h"
	846
	847	INTEGER, INTENT(IN) :: nlon,nlev
	848	REAL, INTENT(IN) :: zcpvap
[4523]	849	REAL, DIMENSION(nlon,nlev), INTENT(IN) :: temp, qv, ql, qs, qbs, uu, vv
[3927]	850	REAL, DIMENSION(nlon,nlev), INTENT(IN) :: zairm
	851	REAL, DIMENSION(nlon), INTENT(OUT) :: zqw_col
	852	REAL, DIMENSION(nlon), INTENT(OUT) :: zql_col
[4523]	853	REAL, DIMENSION(nlon), INTENT(OUT) :: zqs_col, zqbs_col
[3927]	854	REAL, DIMENSION(nlon), INTENT(OUT) :: zek_col
	855	REAL, DIMENSION(nlon), INTENT(OUT) :: zh_dair_col
	856	REAL, DIMENSION(nlon), INTENT(OUT) :: zh_qw_col
	857	REAL, DIMENSION(nlon), INTENT(OUT) :: zh_ql_col
[4523]	858	REAL, DIMENSION(nlon), INTENT(OUT) :: zh_qs_col, zh_qbs_col
[3927]	859	REAL, DIMENSION(nlon), INTENT(OUT) :: zh_col
	860
	861	INTEGER :: i, k
	862
	863
	864	! Reset variables
	865	zqw_col(:) = 0.
	866	zql_col(:) = 0.
	867	zqs_col(:) = 0.
[4523]	868	zqbs_col(:) = 0.
[3927]	869	zek_col(:) = 0.
	870	zh_dair_col(:) = 0.
	871	zh_qw_col(:) = 0.
	872	zh_ql_col(:) = 0.
	873	zh_qs_col(:) = 0.
[4523]	874	zh_qbs_col(:) = 0.
[3927]	875
	876	!JLD write (,) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
	877
	878	! ------------------------------------------------
	879	! Compute vertical sum for each atmospheric column
	880	! ------------------------------------------------
	881	DO k = 1, nlev
	882	DO i = 1, nlon
	883	! Watter mass
	884	zqw_col(i) = zqw_col(i) + qv(i, k)*zairm(i, k)
	885	zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
	886	zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
[4523]	887	zqbs_col(i)= zqbs_col(i) + qbs(i,k)*zairm(i,k)
[3927]	888	! Kinetic Energy
	889	zek_col(i) = zek_col(i) + 0.5(uu(i,k)2+vv(i,k)2)zairm(i, k)
	890	! Air enthalpy : dry air, water vapour, liquid, solid
	891	zh_dair_col(i) = zh_dair_col(i) + rcpd(1.-qv(i,k)-ql(i,k)-qs(i,k)) &
	892	zairm(i, k)*temp(i, k)
	893	zh_qw_col(i) = zh_qw_col(i) + zcpvaptemp(i, k) qv(i, k)*zairm(i, k) !jyg
	894	zh_ql_col(i) = zh_ql_col(i) + (zcpvaptemp(i, k) - rlvtt)ql(i, k)*zairm(i, k) !jyg
	895	zh_qs_col(i) = zh_qs_col(i) + (zcpvaptemp(i, k) - rlstt)qs(i, k)*zairm(i, k) !jyg
[4523]	896	zh_qbs_col(i) = zh_qbs_col(i) + (zcpvaptemp(i, k) - rlstt)qbs(i, k)*zairm(i, k) !jyg
[3927]	897	END DO
	898	END DO
	899	! compute total air enthalpy
[4523]	900	zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:) + zh_qbs_col(:)
[3927]	901
	902	END SUBROUTINE integr_v
	903
	904	SUBROUTINE prt_enerbil (text, itap)
	905	!======================================================================
	906	! Print enenrgy budget diagnotics for the 1D case
	907	!======================================================================
	908
	909	!======================================================================
	910	! Declarations
	911	!======================================================================
	912
	913	USE dimphy, ONLY: klon, klev
	914	USE phys_state_var_mod, ONLY : phys_tstep
[4523]	915	USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con, bs_fall
[3927]	916	USE geometry_mod, ONLY: longitude_deg, latitude_deg
	917	USE print_control_mod, ONLY: prt_level
	918	USE cmp_seri_mod
[4523]	919	USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
	920	& , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
[3927]	921	USE phys_local_var_mod, ONLY: evap, sens
[4523]	922	USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, qbs_seri, q_seri, t_seri &
[3927]	923	& , rain_lsc, snow_lsc
	924	USE climb_hq_mod, ONLY : d_h_col_vdf, f_h_bnd
	925	IMPLICIT none
	926	include "YOMCST.h"
	927
	928	! Arguments :
	929	!------------
	930	CHARACTER() text ! text specifing the involved parametrization
	931	integer itap ! time step number
	932	! local variables
	933	! ---------------
	934	real bilq_seuil, bilh_seuil ! thresold on error in Q and H budget
	935	real bilq_error, bilh_error ! erros in Q and H budget
	936	real bilq_bnd, bilh_bnd ! Q and H budget due to exchange with boundaries
	937	integer bilq_ok, bilh_ok
	938	CHARACTER*(12) status
	939
	940	bilq_seuil = 1.E-10
	941	bilh_seuil = 1.E-1
	942	bilq_ok=0
	943	bilh_ok=0
	944
	945	!!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil
	946	if ( (fl_ebil .GT. 0) .and. (klon .EQ. 1)) then
	947
	948	bilq_bnd = 0.
	949	bilh_bnd = 0.
	950
	951	param: SELECT CASE (text)
	952	CASE("vdf") param
	953	bilq_bnd = evap(1)
	954	bilh_bnd = sens(1)+(rcpv-rcpd)evap(1)t_seri(1,1)
	955	CASE("lsc") param
	956	bilq_bnd = - rain_lsc(1) - snow_lsc(1)
	957	bilh_bnd = (-(rcw-rcpd)t_seri(1,1) + rlvtt) rain_lsc(1) &
	958	& + (-(rcs-rcpd)t_seri(1,1) + rlstt) snow_lsc(1)
[4523]	959	CASE("bs") param
	960	bilq_bnd = - bs_fall(1)
	961	bilh_bnd = (-(rcs-rcpd)t_seri(1,1) + rlstt) bs_fall(1)
[3927]	962	CASE("convection") param
	963	bilq_bnd = - rain_con(1) - snow_con(1)
	964	bilh_bnd = (-(rcw-rcpd)t_seri(1,1) + rlvtt) rain_con(1) &
	965	& + (-(rcs-rcpd)t_seri(1,1) + rlstt) snow_con(1)
	966	CASE("SW") param
	967	bilh_bnd = topsw(1) - solsw(1)
	968	CASE("LW") param
	969	bilh_bnd = -(toplw(1) + sollw(1))
	970	CASE DEFAULT param
	971	bilq_bnd = 0.
	972	bilh_bnd = 0.
	973	END SELECT param
	974
	975	bilq_error = d_qt_col(1) - bilq_bnd
	976	bilh_error = d_h_col(1) - bilh_bnd
	977	! are the errors too large?
	978	if ( abs(bilq_error) .gt. bilq_seuil) bilq_ok=1
	979	if ( abs(bilh_error) .gt. bilh_seuil) bilh_ok=1
	980	!
	981	! Print diagnostics
	982	! =================
	983	if ( (bilq_ok .eq. 0).and.(bilh_ok .eq. 0) ) then
	984	status="enerbil-OK"
	985	else
	986	status="enerbil-PB"
	987	end if
	988
	989	if ( prt_level .GE. 3) then
	990	write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error," H", bilh_error
	991	9010 format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6)
	992	end if
	993	if ( prt_level .GE. 3) then
	994	write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1)
	995	end if
	996	if ( prt_level .GE. 5) then
	997	write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
[4523]	998	write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1)
	999	write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1),d_h_qbs_col(1)
[3927]	1000	end if
	1001
	1002	specific_diag: SELECT CASE (text)
	1003	CASE("vdf") specific_diag
	1004	if ( prt_level .GE. 5) then
	1005	write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1)
	1006	write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1)
	1007	end if
	1008	CASE("lsc") specific_diag
	1009	if ( prt_level .GE. 5) then
	1010	write(,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt rain_lsc(1), -(rcw-rcpd)t_seri(1,1) rain_lsc(1)
	1011	write(,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt snow_lsc(1), -(rcs-rcpd)t_seri(1,1) snow_lsc(1)
	1012	end if
	1013	CASE("convection") specific_diag
	1014	if ( prt_level .GE. 5) then
	1015	write(,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt rain_con(1), -(rcw-rcpd)t_seri(1,1) rain_con(1)
	1016	write(,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt snow_con(1), -(rcs-rcpd)t_seri(1,1) snow_con(1)
	1017	end if
	1018	END SELECT specific_diag
	1019
	1020	9000 format (1x,A8,2x,A35,10E15.6)
	1021
	1022	end if ! end if (fl_ebil .GT. 0)
	1023
	1024	END SUBROUTINE prt_enerbil
	1025
	1026	END MODULE add_phys_tend_mod

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

Download in other formats:

Original Format