Context Navigation

climb_hq_mod.F90 @ 5091

Last change on this file since 5091 was 5081, checked in by abarral, 12 months ago

Commit linked to correcting relevant warnings during gfortran compilation

[minor]
Correct remnants of fixed-form "+" in massbarxy.F90
Correct COMMON alignment in clesphys.h, flux_arp.h, cv*param.h, YOECUMF.h, alpale.h
Correct obsolete logical operators
Restrict use of iso_fortran_env in nf95_abort.f90
Remove redundant save in module declarations
Remove <continue> without labels
Fix nonstandard kind selectors

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

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

Download in other formats:

Original Format