Context Navigation

climb_hq_mod.F90 @ 5116

Last change on this file since 5116 was 5116, checked in by abarral, 4 months ago
rename modules properly lmdz_* move ismin, ismax, minmax into new lmdz_libmath.f90 (lint) uppercase fortran keywords
Property svn:keywords set to `Id`
File size: 28.0 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	#ifdef ISO
8	USE infotrac_phy, ONLY: ntraciso=>ntiso ! ajout C Risi pour isos
9	#endif
10
11	IMPLICIT NONE
12	SAVE
13	PRIVATE
14	PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
15
16	REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
17	!$OMP THREADPRIVATE(gamaq,gamah)
18	REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
19	!$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
20	REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
21	!$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
22	REAL, DIMENSION(:), ALLOCATABLE :: Acoef_Q, Bcoef_Q
23	!$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
24	REAL, DIMENSION(:), ALLOCATABLE :: Acoef_H, Bcoef_H
25	!$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
26	REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
27	!$OMP THREADPRIVATE(Kcoefhq)
28	REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
29	!$OMP THREADPRIVATE(h_old)
30	REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
31	!$OMP THREADPRIVATE(d_h_col_vdf)
32	REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
33	!$OMP THREADPRIVATE(f_h_bnd)
34	#ifdef ISO
35	REAL, DIMENSION(:,:,:), ALLOCATABLE :: gamaxt
36	!$OMP THREADPRIVATE(gamaxt)
37	REAL, DIMENSION(:,:,:), ALLOCATABLE :: Ccoef_XT, Dcoef_XT
38	!$OMP THREADPRIVATE(Ccoef_XT, Dcoef_XT)
39	REAL, DIMENSION(:,:), ALLOCATABLE :: Acoef_XT, Bcoef_XT
40	!$OMP THREADPRIVATE(Acoef_XT, Bcoef_XT)
41	#endif
42
43	CONTAINS
44
45	!****************************************************************************************
46
47	SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
48	delp, temp, q, dtime, &
49	!!! nrlmd le 02/05/2011
50	Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
51	Kcoef_hq_out, gama_q_out, gama_h_out, &
52	!!!
53	Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out &
54	#ifdef ISO
55	,xt, &
56	Ccoef_XT_out,Dcoef_XT_out,gama_xt_out, &
57	Acoef_XT_out, Bcoef_XT_out &
58	#endif
59	)
60	#ifdef ISOVERIF
61	USE isotopes_mod, ONLY: iso_eau,iso_HDO
62	!USE isotopes_verif_mod, ONLY: errmax, errmaxrel
63	USE isotopes_verif_mod
64	#endif
65
66	! This routine calculates recursivly the coefficients C and D
67	! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
68	! the index of the vertical layer.
69
70	! Input arguments
71	!****************************************************************************************
72	INTEGER, INTENT(IN) :: knon
73	REAL, DIMENSION(klon,klev), INTENT(IN) :: coefhq
74	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
75	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
76	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp, delp ! temperature
77	REAL, DIMENSION(klon,klev), INTENT(IN) :: q
78	REAL, INTENT(IN) :: dtime
79	#ifdef ISO
80	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt
81	#endif
82
83
84	! Output arguments
85	!****************************************************************************************
86	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out
87	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out
88	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out
89	REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out
90	#ifdef ISO
91	REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: Acoef_XT_out
92	REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: Bcoef_XT_out
93	#endif
94
95	!!! nrlmd le 02/05/2011
96	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_H_out
97	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_Q_out
98	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_H_out
99	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_Q_out
100	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_hq_out
101	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_q_out
102	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_h_out
103	#ifdef ISO
104	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: Ccoef_XT_out
105	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: Dcoef_XT_out
106	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: gama_xt_out
107	#endif
108	!!!
109
110	! Local variables
111	!****************************************************************************************
112	LOGICAL, SAVE :: first=.TRUE.
113	!$OMP THREADPRIVATE(first)
114	! JLD now renamed h_old and declared in module
115	! REAL, DIMENSION(klon,klev) :: local_H
116	REAL, DIMENSION(klon) :: psref
117	REAL :: delz, pkh
118	INTEGER :: k, i, ierr
119
120	#ifdef ISO
121	real, DIMENSION(klon,2:klev) :: gamaxt_tmp
122	real, DIMENSION(klon,klev) :: xt_tmp
123	real, DIMENSION(klon,klev) :: Ccoef_XT_tmp,Dcoef_XT_tmp
124	real, DIMENSION(klon) :: Acoef_XT_tmp,Bcoef_XT_tmp
125	integer ixt
126	#endif
127	! Include
128	!****************************************************************************************
129	INCLUDE "YOMCST.h"
130	INCLUDE "compbl.h"
131
132	#ifdef ISO
133	#ifdef ISOVERIF
134	if (iso_eau.gt.0) THEN
135	do k = 1, klev
136	DO i = 1, knon
137	CALL iso_verif_egalite_choix( &
138	xt(iso_eau,i,k),q(i,k), &
139	'climb_hq 100',errmax,errmaxrel)
140	enddo
141	enddo
142	endif ! if (iso_eau.gt.0) THEN
143	#endif
144	#endif
145
146	!****************************************************************************************
147	! 1)
148	! Allocation at first time step only
149
150	!****************************************************************************************
151
152	IF (first) THEN
153	first=.FALSE.
154	ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
155	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
156
157	ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
158	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
159
160	ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
161	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
162
163	ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
164	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
165
166	ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
167	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
168
169	ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
170	IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
171
172	ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
173	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
174
175	ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
176	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
177
178	#ifdef ISO
179	ALLOCATE(Ccoef_XT(ntraciso,klon,klev), STAT=ierr)
180	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_XT, ierr=', ierr
181
182	ALLOCATE(Dcoef_XT(ntraciso,klon,klev), STAT=ierr)
183	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_XT, ierr=', ierr
184
185	ALLOCATE(Acoef_XT(ntraciso,klon), Bcoef_XT(ntraciso,klon), STAT=ierr)
186	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_XT and Bcoef_XT, ierr=', ierr
187
188	ALLOCATE(gamaxt(ntraciso,1:klon,2:klev), STAT=ierr)
189	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaxt, ierr=', ierr
190	#endif
191
192	ALLOCATE(h_old(klon,klev), STAT=ierr)
193	IF ( ierr /= 0 ) PRINT*,' pb in allloc h_old, ierr=', ierr
194
195	ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
196	IF ( ierr /= 0 ) PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
197
198	ALLOCATE(f_h_bnd(klon), STAT=ierr)
199	IF ( ierr /= 0 ) PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
200	END IF
201
202	!****************************************************************************************
203	! 2)
204	! Definition of the coeficient K
205
206	!****************************************************************************************
207	Kcoefhq(:,:) = 0.0
208	DO k = 2, klev
209	DO i = 1, knon
210	Kcoefhq(i,k) = &
211	coefhq(i,k)RGRG*dtime /(pplay(i,k-1)-pplay(i,k)) &
212	(paprs(i,k)2/(temp(i,k)+temp(i,k-1))/RD)**2
213	ENDDO
214	ENDDO
215
216	!****************************************************************************************
217	! 3)
218	! Calculation of gama for "Q" and "H"
219
220	!****************************************************************************************
221	! surface pressure is used as reference
222	psref(:) = paprs(:,1)
223
224	! definition of gama
225	IF (iflag_pbl == 1) THEN
226	gamaq(:,:) = 0.0
227	gamah(:,:) = -1.0e-03
228	gamah(:,2) = -2.5e-03
229	#ifdef ISO
230	do ixt=1,ntraciso
231	gamaxt(:,:,:) = 0.0
232	enddo ! do ixt=1,ntraciso
233	#endif
234
235	! conversion de gama
236	DO k = 2, klev
237	DO i = 1, knon
238	delz = RD * (temp(i,k-1)+temp(i,k)) / &
239	2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
240	pkh = (psref(i)/paprs(i,k))**RKAPPA
241
242	! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
243	gamaq(i,k) = gamaq(i,k) * delz
244	! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
245	gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
246	#ifdef ISO
247	do ixt=1,ntraciso
248	gamaxt(ixt,i,k) = gamaxt(ixt,i,k) * delz
249	enddo
250	#endif
251	ENDDO
252	ENDDO
253
254	ELSE
255	gamaq(:,:) = 0.0
256	gamah(:,:) = 0.0
257	#ifdef ISO
258	do ixt=1,ntraciso
259	gamaxt(:,:,:) = 0.0
260	enddo ! do ixt=1,ntraciso
261	#endif
262	ENDIF
263
264	#ifdef ISO
265	#ifdef ISOVERIF
266	do k = 2, klev
267	DO i = 1, knon
268	CALL iso_verif_egalite_choix( &
269	gamaxt(iso_eau,i,k),gamaq(i,k), &
270	'climb_hq 209',errmax,errmaxrel)
271	enddo
272	enddo
273	#endif
274	#endif
275
276
277	!****************************************************************************************
278	! 4)
279	! Calculte the coefficients C and D for specific humidity, q
280
281	!****************************************************************************************
282
283	CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
284	Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
285
286
287	#ifdef ISO
288	do ixt=1,ntraciso
289	! compression
290	do k = 2, klev
291	DO i = 1, knon
292	gamaxt_tmp(i,k)=gamaxt(ixt,i,k)
293	enddo
294	enddo !do k = 2, klev
295	do k = 1, klev
296	DO i = 1, knon
297	xt_tmp(i,k)=xt(ixt,i,k)
298	enddo
299	enddo !do k = 2, klev
300	!appel routine generique
301	CALL calc_coef(knon, Kcoefhq(:,:), gamaxt_tmp(:,:), delp(:,:), xt_tmp(:,:), &
302	Ccoef_XT_tmp(:,:), Dcoef_XT_tmp(:,:), Acoef_XT_tmp, Bcoef_XT_tmp)
303
304	! decompression
305	do k = 1, klev
306	DO i = 1, knon
307	Ccoef_XT(ixt,i,k)=Ccoef_XT_tmp(i,k)
308	Dcoef_XT(ixt,i,k)=Dcoef_XT_tmp(i,k)
309	enddo
310	enddo !do k = 2, klev
311	DO i = 1, knon
312	Acoef_XT(ixt,i)=Acoef_XT_tmp(i)
313	Bcoef_XT(ixt,i)=Bcoef_XT_tmp(i)
314	enddo
315	enddo ! do ixt=1,ntraciso
316	#ifdef ISOVERIF
317	if (iso_eau.gt.0) THEN
318	do k = 1, klev
319	DO i = 1, knon
320	CALL iso_verif_egalite_choix( &
321	Ccoef_XT(iso_eau,i,k),Ccoef_Q(i,k), &
322	'climb_hq 234c',errmax,errmaxrel)
323	CALL iso_verif_egalite_choix( &
324	Dcoef_XT(iso_eau,i,k),Dcoef_Q(i,k), &
325	'climb_hq 234d',errmax,errmaxrel)
326	enddo !DO i = 1, knon
327	enddo !do k = 2, klev
328	DO i = 1, knon
329	CALL iso_verif_egalite_choix( &
330	Acoef_XT(iso_eau,i),Acoef_Q(i), &
331	'climb_hq 234a',errmax,errmaxrel)
332	CALL iso_verif_egalite_choix( &
333	Bcoef_XT(iso_eau,i),Bcoef_Q(i), &
334	'climb_hq 234b',errmax,errmaxrel)
335	enddo !DO i = 1, knon
336	endif !if (iso_eau.gt.0) THEN
337	#endif
338	#endif
339	!****************************************************************************************
340	! 5)
341	! Calculte the coefficients C and D for potentiel entalpie, H
342
343	!****************************************************************************************
344	h_old(:,:) = 0.0
345
346	DO k=1,klev
347	DO i = 1, knon
348	! convertie la temperature en entalpie potentielle
349	h_old(i,k) = RCPD * temp(i,k) * &
350	(psref(i)/pplay(i,k))**RKAPPA
351	ENDDO
352	ENDDO
353
354	CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
355	Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
356
357	!****************************************************************************************
358	! 6)
359	! Return the first layer in output variables
360
361	!****************************************************************************************
362	Acoef_H_out = Acoef_H
363	Bcoef_H_out = Bcoef_H
364	Acoef_Q_out = Acoef_Q
365	Bcoef_Q_out = Bcoef_Q
366	#ifdef ISO
367	Acoef_XT_out = Acoef_XT
368	Bcoef_XT_out = Bcoef_XT
369	#endif
370
371	!****************************************************************************************
372	! 7)
373	! If Pbl is split, return also the other layers in output variables
374
375	!****************************************************************************************
376	!!! jyg le 07/02/2012
377	!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
378	IF (mod(iflag_pbl_split,10) >=1) THEN
379	!!! nrlmd le 02/05/2011
380	DO k= 1, klev
381	DO i= 1, klon
382	Ccoef_H_out(i,k) = Ccoef_H(i,k)
383	Dcoef_H_out(i,k) = Dcoef_H(i,k)
384	Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
385	Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
386	Kcoef_hq_out(i,k) = Kcoefhq(i,k)
387	#ifdef ISO
388	do ixt=1,ntraciso
389	Ccoef_XT_out(ixt,i,k) = Ccoef_XT(ixt,i,k)
390	Dcoef_XT_out(ixt,i,k) = Dcoef_XT(ixt,i,k)
391	enddo
392	#endif
393	IF (k==1) THEN
394	gama_h_out(i,k) = 0.
395	gama_q_out(i,k) = 0.
396	#ifdef ISO
397	do ixt=1,ntraciso
398	gama_xt_out(ixt,i,k) = 0.
399	enddo
400	#endif
401	ELSE
402	gama_h_out(i,k) = gamah(i,k)
403	gama_q_out(i,k) = gamaq(i,k)
404	#ifdef ISO
405	do ixt=1,ntraciso
406	gama_xt_out(ixt,i,k) = gamaxt(ixt,i,k)
407	enddo
408	#endif
409	ENDIF
410	ENDDO
411	ENDDO
412	!!!
413	ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
414	!!!
415
416	END SUBROUTINE climb_hq_down
417
418	!****************************************************************************************
419
420	SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
421
422	! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
423	! where X is H or Q, and k the vertical level k=1,klev
424
425	INCLUDE "YOMCST.h"
426	! Input arguments
427	!****************************************************************************************
428	INTEGER, INTENT(IN) :: knon
429	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp
430	REAL, DIMENSION(klon,klev), INTENT(IN) :: X
431	REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
432
433	! Output arguments
434	!****************************************************************************************
435	REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
436	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef
437
438	! Local variables
439	!****************************************************************************************
440	INTEGER :: k, i
441	REAL :: buf
442
443	!****************************************************************************************
444	! Niveau au sommet, k=klev
445
446	!****************************************************************************************
447	Ccoef(:,:) = 0.0
448	Dcoef(:,:) = 0.0
449
450	DO i = 1, knon
451	buf = delp(i,klev) + Kcoef(i,klev)
452
453	Ccoef(i,klev) = (X(i,klev)delp(i,klev) - Kcoef(i,klev)gama(i,klev))/buf
454	Dcoef(i,klev) = Kcoef(i,klev)/buf
455	END DO
456
457
458	!****************************************************************************************
459	! Niveau (klev-1) <= k <= 2
460
461	!****************************************************************************************
462
463	DO k=(klev-1),2,-1
464	DO i = 1, knon
465	buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
466	Ccoef(i,k) = (X(i,k)delp(i,k) + Kcoef(i,k+1)Ccoef(i,k+1) + &
467	Kcoef(i,k+1)gama(i,k+1) - Kcoef(i,k)gama(i,k))/buf
468	Dcoef(i,k) = Kcoef(i,k)/buf
469	END DO
470	END DO
471
472	!****************************************************************************************
473	! Niveau k=1
474
475	!****************************************************************************************
476
477	DO i = 1, knon
478	buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
479	Acoef(i) = (X(i,1)delp(i,1) + Kcoef(i,2)(gama(i,2)+Ccoef(i,2)))/buf
480	Bcoef(i) = -1. * RG / buf
481	END DO
482
483	END SUBROUTINE calc_coef
484
485	!****************************************************************************************
486
487	SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
488	flx_q1, flx_h1, paprs, pplay, &
489	!!! nrlmd le 02/05/2011
490	Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
491	Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
492	Kcoef_hq_in, gama_q_in, gama_h_in, &
493	!!!
494	flux_q, flux_h, d_q, d_t &
495	#ifdef ISO
496	,xt_old, flx_xt1, &
497	Acoef_XT_in,Bcoef_XT_in,Ccoef_XT_in,Dcoef_XT_in,gama_xt_in, &
498	flux_xt, d_xt &
499	#endif
500	)
501
502	#ifdef ISOVERIF
503	USE infotrac_phy, ONLY: nzone
504	USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18, ridicule
505	USE isotopes_verif_mod
506	#endif
507
508	! This routine calculates the flux and tendency of the specific humidity q and
509	! the potential engergi H.
510	! The quantities q and H are calculated according to
511	! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
512	! C and D are known from before and k is index of the vertical layer.
513
514	! Input arguments
515	!****************************************************************************************
516	INTEGER, INTENT(IN) :: knon
517	REAL, INTENT(IN) :: dtime
518	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_old, q_old
519	REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1
520	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
521	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
522
523	!!! nrlmd le 02/05/2011
524	REAL, DIMENSION(klon), INTENT(IN) :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
525	REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
526	REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_hq_in, gama_q_in, gama_h_in
527	#ifdef ISO
528	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: Acoef_XT_in, Bcoef_XT_in
529	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: Ccoef_XT_in, Dcoef_XT_in
530	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: gama_xt_in
531	#endif
532	!!!
533
534	! Output arguments
535	!****************************************************************************************
536	REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t
537
538	! Local variables
539	!****************************************************************************************
540	LOGICAL, SAVE :: last=.FALSE.
541	REAL, DIMENSION(klon,klev) :: h_new, q_new
542	REAL, DIMENSION(klon) :: psref
543	INTEGER :: k, i, ierr
544
545	#ifdef ISO
546	REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN) :: xt_old
547	REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT) :: flux_xt, d_xt
548	REAL, DIMENSION(ntraciso,klon,klev) :: xt_new
549	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: flx_xt1
550	integer ixt
551	!#ifdef ISOVERIF
552	! integer iso_verif_noNaN_nostop
553	!#endif
554	#endif
555
556	! Include
557	!****************************************************************************************
558	INCLUDE "YOMCST.h"
559	INCLUDE "compbl.h"
560
561	!****************************************************************************************
562	! 1)
563	! Definition of some variables
564	REAL, DIMENSION(klon,klev) :: d_h, zairm
565
566	!****************************************************************************************
567
568
569	#ifdef ISO
570	#ifdef ISOVERIF
571	DO k = 1, klev
572	DO i = 1, knon
573	if (iso_eau.gt.0) THEN
574	CALL iso_verif_egalite(xt_old(iso_eau,i,k), &
575	q_old(i,k),'climb_hq_mod 421')
576	endif
577	#ifdef ISOTRAC
578	IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 422')
579	#endif
580	enddo
581	enddo
582	#endif
583	#endif
584
585	flux_q(:,:) = 0.0
586	flux_h(:,:) = 0.0
587	d_q(:,:) = 0.0
588	d_t(:,:) = 0.0
589	d_h(:,:) = 0.0
590	f_h_bnd(:)= 0.0
591	#ifdef ISO
592	q_new(:,:) = 0.0
593	flux_xt(:,:,:) = 0.0
594	d_xt(:,:,:) = 0.0
595	xt_new(:,:,:) = 0.0
596	#endif
597	psref(1:knon) = paprs(1:knon,1)
598
599	!!! jyg le 07/02/2012
600	!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
601	IF (mod(iflag_pbl_split,10) >=1) THEN
602	!!! nrlmd le 02/05/2011
603	DO i = 1, knon
604	Acoef_H(i)=Acoef_H_in(i)
605	Acoef_Q(i)=Acoef_Q_in(i)
606	Bcoef_H(i)=Bcoef_H_in(i)
607	Bcoef_Q(i)=Bcoef_Q_in(i)
608	#ifdef ISO
609	do ixt=1,ntraciso
610	Acoef_XT(ixt,i)=Acoef_XT_in(ixt,i)
611	Bcoef_XT(ixt,i)=Bcoef_XT_in(ixt,i)
612	enddo
613	#endif
614	ENDDO
615	DO k = 1, klev
616	DO i = 1, knon
617	Ccoef_H(i,k)=Ccoef_H_in(i,k)
618	Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
619	Dcoef_H(i,k)=Dcoef_H_in(i,k)
620	Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
621	Kcoefhq(i,k)=Kcoef_hq_in(i,k)
622	#ifdef ISO
623	do ixt=1,ntraciso
624	Ccoef_XT(ixt,i,k)=Ccoef_XT_in(ixt,i,k)
625	Dcoef_XT(ixt,i,k)=Dcoef_XT_in(ixt,i,k)
626	enddo
627	#endif
628	IF (k>1) THEN
629	gamah(i,k)=gama_h_in(i,k)
630	gamaq(i,k)=gama_q_in(i,k)
631	#ifdef ISO
632	do ixt=1,ntraciso
633	gamaxt(ixt,i,k)=gama_xt_in(ixt,i,k)
634	enddo
635	#endif
636	ENDIF
637	ENDDO
638	ENDDO
639	!!!
640	ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
641	!!!
642
643	!****************************************************************************************
644	! 2)
645	! Calculation of Q and H
646
647	!****************************************************************************************
648
649	!- First layer
650	q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)flx_q1(1:knon)dtime
651	h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)flx_h1(1:knon)dtime
652	f_h_bnd(1:knon) = flx_h1(1:knon)
653	#ifdef ISO
654	do ixt=1,ntraciso
655	xt_new(ixt,1:knon,1) = Acoef_XT(ixt,1:knon) + Bcoef_XT(ixt,1:knon)flx_xt1(ixt,1:knon)dtime
656	enddo ! do ixt=1,ntraciso
657	#endif
658	!- All the other layers
659	DO k = 2, klev
660	DO i = 1, knon
661	q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
662	h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
663	#ifdef ISO
664	do ixt=1,ntraciso
665	xt_new(ixt,i,k) = Ccoef_XT(ixt,i,k) + Dcoef_XT(ixt,i,k)*xt_new(ixt,i,k-1)
666	enddo ! do ixt=1,ntraciso
667	#endif
668	END DO
669	END DO
670
671	#ifdef ISO
672	#ifdef ISOVERIF
673	DO k = 1, klev
674	DO i = 1, knon
675	do ixt=1,ntraciso
676	if (iso_verif_noNaN_nostop(xt_new(ixt,i,k),'climb_hq 507').eq.1) THEN
677	WRITE(,) 'Acoef_XT(ixt,i)=',Acoef_XT(ixt,i)
678	WRITE(,) 'Bcoef_XT(ixt,i)=',Bcoef_XT(ixt,i)
679	WRITE(,) 'flx_xt1(ixt,i)=',flx_xt1(ixt,i)
680	if (k.ge.2) THEN
681	WRITE(,) 'Ccoef_XT(ixt,i,k)=',Ccoef_XT(ixt,i,k)
682	WRITE(,) 'Dcoef_XT(ixt,i,k)=',Dcoef_XT(ixt,i,k)
683	endif
684	stop
685	endif
686	enddo !do ixt=1,ntraciso
687	END DO
688	END DO
689	#endif
690	#ifdef ISOVERIF
691	if (iso_eau.gt.0) THEN
692	CALL iso_verif_egalite_vect2D( &
693	xt_new,q_new, &
694	'climb_hq_mod 504',ntraciso,klon,klev)
695	endif !if (iso_eau.gt.0) THEN
696	if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) THEN
697	do i=1,klon
698	do k=1,klev
699	if (q_new(i,k).gt.ridicule) THEN
700	if (iso_verif_o18_aberrant_nostop( &
701	xt_new(iso_HDO,i,k)/q_new(i,k), &
702	xt_new(iso_O18,i,k)/q_new(i,k), &
703	'climb_hq_mod 690').eq.1) THEN
704	WRITE(,) 'i,k,q_new(i,k)=',i,k,q_new(i,k)
705	stop
706	endif ! if (iso_verif_o18_aberrant_nostop
707	endif !if (q_seri(i,k).gt.errmax) THEN
708	enddo !k=1,klev
709	enddo !i=1,klon
710	endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) THEN
711	#endif
712	#endif
713	!****************************************************************************************
714	! 3)
715	! Calculation of the flux for Q and H
716
717	!****************************************************************************************
718
719	!- The flux at first layer, k=1
720	flux_q(1:knon,1)=flx_q1(1:knon)
721	flux_h(1:knon,1)=flx_h1(1:knon)
722	#ifdef ISO
723	do ixt=1,ntraciso
724	flux_xt(ixt,1:knon,1)=flx_xt1(ixt,1:knon)
725	enddo ! do ixt=1,ntraciso
726	#endif
727
728	!- The flux at all layers above surface
729	DO k = 2, klev
730	DO i = 1, knon
731	flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
732	(q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
733
734	flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
735	(h_new(i,k)-h_new(i,k-1)+gamah(i,k))
736
737	#ifdef ISO
738	do ixt=1,ntraciso
739	flux_xt(ixt,i,k) = (Kcoefhq(i,k)/RG/dtime) * &
740	(xt_new(ixt,i,k)-xt_new(ixt,i,k-1)+gamaxt(ixt,i,k))
741	enddo ! do ixt=1,ntraciso
742	#endif
743	END DO
744	END DO
745
746	!****************************************************************************************
747	! 4)
748	! Calculation of tendency for Q and H
749
750	!****************************************************************************************
751	d_h_col_vdf(:) = 0.0
752	DO k = 1, klev
753	DO i = 1, knon
754	d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
755	d_q(i,k) = q_new(i,k) - q_old(i,k)
756	d_h(i,k) = h_new(i,k) - h_old(i,k)
757	!JLD d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD !correction a venir
758	! layer air mass
759	zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
760	d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
761	#ifdef ISO
762	do ixt=1,ntraciso
763	d_xt(ixt,i,k) = xt_new(ixt,i,k) - xt_old(ixt,i,k)
764	enddo ! do ixt=1,ntraciso
765	#ifdef ISOVERIF
766	do ixt=1,ntraciso
767	CALL iso_verif_noNaN(xt_new(ixt,i,k),'climb_hq 562')
768	CALL iso_verif_noNaN(d_xt(ixt,i,k),'climb_hq 562')
769	enddo ! do ixt=1,ntraciso
770	#endif
771	#ifdef ISOVERIF
772	if (iso_eau.gt.0) THEN
773	CALL iso_verif_egalite(d_xt(iso_eau,i,k), &
774	d_q(i,k),'climb_hq_mod 503')
775	CALL iso_verif_egalite(xt_new(iso_eau,i,k), &
776	q_new(i,k),'climb_hq_mod 503b')
777	endif
778	#ifdef ISOTRAC
779	IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 526')
780	#endif
781	#endif
782	#endif
783	END DO
784	END DO
785
786	#ifdef ISO
787	#ifdef ISOVERIF
788	! WRITE(,) 'climb_hq_mod 758: d_xt,d_q=',d_xt(iso_eau,1,1),d_q(1,1)
789	if (iso_eau.gt.0) THEN
790	CALL iso_verif_egalite_vect2D( &
791	d_xt,d_q, &
792	'climb_hq_mod 761',ntraciso,klon,klev)
793	endif
794	#endif
795	#endif
796
797	!****************************************************************************************
798	! Some deallocations
799
800	!****************************************************************************************
801	IF (last) THEN
802	DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
803	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
804	DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
805	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
806	DEALLOCATE(gamaq, gamah,stat=ierr)
807	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
808	DEALLOCATE(Kcoefhq,stat=ierr)
809	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
810	DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
811	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
812
813	#ifdef ISO
814	DEALLOCATE(Ccoef_XT, Dcoef_XT ,stat=ierr)
815	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_XT, Dcoef_XT, ierr=', ierr
816	DEALLOCATE(Acoef_XT, Bcoef_XT, stat=ierr)
817	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_XT, Bcoef_XT, ierr=', ierr
818	DEALLOCATE(gamaxt, stat=ierr)
819	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaxt, ierr=', ierr
820	#endif
821	END IF
822	END SUBROUTINE climb_hq_up
823
824	!****************************************************************************************
825
826	END MODULE climb_hq_mod
827
828
829
830
831
832

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

Download in other formats:

Original Format