Context Navigation

climb_hq_mod.F90 @ 4500

Last change on this file since 4500 was 4143, checked in by dcugnet, 2 years ago
Some variables are renamed or replaced by direct equivalents: iso_indnum -> tracers(:)%iso_iName niso_possibles -> niso iqiso -> iqIsoPha ; index_trac -> itZonIso ok_iso_verif -> isoCheck ntraceurs_zone -> nzone ; ntraciso -> ntiso qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times. Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18). The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
Property svn:keywords set to `Id`
File size: 28.1 KB

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

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

Download in other formats:

Original Format