Context Navigation

climb_hq_mod.F90 @ 5278

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

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

Download in other formats:

Original Format