Context Navigation

climb_hq_mod.F90 @ 4087

Last change on this file since 4087 was 4004, checked in by Laurent Fairhead, 3 years ago
Added properties on phylmdiso directory LF
Property svn:keywords set to `Id`
File size: 28.2 KB

Rev	Line
[3927]	1	MODULE climb_hq_mod
	2	!
	3	! Module to solve the verctical diffusion of "q" and "H";
	4	! specific humidity and potential energi.
	5	!
	6	USE dimphy
	7	#ifdef ISO
	8	USE infotrac_phy, ONLY: ntraciso ! ajout C Risi pour isos
	9	#endif
	10
	11	IMPLICIT NONE
	12	SAVE
	13	PRIVATE
	14	PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
	15
	16	REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
	17	!$OMP THREADPRIVATE(gamaq,gamah)
	18	REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
	19	!$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
	20	REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
	21	!$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
	22	REAL, DIMENSION(:), ALLOCATABLE :: Acoef_Q, Bcoef_Q
	23	!$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
	24	REAL, DIMENSION(:), ALLOCATABLE :: Acoef_H, Bcoef_H
	25	!$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
	26	REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
	27	!$OMP THREADPRIVATE(Kcoefhq)
	28	REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
	29	!$OMP THREADPRIVATE(h_old)
	30	REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
	31	!$OMP THREADPRIVATE(d_h_col_vdf)
	32	REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
	33	!$OMP THREADPRIVATE(f_h_bnd)
	34	#ifdef ISO
	35	REAL, DIMENSION(:,:,:), ALLOCATABLE :: gamaxt
	36	!$OMP THREADPRIVATE(gamaxt)
	37	REAL, DIMENSION(:,:,:), ALLOCATABLE :: Ccoef_XT, Dcoef_XT
	38	!$OMP THREADPRIVATE(Ccoef_XT, Dcoef_XT)
	39	REAL, DIMENSION(:,:), ALLOCATABLE :: Acoef_XT, Bcoef_XT
	40	!$OMP THREADPRIVATE(Acoef_XT, Bcoef_XT)
	41	#endif
	42
	43	CONTAINS
	44	!
	45	!****************************************************************************************
	46	!
	47	SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
	48	delp, temp, q, dtime, &
	49	!!! nrlmd le 02/05/2011
	50	Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
	51	Kcoef_hq_out, gama_q_out, gama_h_out, &
	52	!!!
	53	Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out &
	54	#ifdef ISO
	55	,xt, &
	56	Ccoef_XT_out,Dcoef_XT_out,gama_xt_out, &
	57	Acoef_XT_out, Bcoef_XT_out &
	58	#endif
	59	)
	60	#ifdef ISOVERIF
	61	!USE infotrac_phy, ONLY: use_iso
	62	USE isotopes_mod, ONLY: iso_eau,iso_HDO
	63	!USE isotopes_verif_mod, ONLY: errmax, errmaxrel
	64	USE isotopes_verif_mod
	65	#endif
	66
	67	! This routine calculates recursivly the coefficients C and D
	68	! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
	69	! the index of the vertical layer.
	70	!
	71	! Input arguments
	72	!****************************************************************************************
	73	INTEGER, INTENT(IN) :: knon
	74	REAL, DIMENSION(klon,klev), INTENT(IN) :: coefhq
	75	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
	76	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
	77	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp, delp ! temperature
	78	REAL, DIMENSION(klon,klev), INTENT(IN) :: q
	79	REAL, INTENT(IN) :: dtime
	80	#ifdef ISO
	81	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt
	82	#endif
	83
	84
	85	! Output arguments
	86	!****************************************************************************************
	87	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out
	88	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out
	89	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out
	90	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out
	91	#ifdef ISO
	92	REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: Acoef_XT_out
	93	REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: Bcoef_XT_out
	94	#endif
	95
	96	!!! nrlmd le 02/05/2011
	97	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_H_out
	98	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_Q_out
	99	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_H_out
	100	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_Q_out
	101	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_hq_out
	102	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_q_out
	103	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_h_out
	104	#ifdef ISO
	105	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: Ccoef_XT_out
	106	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: Dcoef_XT_out
	107	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: gama_xt_out
	108	#endif
	109	!!!
	110
	111	! Local variables
	112	!****************************************************************************************
	113	LOGICAL, SAVE :: first=.TRUE.
	114	!$OMP THREADPRIVATE(first)
	115	! JLD now renamed h_old and declared in module
	116	! REAL, DIMENSION(klon,klev) :: local_H
	117	REAL, DIMENSION(klon) :: psref
	118	REAL :: delz, pkh
	119	INTEGER :: k, i, ierr
	120
	121	#ifdef ISO
	122	real, DIMENSION(klon,2:klev) :: gamaxt_tmp
	123	real, DIMENSION(klon,klev) :: xt_tmp
	124	real, DIMENSION(klon,klev) :: Ccoef_XT_tmp,Dcoef_XT_tmp
	125	real, DIMENSION(klon) :: Acoef_XT_tmp,Bcoef_XT_tmp
	126	integer ixt
	127	#endif
	128	! Include
	129	!****************************************************************************************
	130	INCLUDE "YOMCST.h"
	131	INCLUDE "compbl.h"
	132
	133	#ifdef ISO
	134	#ifdef ISOVERIF
	135	if (iso_eau.gt.0) then
	136	do k = 1, klev
	137	DO i = 1, knon
	138	call iso_verif_egalite_choix( &
	139	xt(iso_eau,i,k),q(i,k), &
	140	'climb_hq 100',errmax,errmaxrel)
	141	enddo
	142	enddo
	143	endif ! if (iso_eau.gt.0) then
	144	#endif
	145	#endif
	146
	147	!****************************************************************************************
	148	! 1)
	149	! Allocation at first time step only
	150	!
	151	!****************************************************************************************
	152
	153	IF (first) THEN
	154	first=.FALSE.
	155	ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
	156	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
	157
	158	ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
	159	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
	160
	161	ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
	162	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
	163
	164	ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
	165	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
	166
	167	ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
	168	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
	169
	170	ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
	171	IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
	172
	173	ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
	174	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
	175
	176	ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
	177	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
	178
	179	#ifdef ISO
	180	ALLOCATE(Ccoef_XT(ntraciso,klon,klev), STAT=ierr)
	181	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_XT, ierr=', ierr
	182
	183	ALLOCATE(Dcoef_XT(ntraciso,klon,klev), STAT=ierr)
	184	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_XT, ierr=', ierr
	185
	186	ALLOCATE(Acoef_XT(ntraciso,klon), Bcoef_XT(ntraciso,klon), STAT=ierr)
	187	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_XT and Bcoef_XT, ierr=', ierr
	188
	189	ALLOCATE(gamaxt(ntraciso,1:klon,2:klev), STAT=ierr)
	190	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaxt, ierr=', ierr
	191	#endif
	192
	193	ALLOCATE(h_old(klon,klev), STAT=ierr)
	194	IF ( ierr /= 0 ) PRINT*,' pb in allloc h_old, ierr=', ierr
	195
	196	ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
	197	IF ( ierr /= 0 ) PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
	198
	199	ALLOCATE(f_h_bnd(klon), STAT=ierr)
	200	IF ( ierr /= 0 ) PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
	201	END IF
	202
	203	!****************************************************************************************
	204	! 2)
	205	! Definition of the coeficient K
	206	!
	207	!****************************************************************************************
	208	Kcoefhq(:,:) = 0.0
	209	DO k = 2, klev
	210	DO i = 1, knon
	211	Kcoefhq(i,k) = &
	212	coefhq(i,k)RGRG*dtime /(pplay(i,k-1)-pplay(i,k)) &
	213	(paprs(i,k)2/(temp(i,k)+temp(i,k-1))/RD)**2
	214	ENDDO
	215	ENDDO
	216
	217	!****************************************************************************************
	218	! 3)
	219	! Calculation of gama for "Q" and "H"
	220	!
	221	!****************************************************************************************
	222	! surface pressure is used as reference
	223	psref(:) = paprs(:,1)
	224
	225	! definition of gama
	226	IF (iflag_pbl == 1) THEN
	227	gamaq(:,:) = 0.0
	228	gamah(:,:) = -1.0e-03
	229	gamah(:,2) = -2.5e-03
	230	#ifdef ISO
	231	do ixt=1,ntraciso
	232	gamaxt(:,:,:) = 0.0
	233	enddo ! do ixt=1,ntraciso
	234	#endif
	235
	236	! conversion de gama
	237	DO k = 2, klev
	238	DO i = 1, knon
	239	delz = RD * (temp(i,k-1)+temp(i,k)) / &
	240	2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
	241	pkh = (psref(i)/paprs(i,k))**RKAPPA
	242
	243	! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
	244	gamaq(i,k) = gamaq(i,k) * delz
	245	! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
	246	gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
	247	#ifdef ISO
	248	do ixt=1,ntraciso
	249	gamaxt(ixt,i,k) = gamaxt(ixt,i,k) * delz
	250	enddo
	251	#endif
	252	ENDDO
	253	ENDDO
	254
	255	ELSE
	256	gamaq(:,:) = 0.0
	257	gamah(:,:) = 0.0
	258	#ifdef ISO
	259	do ixt=1,ntraciso
	260	gamaxt(:,:,:) = 0.0
	261	enddo ! do ixt=1,ntraciso
	262	#endif
	263	ENDIF
	264
	265	#ifdef ISO
	266	#ifdef ISOVERIF
	267	do k = 2, klev
	268	DO i = 1, knon
	269	call iso_verif_egalite_choix( &
	270	gamaxt(iso_eau,i,k),gamaq(i,k), &
	271	'climb_hq 209',errmax,errmaxrel)
	272	enddo
	273	enddo
	274	#endif
	275	#endif
	276
	277
	278	!****************************************************************************************
	279	! 4)
	280	! Calculte the coefficients C and D for specific humidity, q
	281	!
	282	!****************************************************************************************
	283
	284	CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
	285	Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
	286
	287
	288	#ifdef ISO
	289	do ixt=1,ntraciso
	290	! compression
	291	do k = 2, klev
	292	DO i = 1, knon
	293	gamaxt_tmp(i,k)=gamaxt(ixt,i,k)
	294	enddo
	295	enddo !do k = 2, klev
	296	do k = 1, klev
	297	DO i = 1, knon
	298	xt_tmp(i,k)=xt(ixt,i,k)
	299	enddo
	300	enddo !do k = 2, klev
	301	!appel routine generique
	302	CALL calc_coef(knon, Kcoefhq(:,:), gamaxt_tmp(:,:), delp(:,:), xt_tmp(:,:), &
	303	Ccoef_XT_tmp(:,:), Dcoef_XT_tmp(:,:), Acoef_XT_tmp, Bcoef_XT_tmp)
	304
	305	! decompression
	306	do k = 1, klev
	307	DO i = 1, knon
	308	Ccoef_XT(ixt,i,k)=Ccoef_XT_tmp(i,k)
	309	Dcoef_XT(ixt,i,k)=Dcoef_XT_tmp(i,k)
	310	enddo
	311	enddo !do k = 2, klev
	312	DO i = 1, knon
	313	Acoef_XT(ixt,i)=Acoef_XT_tmp(i)
	314	Bcoef_XT(ixt,i)=Bcoef_XT_tmp(i)
	315	enddo
	316	enddo ! do ixt=1,ntraciso
	317	#ifdef ISOVERIF
	318	if (iso_eau.gt.0) then
	319	do k = 1, klev
	320	DO i = 1, knon
	321	call iso_verif_egalite_choix( &
	322	Ccoef_XT(iso_eau,i,k),Ccoef_Q(i,k), &
	323	'climb_hq 234c',errmax,errmaxrel)
	324	call iso_verif_egalite_choix( &
	325	Dcoef_XT(iso_eau,i,k),Dcoef_Q(i,k), &
	326	'climb_hq 234d',errmax,errmaxrel)
	327	enddo !DO i = 1, knon
	328	enddo !do k = 2, klev
	329	DO i = 1, knon
	330	call iso_verif_egalite_choix( &
	331	Acoef_XT(iso_eau,i),Acoef_Q(i), &
	332	'climb_hq 234a',errmax,errmaxrel)
	333	call iso_verif_egalite_choix( &
	334	Bcoef_XT(iso_eau,i),Bcoef_Q(i), &
	335	'climb_hq 234b',errmax,errmaxrel)
	336	enddo !DO i = 1, knon
	337	endif !if (iso_eau.gt.0) then
	338	#endif
	339	#endif
	340	!****************************************************************************************
	341	! 5)
	342	! Calculte the coefficients C and D for potentiel entalpie, H
	343	!
	344	!****************************************************************************************
	345	h_old(:,:) = 0.0
	346
	347	DO k=1,klev
	348	DO i = 1, knon
	349	! convertie la temperature en entalpie potentielle
	350	h_old(i,k) = RCPD * temp(i,k) * &
	351	(psref(i)/pplay(i,k))**RKAPPA
	352	ENDDO
	353	ENDDO
	354
	355	CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
	356	Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
	357
	358	!****************************************************************************************
	359	! 6)
	360	! Return the first layer in output variables
	361	!
	362	!****************************************************************************************
	363	Acoef_H_out = Acoef_H
	364	Bcoef_H_out = Bcoef_H
	365	Acoef_Q_out = Acoef_Q
	366	Bcoef_Q_out = Bcoef_Q
	367	#ifdef ISO
	368	Acoef_XT_out = Acoef_XT
	369	Bcoef_XT_out = Bcoef_XT
	370	#endif
	371
	372	!****************************************************************************************
	373	! 7)
	374	! If Pbl is split, return also the other layers in output variables
	375	!
	376	!****************************************************************************************
	377	!!! jyg le 07/02/2012
	378	!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
	379	IF (mod(iflag_pbl_split,10) .ge.1) THEN
	380	!!! nrlmd le 02/05/2011
	381	DO k= 1, klev
	382	DO i= 1, klon
	383	Ccoef_H_out(i,k) = Ccoef_H(i,k)
	384	Dcoef_H_out(i,k) = Dcoef_H(i,k)
	385	Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
	386	Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
	387	Kcoef_hq_out(i,k) = Kcoefhq(i,k)
	388	#ifdef ISO
	389	do ixt=1,ntraciso
	390	Ccoef_XT_out(ixt,i,k) = Ccoef_XT(ixt,i,k)
	391	Dcoef_XT_out(ixt,i,k) = Dcoef_XT(ixt,i,k)
	392	enddo
	393	#endif
	394	IF (k.eq.1) THEN
	395	gama_h_out(i,k) = 0.
	396	gama_q_out(i,k) = 0.
	397	#ifdef ISO
	398	do ixt=1,ntraciso
	399	gama_xt_out(ixt,i,k) = 0.
	400	enddo
	401	#endif
	402	ELSE
	403	gama_h_out(i,k) = gamah(i,k)
	404	gama_q_out(i,k) = gamaq(i,k)
	405	#ifdef ISO
	406	do ixt=1,ntraciso
	407	gama_xt_out(ixt,i,k) = gamaxt(ixt,i,k)
	408	enddo
	409	#endif
	410	ENDIF
	411	ENDDO
	412	ENDDO
	413	!!!
	414	ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
	415	!!!
	416
	417	END SUBROUTINE climb_hq_down
	418	!
	419	!****************************************************************************************
	420	!
	421	SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
	422	!
	423	! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
	424	! where X is H or Q, and k the vertical level k=1,klev
	425	!
	426	INCLUDE "YOMCST.h"
	427	! Input arguments
	428	!****************************************************************************************
	429	INTEGER, INTENT(IN) :: knon
	430	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp
	431	REAL, DIMENSION(klon,klev), INTENT(IN) :: X
	432	REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
	433
	434	! Output arguments
	435	!****************************************************************************************
	436	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
	437	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef
	438
	439	! Local variables
	440	!****************************************************************************************
	441	INTEGER :: k, i
	442	REAL :: buf
	443
	444	!****************************************************************************************
	445	! Niveau au sommet, k=klev
	446	!
	447	!****************************************************************************************
	448	Ccoef(:,:) = 0.0
	449	Dcoef(:,:) = 0.0
	450
	451	DO i = 1, knon
	452	buf = delp(i,klev) + Kcoef(i,klev)
	453
	454	Ccoef(i,klev) = (X(i,klev)delp(i,klev) - Kcoef(i,klev)gama(i,klev))/buf
	455	Dcoef(i,klev) = Kcoef(i,klev)/buf
	456	END DO
	457
	458
	459	!****************************************************************************************
	460	! Niveau (klev-1) <= k <= 2
	461	!
	462	!****************************************************************************************
	463
	464	DO k=(klev-1),2,-1
	465	DO i = 1, knon
	466	buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
	467	Ccoef(i,k) = (X(i,k)delp(i,k) + Kcoef(i,k+1)Ccoef(i,k+1) + &
	468	Kcoef(i,k+1)gama(i,k+1) - Kcoef(i,k)gama(i,k))/buf
	469	Dcoef(i,k) = Kcoef(i,k)/buf
	470	END DO
	471	END DO
	472
	473	!****************************************************************************************
	474	! Niveau k=1
	475	!
	476	!****************************************************************************************
	477
	478	DO i = 1, knon
	479	buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
	480	Acoef(i) = (X(i,1)delp(i,1) + Kcoef(i,2)(gama(i,2)+Ccoef(i,2)))/buf
	481	Bcoef(i) = -1. * RG / buf
	482	END DO
	483
	484	END SUBROUTINE calc_coef
	485	!
	486	!****************************************************************************************
	487	!
	488	SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
	489	flx_q1, flx_h1, paprs, pplay, &
	490	!!! nrlmd le 02/05/2011
	491	Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
	492	Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
	493	Kcoef_hq_in, gama_q_in, gama_h_in, &
	494	!!!
	495	flux_q, flux_h, d_q, d_t &
	496	#ifdef ISO
	497	,xt_old, flx_xt1, &
	498	Acoef_XT_in,Bcoef_XT_in,Ccoef_XT_in,Dcoef_XT_in,gama_xt_in, &
	499	flux_xt, d_xt &
	500	#endif
	501	)
	502
	503	#ifdef ISOVERIF
	504	USE infotrac_phy, ONLY: ok_isotrac
	505	USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18, ridicule
	506	USE isotopes_verif_mod
	507	#endif
	508	!
	509	! This routine calculates the flux and tendency of the specific humidity q and
	510	! the potential engergi H.
	511	! The quantities q and H are calculated according to
	512	! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
	513	! C and D are known from before and k is index of the vertical layer.
	514	!
	515
	516	! Input arguments
	517	!****************************************************************************************
	518	INTEGER, INTENT(IN) :: knon
	519	REAL, INTENT(IN) :: dtime
	520	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_old, q_old
	521	REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1
	522	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
	523	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
	524
	525	!!! nrlmd le 02/05/2011
	526	REAL, DIMENSION(klon), INTENT(IN) :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
	527	REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
	528	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_hq_in, gama_q_in, gama_h_in
	529	#ifdef ISO
	530	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: Acoef_XT_in, Bcoef_XT_in
	531	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: Ccoef_XT_in, Dcoef_XT_in
	532	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: gama_xt_in
	533	#endif
	534	!!!
	535
	536	! Output arguments
	537	!****************************************************************************************
	538	REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t
	539
	540	! Local variables
	541	!****************************************************************************************
	542	LOGICAL, SAVE :: last=.FALSE.
	543	REAL, DIMENSION(klon,klev) :: h_new, q_new
	544	REAL, DIMENSION(klon) :: psref
	545	INTEGER :: k, i, ierr
	546
	547	#ifdef ISO
	548	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt_old
	549	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: flux_xt, d_xt
	550	REAL, DIMENSION(ntraciso,klon,klev) :: xt_new
	551	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: flx_xt1
	552	integer ixt
	553	!#ifdef ISOVERIF
	554	! integer iso_verif_noNaN_nostop
	555	!#endif
	556	#endif
	557
	558	! Include
	559	!****************************************************************************************
	560	INCLUDE "YOMCST.h"
	561	INCLUDE "compbl.h"
	562
	563	!****************************************************************************************
	564	! 1)
	565	! Definition of some variables
	566	REAL, DIMENSION(klon,klev) :: d_h, zairm
	567	!
	568	!****************************************************************************************
	569
	570
	571	#ifdef ISO
	572	#ifdef ISOVERIF
	573	DO k = 1, klev
	574	DO i = 1, knon
	575	if (iso_eau.gt.0) then
	576	call iso_verif_egalite(xt_old(iso_eau,i,k), &
	577	& q_old(i,k),'climb_hq_mod 421')
	578	endif
	579	#ifdef ISOTRAC
	580	if (ok_isotrac) then
	581	call iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 422')
	582	endif
	583	#endif
	584	enddo
	585	enddo
	586	#endif
	587	#endif
	588
	589	flux_q(:,:) = 0.0
	590	flux_h(:,:) = 0.0
	591	d_q(:,:) = 0.0
	592	d_t(:,:) = 0.0
	593	d_h(:,:) = 0.0
	594	f_h_bnd(:)= 0.0
	595	#ifdef ISO
	596	q_new(:,:) = 0.0
	597	flux_xt(:,:,:) = 0.0
	598	d_xt(:,:,:) = 0.0
	599	xt_new(:,:,:) = 0.0
	600	#endif
	601	psref(1:knon) = paprs(1:knon,1)
	602
	603	!!! jyg le 07/02/2012
	604	!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
	605	IF (mod(iflag_pbl_split,10) .ge.1) THEN
	606	!!! nrlmd le 02/05/2011
	607	DO i = 1, knon
	608	Acoef_H(i)=Acoef_H_in(i)
	609	Acoef_Q(i)=Acoef_Q_in(i)
	610	Bcoef_H(i)=Bcoef_H_in(i)
	611	Bcoef_Q(i)=Bcoef_Q_in(i)
	612	#ifdef ISO
	613	do ixt=1,ntraciso
	614	Acoef_XT(ixt,i)=Acoef_XT_in(ixt,i)
	615	Bcoef_XT(ixt,i)=Bcoef_XT_in(ixt,i)
	616	enddo
	617	#endif
	618	ENDDO
	619	DO k = 1, klev
	620	DO i = 1, knon
	621	Ccoef_H(i,k)=Ccoef_H_in(i,k)
	622	Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
	623	Dcoef_H(i,k)=Dcoef_H_in(i,k)
	624	Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
	625	Kcoefhq(i,k)=Kcoef_hq_in(i,k)
	626	#ifdef ISO
	627	do ixt=1,ntraciso
	628	Ccoef_XT(ixt,i,k)=Ccoef_XT_in(ixt,i,k)
	629	Dcoef_XT(ixt,i,k)=Dcoef_XT_in(ixt,i,k)
	630	enddo
	631	#endif
	632	IF (k.gt.1) THEN
	633	gamah(i,k)=gama_h_in(i,k)
	634	gamaq(i,k)=gama_q_in(i,k)
	635	#ifdef ISO
	636	do ixt=1,ntraciso
	637	gamaxt(ixt,i,k)=gama_xt_in(ixt,i,k)
	638	enddo
	639	#endif
	640	ENDIF
	641	ENDDO
	642	ENDDO
	643	!!!
	644	ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
	645	!!!
	646
	647	!****************************************************************************************
	648	! 2)
	649	! Calculation of Q and H
	650	!
	651	!****************************************************************************************
	652
	653	!- First layer
	654	q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)flx_q1(1:knon)dtime
	655	h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)flx_h1(1:knon)dtime
	656	f_h_bnd(1:knon) = flx_h1(1:knon)
	657	#ifdef ISO
	658	do ixt=1,ntraciso
	659	xt_new(ixt,1:knon,1) = Acoef_XT(ixt,1:knon) + Bcoef_XT(ixt,1:knon)flx_xt1(ixt,1:knon)dtime
	660	enddo ! do ixt=1,ntraciso
	661	#endif
	662	!- All the other layers
	663	DO k = 2, klev
	664	DO i = 1, knon
	665	q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
	666	h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
	667	#ifdef ISO
	668	do ixt=1,ntraciso
	669	xt_new(ixt,i,k) = Ccoef_XT(ixt,i,k) + Dcoef_XT(ixt,i,k)*xt_new(ixt,i,k-1)
	670	enddo ! do ixt=1,ntraciso
	671	#endif
	672	END DO
	673	END DO
	674
	675	#ifdef ISO
	676	#ifdef ISOVERIF
	677	DO k = 1, klev
	678	DO i = 1, knon
	679	do ixt=1,ntraciso
	680	if (iso_verif_noNaN_nostop(xt_new(ixt,i,k),'climb_hq 507').eq.1) then
	681	write(,) 'Acoef_XT(ixt,i)=',Acoef_XT(ixt,i)
	682	write(,) 'Bcoef_XT(ixt,i)=',Bcoef_XT(ixt,i)
	683	write(,) 'flx_xt1(ixt,i)=',flx_xt1(ixt,i)
	684	if (k.ge.2) then
	685	write(,) 'Ccoef_XT(ixt,i,k)=',Ccoef_XT(ixt,i,k)
	686	write(,) 'Dcoef_XT(ixt,i,k)=',Dcoef_XT(ixt,i,k)
	687	endif
	688	stop
	689	endif
	690	enddo !do ixt=1,ntraciso
	691	END DO
	692	END DO
	693	#endif
	694	#ifdef ISOVERIF
	695	if (iso_eau.gt.0) then
	696	call iso_verif_egalite_vect2D( &
	697	& xt_new,q_new, &
	698	& 'climb_hq_mod 504',ntraciso,klon,klev)
	699	endif !if (iso_eau.gt.0) then
	700	if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
	701	do i=1,klon
	702	do k=1,klev
	703	if (q_new(i,k).gt.ridicule) then
	704	if (iso_verif_o18_aberrant_nostop( &
	705	& xt_new(iso_HDO,i,k)/q_new(i,k), &
	706	& xt_new(iso_O18,i,k)/q_new(i,k), &
	707	& 'climb_hq_mod 690').eq.1) then
	708	write(,) 'i,k,q_new(i,k)=',i,k,q_new(i,k)
	709	stop
	710	endif ! if (iso_verif_o18_aberrant_nostop
	711	endif !if (q_seri(i,k).gt.errmax) then
	712	enddo !k=1,klev
	713	enddo !i=1,klon
	714	endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
	715	#endif
	716	#endif
	717	!****************************************************************************************
	718	! 3)
	719	! Calculation of the flux for Q and H
	720	!
	721	!****************************************************************************************
	722
	723	!- The flux at first layer, k=1
	724	flux_q(1:knon,1)=flx_q1(1:knon)
	725	flux_h(1:knon,1)=flx_h1(1:knon)
	726	#ifdef ISO
	727	do ixt=1,ntraciso
	728	flux_xt(ixt,1:knon,1)=flx_xt1(ixt,1:knon)
	729	enddo ! do ixt=1,ntraciso
	730	#endif
	731
	732	!- The flux at all layers above surface
	733	DO k = 2, klev
	734	DO i = 1, knon
	735	flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
	736	(q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
	737
	738	flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
	739	(h_new(i,k)-h_new(i,k-1)+gamah(i,k))
	740
	741	#ifdef ISO
	742	do ixt=1,ntraciso
	743	flux_xt(ixt,i,k) = (Kcoefhq(i,k)/RG/dtime) * &
	744	(xt_new(ixt,i,k)-xt_new(ixt,i,k-1)+gamaxt(ixt,i,k))
	745	enddo ! do ixt=1,ntraciso
	746	#endif
	747	END DO
	748	END DO
	749
	750	!****************************************************************************************
	751	! 4)
	752	! Calculation of tendency for Q and H
	753	!
	754	!****************************************************************************************
	755	d_h_col_vdf(:) = 0.0
	756	DO k = 1, klev
	757	DO i = 1, knon
	758	d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
	759	d_q(i,k) = q_new(i,k) - q_old(i,k)
	760	d_h(i,k) = h_new(i,k) - h_old(i,k)
	761	!JLD d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD !correction a venir
	762	! layer air mass
	763	zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
	764	d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
	765	#ifdef ISO
	766	do ixt=1,ntraciso
	767	d_xt(ixt,i,k) = xt_new(ixt,i,k) - xt_old(ixt,i,k)
	768	enddo ! do ixt=1,ntraciso
	769	#ifdef ISOVERIF
	770	do ixt=1,ntraciso
	771	call iso_verif_noNaN(xt_new(ixt,i,k),'climb_hq 562')
	772	call iso_verif_noNaN(d_xt(ixt,i,k),'climb_hq 562')
	773	enddo ! do ixt=1,ntraciso
	774	#endif
	775	#ifdef ISOVERIF
	776	if (iso_eau.gt.0) then
	777	call iso_verif_egalite(d_xt(iso_eau,i,k), &
	778	& d_q(i,k),'climb_hq_mod 503')
	779	call iso_verif_egalite(xt_new(iso_eau,i,k), &
	780	& q_new(i,k),'climb_hq_mod 503b')
	781	endif
	782	#ifdef ISOTRAC
	783	if (ok_isotrac) then
	784	call iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 526')
	785	endif
	786	#endif
	787	#endif
	788	#endif
	789	END DO
	790	END DO
	791
	792	#ifdef ISO
	793	#ifdef ISOVERIF
	794	! write(,) 'climb_hq_mod 758: d_xt,d_q=',d_xt(iso_eau,1,1),d_q(1,1)
	795	if (iso_eau.gt.0) then
	796	call iso_verif_egalite_vect2D( &
	797	d_xt,d_q, &
	798	'climb_hq_mod 761',ntraciso,klon,klev)
	799	endif
	800	#endif
	801	#endif
	802
	803	!****************************************************************************************
	804	! Some deallocations
	805	!
	806	!****************************************************************************************
	807	IF (last) THEN
	808	DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
	809	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
	810	DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
	811	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
	812	DEALLOCATE(gamaq, gamah,stat=ierr)
	813	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
	814	DEALLOCATE(Kcoefhq,stat=ierr)
	815	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
	816	DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
	817	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
	818
	819	#ifdef ISO
	820	DEALLOCATE(Ccoef_XT, Dcoef_XT ,stat=ierr)
	821	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_XT, Dcoef_XT, ierr=', ierr
	822	DEALLOCATE(Acoef_XT, Bcoef_XT, stat=ierr)
	823	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_XT, Bcoef_XT, ierr=', ierr
	824	DEALLOCATE(gamaxt, stat=ierr)
	825	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaxt, ierr=', ierr
	826	#endif
	827	END IF
	828	END SUBROUTINE climb_hq_up
	829	!
	830	!****************************************************************************************
	831	!
	832	END MODULE climb_hq_mod
	833
	834
	835
	836
	837
	838

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

Download in other formats:

Original Format