Context Navigation

climb_hq_mod.F90

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

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

Download in other formats:

Original Format