Context Navigation

climb_hq_mod.F90

Last change on this file was 3331, checked in by acozic, 7 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 26.8 KB

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

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

Download in other formats:

Original Format