Context Navigation

climb_hq_mod.F90 @ 5276

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

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

Download in other formats:

Original Format