Context Navigation

climb_hq_mod.f90 @ 5300

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

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

Download in other formats:

Original Format