Context Navigation

climb_hq_mod.f90 @ 5289

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

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

Download in other formats:

Original Format