Context Navigation

climb_hq_mod.F90 @ 5209

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

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

Download in other formats:

Original Format