Context Navigation

climb_hq_mod.F90 @ 5444

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

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

Download in other formats:

Original Format