Context Navigation

climb_hq_mod.F90 @ 845

Last change on this file since 845 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

Line
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
30	INCLUDE "YOMCST.h"
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	!
185	INCLUDE "YOMCST.h"
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	!
257	INCLUDE "YOMCST.h"
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