Context Navigation

add_phys_tend_mod.F90 @ 5454

Last change on this file since 5454 was 5285, checked in by abarral, 2 months ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
Property svn:keywords set to `Id`
File size: 36.0 KB

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

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

Download in other formats:

Original Format