Context Navigation

climb_hq_mod.F90 @ 963

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

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

Download in other formats:

Original Format