Context Navigation

climb_hq_mod.F90 @ 5322

Last change on this file since 5322 was 5296, checked in by abarral, 2 weeks ago
Turn compbl.h into a module Move calcul_REGDYN.h to obsolete Create phys_constants_mod.f90
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	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	!****************************************************************************************
557	! 1)
558	! Definition of some variables
559	REAL, DIMENSION(klon,klev) :: d_h, zairm
560	!
561	!****************************************************************************************
562
563
564	#ifdef ISO
565	#ifdef ISOVERIF
566	DO k = 1, klev
567	DO i = 1, knon
568	if (iso_eau.gt.0) then
569	call iso_verif_egalite(xt_old(iso_eau,i,k), &
570	& q_old(i,k),'climb_hq_mod 421')
571	endif
572	#ifdef ISOTRAC
573	IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 422')
574	#endif
575	enddo
576	enddo
577	#endif
578	#endif
579
580	flux_q(:,:) = 0.0
581	flux_h(:,:) = 0.0
582	d_q(:,:) = 0.0
583	d_t(:,:) = 0.0
584	d_h(:,:) = 0.0
585	f_h_bnd(:)= 0.0
586	#ifdef ISO
587	q_new(:,:) = 0.0
588	flux_xt(:,:,:) = 0.0
589	d_xt(:,:,:) = 0.0
590	xt_new(:,:,:) = 0.0
591	#endif
592	psref(1:knon) = paprs(1:knon,1)
593
594	!!! jyg le 07/02/2012
595	!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
596	IF (mod(iflag_pbl_split,10) .ge.1) THEN
597	!!! nrlmd le 02/05/2011
598	DO i = 1, knon
599	Acoef_H(i)=Acoef_H_in(i)
600	Acoef_Q(i)=Acoef_Q_in(i)
601	Bcoef_H(i)=Bcoef_H_in(i)
602	Bcoef_Q(i)=Bcoef_Q_in(i)
603	#ifdef ISO
604	do ixt=1,ntraciso
605	Acoef_XT(ixt,i)=Acoef_XT_in(ixt,i)
606	Bcoef_XT(ixt,i)=Bcoef_XT_in(ixt,i)
607	enddo
608	#endif
609	ENDDO
610	DO k = 1, klev
611	DO i = 1, knon
612	Ccoef_H(i,k)=Ccoef_H_in(i,k)
613	Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
614	Dcoef_H(i,k)=Dcoef_H_in(i,k)
615	Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
616	Kcoefhq(i,k)=Kcoef_hq_in(i,k)
617	#ifdef ISO
618	do ixt=1,ntraciso
619	Ccoef_XT(ixt,i,k)=Ccoef_XT_in(ixt,i,k)
620	Dcoef_XT(ixt,i,k)=Dcoef_XT_in(ixt,i,k)
621	enddo
622	#endif
623	IF (k.gt.1) THEN
624	gamah(i,k)=gama_h_in(i,k)
625	gamaq(i,k)=gama_q_in(i,k)
626	#ifdef ISO
627	do ixt=1,ntraciso
628	gamaxt(ixt,i,k)=gama_xt_in(ixt,i,k)
629	enddo
630	#endif
631	ENDIF
632	ENDDO
633	ENDDO
634	!!!
635	ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
636	!!!
637
638	!****************************************************************************************
639	! 2)
640	! Calculation of Q and H
641	!
642	!****************************************************************************************
643
644	!- First layer
645	q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)flx_q1(1:knon)dtime
646	h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)flx_h1(1:knon)dtime
647	f_h_bnd(1:knon) = flx_h1(1:knon)
648	#ifdef ISO
649	do ixt=1,ntraciso
650	xt_new(ixt,1:knon,1) = Acoef_XT(ixt,1:knon) + Bcoef_XT(ixt,1:knon)flx_xt1(ixt,1:knon)dtime
651	enddo ! do ixt=1,ntraciso
652	#endif
653	!- All the other layers
654	DO k = 2, klev
655	DO i = 1, knon
656	q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
657	h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
658	#ifdef ISO
659	do ixt=1,ntraciso
660	xt_new(ixt,i,k) = Ccoef_XT(ixt,i,k) + Dcoef_XT(ixt,i,k)*xt_new(ixt,i,k-1)
661	enddo ! do ixt=1,ntraciso
662	#endif
663	END DO
664	END DO
665
666	#ifdef ISO
667	#ifdef ISOVERIF
668	DO k = 1, klev
669	DO i = 1, knon
670	do ixt=1,ntraciso
671	if (iso_verif_noNaN_nostop(xt_new(ixt,i,k),'climb_hq 507').eq.1) then
672	write(,) 'Acoef_XT(ixt,i)=',Acoef_XT(ixt,i)
673	write(,) 'Bcoef_XT(ixt,i)=',Bcoef_XT(ixt,i)
674	write(,) 'flx_xt1(ixt,i)=',flx_xt1(ixt,i)
675	if (k.ge.2) then
676	write(,) 'Ccoef_XT(ixt,i,k)=',Ccoef_XT(ixt,i,k)
677	write(,) 'Dcoef_XT(ixt,i,k)=',Dcoef_XT(ixt,i,k)
678	endif
679	stop
680	endif
681	enddo !do ixt=1,ntraciso
682	END DO
683	END DO
684	#endif
685	#ifdef ISOVERIF
686	if (iso_eau.gt.0) then
687	call iso_verif_egalite_vect2D( &
688	& xt_new,q_new, &
689	& 'climb_hq_mod 504',ntraciso,klon,klev)
690	endif !if (iso_eau.gt.0) then
691	if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
692	do i=1,klon
693	do k=1,klev
694	if (q_new(i,k).gt.ridicule) then
695	if (iso_verif_o18_aberrant_nostop( &
696	& xt_new(iso_HDO,i,k)/q_new(i,k), &
697	& xt_new(iso_O18,i,k)/q_new(i,k), &
698	& 'climb_hq_mod 690').eq.1) then
699	write(,) 'i,k,q_new(i,k)=',i,k,q_new(i,k)
700	stop
701	endif ! if (iso_verif_o18_aberrant_nostop
702	endif !if (q_seri(i,k).gt.errmax) then
703	enddo !k=1,klev
704	enddo !i=1,klon
705	endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
706	#endif
707	#endif
708	!****************************************************************************************
709	! 3)
710	! Calculation of the flux for Q and H
711	!
712	!****************************************************************************************
713
714	!- The flux at first layer, k=1
715	flux_q(1:knon,1)=flx_q1(1:knon)
716	flux_h(1:knon,1)=flx_h1(1:knon)
717	#ifdef ISO
718	do ixt=1,ntraciso
719	flux_xt(ixt,1:knon,1)=flx_xt1(ixt,1:knon)
720	enddo ! do ixt=1,ntraciso
721	#endif
722
723	!- The flux at all layers above surface
724	DO k = 2, klev
725	DO i = 1, knon
726	flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
727	(q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
728
729	flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
730	(h_new(i,k)-h_new(i,k-1)+gamah(i,k))
731
732	#ifdef ISO
733	do ixt=1,ntraciso
734	flux_xt(ixt,i,k) = (Kcoefhq(i,k)/RG/dtime) * &
735	(xt_new(ixt,i,k)-xt_new(ixt,i,k-1)+gamaxt(ixt,i,k))
736	enddo ! do ixt=1,ntraciso
737	#endif
738	END DO
739	END DO
740
741	!****************************************************************************************
742	! 4)
743	! Calculation of tendency for Q and H
744	!
745	!****************************************************************************************
746	d_h_col_vdf(:) = 0.0
747	DO k = 1, klev
748	DO i = 1, knon
749	d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
750	d_q(i,k) = q_new(i,k) - q_old(i,k)
751	d_h(i,k) = h_new(i,k) - h_old(i,k)
752	!JLD d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD !correction a venir
753	! layer air mass
754	zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
755	d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
756	#ifdef ISO
757	do ixt=1,ntraciso
758	d_xt(ixt,i,k) = xt_new(ixt,i,k) - xt_old(ixt,i,k)
759	enddo ! do ixt=1,ntraciso
760	#ifdef ISOVERIF
761	do ixt=1,ntraciso
762	call iso_verif_noNaN(xt_new(ixt,i,k),'climb_hq 562')
763	call iso_verif_noNaN(d_xt(ixt,i,k),'climb_hq 562')
764	enddo ! do ixt=1,ntraciso
765	#endif
766	#ifdef ISOVERIF
767	if (iso_eau.gt.0) then
768	call iso_verif_egalite(d_xt(iso_eau,i,k), &
769	& d_q(i,k),'climb_hq_mod 503')
770	call iso_verif_egalite(xt_new(iso_eau,i,k), &
771	& q_new(i,k),'climb_hq_mod 503b')
772	endif
773	#ifdef ISOTRAC
774	IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 526')
775	#endif
776	#endif
777	#endif
778	END DO
779	END DO
780
781	#ifdef ISO
782	#ifdef ISOVERIF
783	! write(,) 'climb_hq_mod 758: d_xt,d_q=',d_xt(iso_eau,1,1),d_q(1,1)
784	if (iso_eau.gt.0) then
785	call iso_verif_egalite_vect2D( &
786	d_xt,d_q, &
787	'climb_hq_mod 761',ntraciso,klon,klev)
788	endif
789	#endif
790	#endif
791
792	!****************************************************************************************
793	! Some deallocations
794	!
795	!****************************************************************************************
796	IF (last) THEN
797	DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
798	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
799	DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
800	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
801	DEALLOCATE(gamaq, gamah,stat=ierr)
802	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
803	DEALLOCATE(Kcoefhq,stat=ierr)
804	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
805	DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
806	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
807
808	#ifdef ISO
809	DEALLOCATE(Ccoef_XT, Dcoef_XT ,stat=ierr)
810	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_XT, Dcoef_XT, ierr=', ierr
811	DEALLOCATE(Acoef_XT, Bcoef_XT, stat=ierr)
812	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_XT, Bcoef_XT, ierr=', ierr
813	DEALLOCATE(gamaxt, stat=ierr)
814	IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaxt, ierr=', ierr
815	#endif
816	END IF
817	END SUBROUTINE climb_hq_up
818	!
819	!****************************************************************************************
820	!
821	END MODULE climb_hq_mod
822
823
824
825
826
827

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

Download in other formats:

Original Format