Context Navigation

climb_hq_mod.F90 @ 1118

Last change on this file since 1118 was 1085, checked in by jghattas, 16 years ago
Bug fix pour OpenMP. Anne Cozic
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 13.6 KB

Rev	Line
[782]	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
	8	IMPLICIT NONE
	9	SAVE
	10	PRIVATE
	11	PUBLIC :: climb_hq_down, climb_hq_up
	12
	13	REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
	14	!$OMP THREADPRIVATE(gamaq,gamah)
	15	REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
	16	!$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
	17	REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
	18	!$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
[1067]	19	REAL, DIMENSION(:), ALLOCATABLE :: Acoef_Q, Bcoef_Q
	20	!$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
	21	REAL, DIMENSION(:), ALLOCATABLE :: Acoef_H, Bcoef_H
	22	!$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
[782]	23	REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
	24	!$OMP THREADPRIVATE(Kcoefhq)
	25
	26	CONTAINS
	27	!
	28	!****************************************************************************************
	29	!
	30	SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
	31	delp, temp, q, dtime, &
[1067]	32	Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)
[782]	33
[793]	34	INCLUDE "YOMCST.h"
[782]	35	! This routine calculates recursivly the coefficients C and D
	36	! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
	37	! the index of the vertical layer.
	38	!
	39	! Input arguments
	40	!****************************************************************************************
	41	INTEGER, INTENT(IN) :: knon
	42	REAL, DIMENSION(klon,klev), INTENT(IN) :: coefhq
	43	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
	44	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
	45	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp, delp ! temperature
	46	REAL, DIMENSION(klon,klev), INTENT(IN) :: q
	47	REAL, INTENT(IN) :: dtime
	48
	49	! Output arguments
	50	!****************************************************************************************
[1067]	51	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out
	52	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out
	53	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out
	54	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out
[782]	55
	56	! Local variables
	57	!****************************************************************************************
[1067]	58	LOGICAL, SAVE :: first=.TRUE.
[1085]	59	!$OMP THREADPRIVATE(first)
[1067]	60	REAL, DIMENSION(klon,klev) :: local_H
	61	REAL, DIMENSION(klon) :: psref
	62	REAL :: delz, pkh
[782]	63	INTEGER :: k, i, ierr
	64
	65	! Include
	66	!****************************************************************************************
[1067]	67	INCLUDE "compbl.h"
[782]	68
	69
	70	!****************************************************************************************
	71	! 1)
[1067]	72	! Allocation at first time step only
[782]	73	!
	74	!****************************************************************************************
	75
[1067]	76	IF (first) THEN
	77	first=.FALSE.
	78	ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
	79	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
[782]	80
[1067]	81	ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
	82	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
	83
	84	ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
	85	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
	86
	87	ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
	88	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
	89
	90	ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
	91	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
	92
	93	ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
	94	IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
	95
	96	ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
	97	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
	98
	99	ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
	100	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
	101	END IF
[782]	102
	103	!****************************************************************************************
	104	! 2)
	105	! Definition of the coeficient K
	106	!
	107	!****************************************************************************************
	108	Kcoefhq(:,:) = 0.0
	109	DO k = 2, klev
	110	DO i = 1, knon
	111	Kcoefhq(i,k) = &
[1067]	112	coefhq(i,k)RGRG*dtime /(pplay(i,k-1)-pplay(i,k)) &
[782]	113	(paprs(i,k)2/(temp(i,k)+temp(i,k-1))/RD)**2
	114	ENDDO
	115	ENDDO
	116
	117	!****************************************************************************************
	118	! 3)
	119	! Calculation of gama for "Q" and "H"
	120	!
	121	!****************************************************************************************
	122	! surface pressure is used as reference
	123	psref(:) = paprs(:,1)
	124
[1067]	125	! definition of gama
[782]	126	IF (iflag_pbl == 1) THEN
	127	gamaq(:,:) = 0.0
	128	gamah(:,:) = -1.0e-03
	129	gamah(:,2) = -2.5e-03
