Context Navigation

climb_hq_mod.F90 @ 3927

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

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

Download in other formats:

Original Format