Context Navigation

climb_hq_mod.F90 @ 5452

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

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

Download in other formats:

Original Format