[1067]	130
	131	! conversion de gama
	132	DO k = 2, klev
	133	DO i = 1, knon
	134	delz = RD * (temp(i,k-1)+temp(i,k)) / &
	135	2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
	136	pkh = (psref(i)/paprs(i,k))**RKAPPA
	137
	138	! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
	139	gamaq(i,k) = gamaq(i,k) * delz
	140	! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
	141	gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
	142	ENDDO
	143	ENDDO
	144
[782]	145	ELSE
	146	gamaq(:,:) = 0.0
	147	gamah(:,:) = 0.0
	148	ENDIF
	149
	150
	151	!****************************************************************************************
	152	! 4)
	153	! Calculte the coefficients C and D for specific humidity, q
	154	!
	155	!****************************************************************************************
	156
[1067]	157	CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
	158	Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
[782]	159
	160	!****************************************************************************************
	161	! 5)
	162	! Calculte the coefficients C and D for potentiel entalpie, H
	163	!
	164	!****************************************************************************************
	165	local_H(:,:) = 0.0
	166
	167	DO k=1,klev
	168	DO i = 1, knon
[1067]	169	! convertie la temperature en entalpie potentielle
[782]	170	local_H(i,k) = RCPD * temp(i,k) * &
	171	(psref(i)/pplay(i,k))**RKAPPA
	172	ENDDO
	173	ENDDO
	174
[1067]	175	CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), local_H(:,:), &
	176	Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
[782]	177
	178	!****************************************************************************************
	179	! 6)
	180	! Return the first layer in output variables
	181	!
	182	!****************************************************************************************
[1067]	183	Acoef_H_out = Acoef_H
	184	Bcoef_H_out = Bcoef_H
	185	Acoef_Q_out = Acoef_Q
	186	Bcoef_Q_out = Bcoef_Q
[782]	187
	188	END SUBROUTINE climb_hq_down
	189	!
	190	!****************************************************************************************
	191	!
[1067]	192	SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
[782]	193	!
	194	! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
	195	! where X is H or Q, and k the vertical level k=1,klev
	196	!
[793]	197	INCLUDE "YOMCST.h"
[782]	198	! Input arguments
	199	!****************************************************************************************
	200	INTEGER, INTENT(IN) :: knon
[1067]	201	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp
[782]	202	REAL, DIMENSION(klon,klev), INTENT(IN) :: X
	203	REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
	204
	205	! Output arguments
	206	!****************************************************************************************
[1067]	207	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
	208	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef
[782]	209
	210	! Local variables
	211	!****************************************************************************************
	212	INTEGER :: k, i
	213	REAL :: buf
	214
	215	!****************************************************************************************
	216	! Niveau au sommet, k=klev
	217	!
	218	!****************************************************************************************
	219	Ccoef(:,:) = 0.0
	220	Dcoef(:,:) = 0.0
	221
	222	DO i = 1, knon
[1067]	223	buf = delp(i,klev) + Kcoef(i,klev)
[782]	224
[1067]	225	Ccoef(i,klev) = (X(i,klev)delp(i,klev) - Kcoef(i,klev)gama(i,klev))/buf
[782]	226	Dcoef(i,klev) = Kcoef(i,klev)/buf
	227	END DO
	228
	229
	230	!****************************************************************************************
	231	! Niveau (klev-1) <= k <= 2
	232	!
	233	!****************************************************************************************
	234
	235	DO k=(klev-1),2,-1
	236	DO i = 1, knon
[1067]	237	buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
	238	Ccoef(i,k) = (X(i,k)delp(i,k) + Kcoef(i,k+1)Ccoef(i,k+1) + &
[782]	239	Kcoef(i,k+1)gama(i,k+1) - Kcoef(i,k)gama(i,k))/buf
	240	Dcoef(i,k) = Kcoef(i,k)/buf
	241	END DO
	242	END DO
	243
	244	!****************************************************************************************
	245	! Niveau k=1
	246	!
	247	!****************************************************************************************
	248
	249	DO i = 1, knon
