Context Navigation

climb_hq_mod.F90 @ 2791

Last change on this file since 2791 was 2187, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2158:2186 into testing branch.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 16.3 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 :: 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)
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, &
32	!!! nrlmd le 02/05/2011
33	Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
34	Kcoef_hq_out, gama_q_out, gama_h_out, &
35	!!!
36	Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)
37
38	! This routine calculates recursivly the coefficients C and D
39	! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
40	! the index of the vertical layer.
41	!
42	! Input arguments
43	!****************************************************************************************
44	INTEGER, INTENT(IN) :: knon
45	REAL, DIMENSION(klon,klev), INTENT(IN) :: coefhq
46	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
47	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
48	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp, delp ! temperature
49	REAL, DIMENSION(klon,klev), INTENT(IN) :: q
50	REAL, INTENT(IN) :: dtime
51
52	! Output arguments
53	!****************************************************************************************
54	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out
55	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out
56	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out
57	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out
58
59	!!! nrlmd le 02/05/2011
60	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_H_out
61	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_Q_out
62	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_H_out
63	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_Q_out
64	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_hq_out
65	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_q_out
66	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_h_out
67	!!!
68
69	! Local variables
70	!****************************************************************************************
71	LOGICAL, SAVE :: first=.TRUE.
72	!$OMP THREADPRIVATE(first)
73	REAL, DIMENSION(klon,klev) :: local_H
74	REAL, DIMENSION(klon) :: psref
75	REAL :: delz, pkh
76	INTEGER :: k, i, ierr
77
78	! Include
79	!****************************************************************************************
80	INCLUDE "YOMCST.h"
81	INCLUDE "compbl.h"
82
83
84	!****************************************************************************************
85	! 1)
86	! Allocation at first time step only
87	!
88	!****************************************************************************************
89
90	IF (first) THEN
91	first=.FALSE.
92	ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
93	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
94
95	ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
96	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
97
98	ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
99	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
100
101	ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
102	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
103
104	ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
105	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
106
107	ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
108	IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
109
110	ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
111	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
112
113	ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
114	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
115	END IF
116
117	!****************************************************************************************
118	! 2)
119	! Definition of the coeficient K
120	!
121	!****************************************************************************************
122	Kcoefhq(:,:) = 0.0
123	DO k = 2, klev
124	DO i = 1, knon
125	Kcoefhq(i,k) = &
126	coefhq(i,k)RGRG*dtime /(pplay(i,k-1)-pplay(i,k)) &
127	(paprs(i,k)2/(temp(i,k)+temp(i,k-1))/RD)**2
128	ENDDO
129	ENDDO
130
131	!****************************************************************************************
132	! 3)
133	! Calculation of gama for "Q" and "H"
134	!
135	!****************************************************************************************
136	! surface pressure is used as reference
137	psref(:) = paprs(:,1)
138
139	! definition of gama
140	IF (iflag_pbl == 1) THEN
141	gamaq(:,:) = 0.0
142	gamah(:,:) = -1.0e-03
143	gamah(:,2) = -2.5e-03
144
145	! conversion de gama
146	DO k = 2, klev
147	DO i = 1, knon
148	delz = RD * (temp(i,k-1)+temp(i,k)) / &
149	2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
150	pkh = (psref(i)/paprs(i,k))**RKAPPA
151
152	! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
153	gamaq(i,k) = gamaq(i,k) * delz
154	! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
155	gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
156	ENDDO
157	ENDDO
158
159	ELSE
160	gamaq(:,:) = 0.0
161	gamah(:,:) = 0.0
162	ENDIF
163
164
165	!****************************************************************************************
166	! 4)
167	! Calculte the coefficients C and D for specific humidity, q
168	!
169	!****************************************************************************************
170
171	CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
172	Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
173
174	!****************************************************************************************
175	! 5)
176	! Calculte the coefficients C and D for potentiel entalpie, H
177	!
178	!****************************************************************************************
179	local_H(:,:) = 0.0
180
181	DO k=1,klev
182	DO i = 1, knon
183	! convertie la temperature en entalpie potentielle
184	local_H(i,k) = RCPD * temp(i,k) * &
185	(psref(i)/pplay(i,k))**RKAPPA
186	ENDDO
187	ENDDO
188
189	CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), local_H(:,:), &
190	Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
191
192	!****************************************************************************************
193	! 6)
194	! Return the first layer in output variables
195	!
196	!****************************************************************************************
197	Acoef_H_out = Acoef_H
198	Bcoef_H_out = Bcoef_H
199	Acoef_Q_out = Acoef_Q
200	Bcoef_Q_out = Bcoef_Q
201
202	!****************************************************************************************
203	! 7)
204	! If Pbl is split, return also the other layers in output variables
205	!
206	!****************************************************************************************
207	!!! jyg le 07/02/2012
208	IF (mod(iflag_pbl_split,2) .eq.1) THEN
209	!!! nrlmd le 02/05/2011
210	DO k= 1, klev
211	DO i= 1, klon
212	Ccoef_H_out(i,k) = Ccoef_H(i,k)
213	Dcoef_H_out(i,k) = Dcoef_H(i,k)
214	Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
215	Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
216	Kcoef_hq_out(i,k) = Kcoefhq(i,k)
217	IF (k.eq.1) THEN
218	gama_h_out(i,k) = 0.
219	gama_q_out(i,k) = 0.
220	ELSE
221	gama_h_out(i,k) = gamah(i,k)
222	gama_q_out(i,k) = gamaq(i,k)
223	ENDIF
224	ENDDO
225	ENDDO
226	!!!
227	ENDIF ! (mod(iflag_pbl_split,2) .eq.1)
228	!!!
229
230	END SUBROUTINE climb_hq_down
231	!
232	!****************************************************************************************
233	!
234	SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
235	!
236	! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
237	! where X is H or Q, and k the vertical level k=1,klev
238	!
239	INCLUDE "YOMCST.h"
240	! Input arguments
241	!****************************************************************************************
242	INTEGER, INTENT(IN) :: knon
243	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp
244	REAL, DIMENSION(klon,klev), INTENT(IN) :: X
245	REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
246
247	! Output arguments
248	!****************************************************************************************
249	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
250	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef
251
252	! Local variables
253	!****************************************************************************************
254	INTEGER :: k, i
255	REAL :: buf
256
257	!****************************************************************************************
258	! Niveau au sommet, k=klev
259	!
260	!****************************************************************************************
261	Ccoef(:,:) = 0.0
262	Dcoef(:,:) = 0.0
263
264	DO i = 1, knon
265	buf = delp(i,klev) + Kcoef(i,klev)
266
267	Ccoef(i,klev) = (X(i,klev)delp(i,klev) - Kcoef(i,klev)gama(i,klev))/buf
268	Dcoef(i,klev) = Kcoef(i,klev)/buf
269	END DO
270
271
272	!****************************************************************************************
273	! Niveau (klev-1) <= k <= 2
274	!
275	!****************************************************************************************
276
277	DO k=(klev-1),2,-1
278	DO i = 1, knon
279	buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
280	Ccoef(i,k) = (X(i,k)delp(i,k) + Kcoef(i,k+1)Ccoef(i,k+1) + &
281	Kcoef(i,k+1)gama(i,k+1) - Kcoef(i,k)gama(i,k))/buf
282	Dcoef(i,k) = Kcoef(i,k)/buf
283	END DO
284	END DO
285
286	!****************************************************************************************
287	! Niveau k=1
288	!
289	!****************************************************************************************
290
291	DO i = 1, knon
292	buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
293	Acoef(i) = (X(i,1)delp(i,1) + Kcoef(i,2)(gama(i,2)+Ccoef(i,2)))/buf
294	Bcoef(i) = -1. * RG / buf
295	END DO
296
297	END SUBROUTINE calc_coef
298	!
299	!****************************************************************************************
300	!
301	SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
302	flx_q1, flx_h1, paprs, pplay, &
303	!!! nrlmd le 02/05/2011
304	Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
305	Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
306	Kcoef_hq_in, gama_q_in, gama_h_in, &
307	!!!
308	flux_q, flux_h, d_q, d_t)
309	!
310	! This routine calculates the flux and tendency of the specific humidity q and
311	! the potential engergi H.
312	! The quantities q and H are calculated according to
313	! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
314	! C and D are known from before and k is index of the vertical layer.
315	!
316
317	! Input arguments
318	!****************************************************************************************
319	INTEGER, INTENT(IN) :: knon
320	REAL, INTENT(IN) :: dtime
321	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_old, q_old
322	REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1
323	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
324	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
325
326	!!! nrlmd le 02/05/2011
327	REAL, DIMENSION(klon), INTENT(IN) :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
328	REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
329	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_hq_in, gama_q_in, gama_h_in
330	!!!
331
332	! Output arguments
333	!****************************************************************************************
334	REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t
335
336	! Local variables
337	!****************************************************************************************
338	LOGICAL, SAVE :: last=.FALSE.
339	REAL, DIMENSION(klon,klev) :: h_new, q_new
340	REAL, DIMENSION(klon) :: psref
341	INTEGER :: k, i, ierr
342
343	! Include
344	!****************************************************************************************
345	INCLUDE "YOMCST.h"
346	INCLUDE "compbl.h"
347
348	!****************************************************************************************
349	! 1)
350	! Definition of some variables
351	!
352	!****************************************************************************************
353	flux_q(:,:) = 0.0
354	flux_h(:,:) = 0.0
355	d_q(:,:) = 0.0
356	d_t(:,:) = 0.0
357
358	psref(1:knon) = paprs(1:knon,1)
359
360	!!! jyg le 07/02/2012
361	IF (mod(iflag_pbl_split,2) .eq.1) THEN
362	!!! nrlmd le 02/05/2011
363	DO i = 1, knon
364	Acoef_H(i)=Acoef_H_in(i)
365	Acoef_Q(i)=Acoef_Q_in(i)
366	Bcoef_H(i)=Bcoef_H_in(i)
367	Bcoef_Q(i)=Bcoef_Q_in(i)
368	ENDDO
369	DO k = 1, klev
370	DO i = 1, knon
371	Ccoef_H(i,k)=Ccoef_H_in(i,k)
372	Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
373	Dcoef_H(i,k)=Dcoef_H_in(i,k)
374	Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
375	Kcoefhq(i,k)=Kcoef_hq_in(i,k)
376	IF (k.gt.1) THEN
377	gamah(i,k)=gama_h_in(i,k)
378	gamaq(i,k)=gama_q_in(i,k)
379	ENDIF
380	ENDDO
381	ENDDO
382	!!!
383	ENDIF ! (mod(iflag_pbl_split,2) .eq.1)
384	!!!
385
386	!****************************************************************************************
387	! 2)
388	! Calculation of Q and H
389	!
390	!****************************************************************************************
391
392	!- First layer
393	q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)flx_q1(1:knon)dtime
394	h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)flx_h1(1:knon)dtime
395
396	!- All the other layers
397	DO k = 2, klev
398	DO i = 1, knon
399	q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
400	h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
401	END DO
402	END DO
403	!****************************************************************************************
404	! 3)
405	! Calculation of the flux for Q and H
406	!
407	!****************************************************************************************
408
409	!- The flux at first layer, k=1
410	flux_q(1:knon,1)=flx_q1(1:knon)
411	flux_h(1:knon,1)=flx_h1(1:knon)
412
413	!- The flux at all layers above surface
414	DO k = 2, klev
415	DO i = 1, knon
416	flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
417	(q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
418
419	flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
420	(h_new(i,k)-h_new(i,k-1)+gamah(i,k))
421	END DO
422	END DO
423
424	!****************************************************************************************
425	! 4)
426	! Calculation of tendency for Q and H
427	!
428	!****************************************************************************************
429
430	DO k = 1, klev
431	DO i = 1, knon
432	d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
433	d_q(i,k) = q_new(i,k) - q_old(i,k)
434	END DO
435	END DO
436
437	!****************************************************************************************
438	! Some deallocations
439	!
440	!****************************************************************************************
441	IF (last) THEN
442	DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
443	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
444	DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
445	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
446	DEALLOCATE(gamaq, gamah,stat=ierr)
447	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
448	DEALLOCATE(Kcoefhq,stat=ierr)
449	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
450	END IF
451	END SUBROUTINE climb_hq_up
452	!
453	!****************************************************************************************
454	!
455	END MODULE climb_hq_mod
456
457
458
459
460
461

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

Download in other formats:

Original Format