Context Navigation

climb_hq_mod.F90 @ 3774

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

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

Download in other formats:

Original Format