[1067]	250	buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
	251	Acoef(i) = (X(i,1)delp(i,1) + Kcoef(i,2)(gama(i,2)+Ccoef(i,2)))/buf
	252	Bcoef(i) = -1. * RG / buf
[782]	253	END DO
	254
	255	END SUBROUTINE calc_coef
	256	!
	257	!****************************************************************************************
	258	!
	259	SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
[1067]	260	flx_q1, flx_h1, paprs, pplay, &
[782]	261	flux_q, flux_h, d_q, d_t)
	262	!
	263	! This routine calculates the flux and tendency of the specific humidity q and
	264	! the potential engergi H.
	265	! The quantities q and H are calculated according to
	266	! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
	267	! C and D are known from before and k is index of the vertical layer.
	268	!
[793]	269	INCLUDE "YOMCST.h"
[782]	270	! Input arguments
	271	!****************************************************************************************
	272	INTEGER, INTENT(IN) :: knon
	273	REAL, INTENT(IN) :: dtime
	274	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_old, q_old
[1067]	275	REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1
[782]	276	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
	277	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
	278
	279	! Output arguments
	280	!****************************************************************************************
	281	REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t
	282
	283	! Local variables
	284	!****************************************************************************************
[1067]	285	LOGICAL, SAVE :: last=.FALSE.
[782]	286	REAL, DIMENSION(klon,klev) :: h_new, q_new
	287	REAL, DIMENSION(klon) :: psref
[996]	288	INTEGER :: k, i, ierr
[782]	289
	290	!****************************************************************************************
	291	! 1)
	292	! Definition of some variables
	293	!
	294	!****************************************************************************************
	295	flux_q(:,:) = 0.0
	296	flux_h(:,:) = 0.0
	297	d_q(:,:) = 0.0
	298	d_t(:,:) = 0.0
	299
	300	psref(1:knon) = paprs(1:knon,1)
	301
	302	!****************************************************************************************
	303	! 2)
	304	! Calculation of Q and H
	305	!
	306	!****************************************************************************************
	307
	308	!- First layer
[1067]	309	q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)flx_q1(1:knon)dtime
	310	h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)flx_h1(1:knon)dtime
[782]	311
[1067]	312	!- All the other layers
[782]	313	DO k = 2, klev
	314	DO i = 1, knon
	315	q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
	316	h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
	317	END DO
	318	END DO
	319	!****************************************************************************************
	320	! 3)
	321	! Calculation of the flux for Q and H
	322	!
	323	!****************************************************************************************
	324
	325	!- The flux at first layer, k=1
	326	flux_q(1:knon,1)=flx_q1(1:knon)
[1067]	327	flux_h(1:knon,1)=flx_h1(1:knon)
[782]	328
	329	!- The flux at all layers above surface
	330	DO k = 2, klev
	331	DO i = 1, knon
	332	flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
	333	(q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
	334
	335	flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
[1067]	336	(h_new(i,k)-h_new(i,k-1)+gamah(i,k))
[782]	337	END DO
	338	END DO
	339
	340	!****************************************************************************************
	341	! 4)
	342	! Calculation of tendency for Q and H
	343	!
	344	!****************************************************************************************
	345
	346	DO k = 1, klev
	347	DO i = 1, knon
[1067]	348	d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
[782]	349	d_q(i,k) = q_new(i,k) - q_old(i,k)
	350	END DO
	351	END DO
	352
	353	!****************************************************************************************
	354	! Some deallocations
	355	!
	356	!****************************************************************************************
[1067]	357	IF (last) THEN
	358	DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
	359	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
	360	DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
	361	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
	362	DEALLOCATE(gamaq, gamah,stat=ierr)
	363	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
	364	DEALLOCATE(Kcoefhq,stat=ierr)
	365	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
	366	END IF
[782]	367	END SUBROUTINE climb_hq_up
	368	!
	369	!****************************************************************************************
	370	!
	371	END MODULE climb_hq_mod
	372
[1067]	373
[782]	374
	375
	376
	377

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

Download in other formats:

Original Format