Context Navigation

pentes_ini.F @ 2757

Last change on this file since 2757 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 13.4 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
5
6	USE comconst_mod, ONLY: dtvr,pi
7
8	IMPLICIT NONE
9
10	c=======================================================================
11	c Adaptation LMDZ: A.Armengaud (LGGE)
12	c ----------------
13	c
14	c ********************************************************************
15	c Transport des traceurs par la methode des pentes
16	c ********************************************************************
17	c Reference possible : Russel. G.L., Lerner J.A.:
18	c A new Finite-Differencing Scheme for Traceur Transport
19	c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
20	c ********************************************************************
21	c q,w,masse,pbaru et pbarv
22	c sont des arguments d'entree pour le s-pg ....
23	c
24	c=======================================================================
25
26
27	#include "dimensions.h"
28	#include "paramet.h"
29	#include "comgeom2.h"
30
31	c Arguments:
32	c ----------
33	integer mode
34	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
35	REAL q( iip1,jjp1,llm,0:3)
36	REAL w( ip1jmp1,llm )
37	REAL masse( iip1,jjp1,llm)
38	c Local:
39	c ------
40	LOGICAL limit
41	REAL sm ( iip1,jjp1, llm )
42	REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm )
43	REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm )
44	real masn,mass,zz
45	INTEGER i,j,l,iq
46
47	c modif Fred 24 03 96
48
49	real sinlon(iip1),sinlondlon(iip1)
50	real coslon(iip1),coslondlon(iip1)
51	save sinlon,coslon,sinlondlon,coslondlon
52	real dyn1,dyn2,dys1,dys2
53	real qpn,qps,dqzpn,dqzps
54	real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
55	real qmin,zq,pente_max
56	c
57	REAL SSUM
58	integer ismax,ismin,lati,latf
59	EXTERNAL SSUM, convflu,ismin,ismax
60	logical first
61	save first
62	c fin modif
63
64	c EXTERNAL masskg
65	EXTERNAL advx
66	EXTERNAL advy
67	EXTERNAL advz
68
69	c modif Fred 24 03 96
70	data first/.true./
71
72	limit = .TRUE.
73	pente_max=2
74	c if (mode.eq.1.or.mode.eq.3) then
75	c if (mode.eq.1) then
76	if (mode.ge.1) then
77	lati=2
78	latf=jjm
79	else
80	lati=1
81	latf=jjp1
82	endif
83
84	qmin=0.4995
85	qmin=0.
86	if(first) then
87	print*,'SCHEMA AMONT NOUVEAU'
88	first=.false.
89	do i=2,iip1
90	coslon(i)=cos(rlonv(i))
91	sinlon(i)=sin(rlonv(i))
92	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
93	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
94	print*,coslondlon(i),sinlondlon(i)
95	enddo
96	coslon(1)=coslon(iip1)
97	coslondlon(1)=coslondlon(iip1)
98	sinlon(1)=sinlon(iip1)
99	sinlondlon(1)=sinlondlon(iip1)
100	print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)
101	print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)
102	DO l = 1,llm
103	DO j = 1,jjp1
104	DO i = 1,iip1
105	q ( i,j,l,1 )=0.
106	q ( i,j,l,2 )=0.
107	q ( i,j,l,3 )=0.
108	ENDDO
109	ENDDO
110	ENDDO
111
112	endif
113	c Fin modif Fred
114
115	c *** q contient les qqtes de traceur avant l'advection
116
117	c *** Affectation des tableaux S a partir de Q
118	c *** Rem : utilisation de SCOPY ulterieurement
119
120	DO l = 1,llm
121	DO j = 1,jjp1
122	DO i = 1,iip1
123	s0( i,j,llm+1-l ) = q ( i,j,l,0 )
124	sx( i,j,llm+1-l ) = q ( i,j,l,1 )
125	sy( i,j,llm+1-l ) = q ( i,j,l,2 )
126	sz( i,j,llm+1-l ) = q ( i,j,l,3 )
127	ENDDO
128	ENDDO
129	ENDDO
130
131	c PRINT*,'----- S0 just before conversion -------'
132	c PRINT*,'S0(16,12,1)=',s0(16,12,1)
133	c PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
134
135	c *** On calcule la masse d'air en kg
136
137	DO l = 1,llm
138	DO j = 1,jjp1
139	DO i = 1,iip1
140	sm ( i,j,llm+1-l)=masse( i,j,l )
141	ENDDO
142	ENDDO
143	ENDDO
144
145	c *** On converti les champs S en atome (resp. kg)
146	c *** Les routines d'advection traitent les champs
147	c *** a advecter si ces derniers sont en atome (resp. kg)
148	c *** A optimiser !!!
149
150	DO l = 1,llm
151	DO j = 1,jjp1
152	DO i = 1,iip1
153	s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )
154	sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )
155	sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )
156	sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )
157	ENDDO
158	ENDDO
159	ENDDO
160
161	c ss0 = 0.
162	c DO l = 1,llm
163	c DO j = 1,jjp1
164	c DO i = 1,iim
165	c ss0 = ss0 + s0 ( i,j,l )
166	c ENDDO
167	c ENDDO
168	c ENDDO
169	c PRINT*, 'valeur tot s0 avant advection=',ss0
170
171	c *** Appel des subroutines d'advection en X, en Y et en Z
172	c *** Advection avec "time-splitting"
173
174	c-----------------------------------------------------------
175	c PRINT*,'----- S0 just before ADVX -------'
176	c PRINT*,'S0(16,12,1)=',s0(16,12,1)
177
178	c-----------------------------------------------------------
179	c do l=1,llm
180	c do j=1,jjp1
181	c do i=1,iip1
182	c zq=s0(i,j,l)/sm(i,j,l)
183	c if(zq.lt.qmin)
184	c , print*,'avant advx1, s0(',i,',',j,',',l,')=',zq
185	c enddo
186	c enddo
187	c enddo
188	CCC
189	if(mode.eq.2) then
190	do l=1,llm
191	s0s=0.
192	s0n=0.
193	dyn1=0.
194	dys1=0.
195	dyn2=0.
196	dys2=0.
197	smn=0.
198	sms=0.
199	do i=1,iim
200	smn=smn+sm(i,1,l)
201	sms=sms+sm(i,jjp1,l)
202	s0n=s0n+s0(i,1,l)
203	s0s=s0s+s0(i,jjp1,l)
204	zz=sy(i,1,l)/sm(i,1,l)
205	dyn1=dyn1+sinlondlon(i)*zz
206	dyn2=dyn2+coslondlon(i)*zz
207	zz=sy(i,jjp1,l)/sm(i,jjp1,l)
208	dys1=dys1+sinlondlon(i)*zz
209	dys2=dys2+coslondlon(i)*zz
210	enddo
211	do i=1,iim
212	sy(i,1,l)=dyn1sinlon(i)+dyn2coslon(i)
213	sy(i,jjp1,l)=dys1sinlon(i)+dys2coslon(i)
214	enddo
215	do i=1,iim
216	s0(i,1,l)=s0n/smn+sy(i,1,l)
217	s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
218	enddo
219
220	s0(iip1,1,l)=s0(1,1,l)
221	s0(iip1,jjp1,l)=s0(1,jjp1,l)
222
223	do i=1,iim
224	sxn(i)=s0(i+1,1,l)-s0(i,1,l)
225	sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
226	c on rerentre les masses
227	enddo
228	do i=1,iim
229	sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
230	sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
231	s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
232	s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
233	enddo
234	sxn(iip1)=sxn(1)
235	sxs(iip1)=sxs(1)
236	do i=1,iim
237	sx(i+1,1,l)=0.25(sxn(i)+sxn(i+1))sm(i+1,1,l)
238	sx(i+1,jjp1,l)=0.25(sxs(i)+sxs(i+1))sm(i+1,jjp1,l)
239	enddo
240	s0(iip1,1,l)=s0(1,1,l)
241	s0(iip1,jjp1,l)=s0(1,jjp1,l)
242	sy(iip1,1,l)=sy(1,1,l)
243	sy(iip1,jjp1,l)=sy(1,jjp1,l)
244	sx(1,1,l)=sx(iip1,1,l)
245	sx(1,jjp1,l)=sx(iip1,jjp1,l)
246	enddo
247	endif
248
249	if (mode.eq.4) then
250	do l=1,llm
251	do i=1,iip1
252	sx(i,1,l)=0.
253	sx(i,jjp1,l)=0.
254	sy(i,1,l)=0.
255	sy(i,jjp1,l)=0.
256	enddo
257	enddo
258	endif
259	call limx(s0,sx,sm,pente_max)
260	c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
261	call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
262	c call minmaxq(zq,1.e33,-1.e33,'avant advy ')
263	if (mode.eq.4) then
264	do l=1,llm
265	do i=1,iip1
266	sx(i,1,l)=0.
267	sx(i,jjp1,l)=0.
268	sy(i,1,l)=0.
269	sy(i,jjp1,l)=0.
270	enddo
271	enddo
272	endif
273	call limy(s0,sy,sm,pente_max)
274	call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
275	c call minmaxq(zq,1.e33,-1.e33,'avant advz ')
276	do j=1,jjp1
277	do i=1,iip1
278	sz(i,j,1)=0.
279	sz(i,j,llm)=0.
280	enddo
281	enddo
282	call limz(s0,sz,sm,pente_max)
283	call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
284	if (mode.eq.4) then
285	do l=1,llm
286	do i=1,iip1
287	sx(i,1,l)=0.
288	sx(i,jjp1,l)=0.
289	sy(i,1,l)=0.
290	sy(i,jjp1,l)=0.
291	enddo
292	enddo
293	endif
294	call limy(s0,sy,sm,pente_max)
295	call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
296	do l=1,llm
297	do j=1,jjp1
298	sm(iip1,j,l)=sm(1,j,l)
299	s0(iip1,j,l)=s0(1,j,l)
300	sx(iip1,j,l)=sx(1,j,l)
301	sy(iip1,j,l)=sy(1,j,l)
302	sz(iip1,j,l)=sz(1,j,l)
303	enddo
304	enddo
305
306
307	c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
308	if (mode.eq.4) then
309	do l=1,llm
310	do i=1,iip1
311	sx(i,1,l)=0.
312	sx(i,jjp1,l)=0.
313	sy(i,1,l)=0.
314	sy(i,jjp1,l)=0.
315	enddo
316	enddo
317	endif
318	call limx(s0,sx,sm,pente_max)
319	call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
320	c call minmaxq(zq,1.e33,-1.e33,'apres advx ')
321	c do l=1,llm
322	c do j=1,jjp1
323	c do i=1,iip1
324	c zq=s0(i,j,l)/sm(i,j,l)
325	c if(zq.lt.qmin)
326	c , print*,'apres advx2, s0(',i,',',j,',',l,')=',zq
327	c enddo
328	c enddo
329	c enddo
330	c *** On repasse les S dans la variable q directement 14/10/94
331	c On revient a des rapports de melange en divisant par la masse
332
333	c En dehors des poles:
334
335	DO l = 1,llm
336	DO j = 1,jjp1
337	DO i = 1,iim
338	q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
339	q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
340	q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
341	q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
342	ENDDO
343	ENDDO
344	ENDDO
345
346	c Traitements specifiques au pole
347
348	if(mode.ge.1) then
349	DO l=1,llm
350	c filtrages aux poles
351	masn=ssum(iim,sm(1,1,l),1)
352	mass=ssum(iim,sm(1,jjp1,l),1)
353	qpn=ssum(iim,s0(1,1,l),1)/masn
354	qps=ssum(iim,s0(1,jjp1,l),1)/mass
355	dqzpn=ssum(iim,sz(1,1,l),1)/masn
356	dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
357	do i=1,iip1
358	q( i,1,llm+1-l,3)=dqzpn
359	q( i,jjp1,llm+1-l,3)=dqzps
360	q( i,1,llm+1-l,0)=qpn
361	q( i,jjp1,llm+1-l,0)=qps
362	enddo
363	if(mode.eq.3) then
364	dyn1=0.
365	dys1=0.
366	dyn2=0.
367	dys2=0.
368	do i=1,iim
369	dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
370	dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
371	dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
372	dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
373	enddo
374	do i=1,iim
375	q(i,1,llm+1-l,2)=
376	s (sinlon(i)dyn1+coslon(i)dyn2)
377	q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
378	q(i,jjp1,llm+1-l,2)=
379	s (sinlon(i)dys1+coslon(i)dys2)
380	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
381	s -q(i,jjp1,llm+1-l,2)
382	enddo
383	endif
384	if(mode.eq.1) then
385	c on filtre les valeurs au bord de la "grande maille pole"
386	dyn1=0.
387	dys1=0.
388	dyn2=0.
389	dys2=0.
390	do i=1,iim
391	zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
392	dyn1=dyn1+sinlondlon(i)*zz
393	dyn2=dyn2+coslondlon(i)*zz
394	zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
395	dys1=dys1+sinlondlon(i)*zz
396	dys2=dys2+coslondlon(i)*zz
397	enddo
398	do i=1,iim
399	q(i,1,llm+1-l,2)=
400	s (sinlon(i)dyn1+coslon(i)dyn2)/2.
401	q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
402	q(i,jjp1,llm+1-l,2)=
403	s (sinlon(i)dys1+coslon(i)dys2)/2.
404	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
405	s -q(i,jjp1,llm+1-l,2)
406	enddo
407	q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
408	q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
409
410	do i=1,iim
411	sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
412	sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
413	enddo
414	sxn(iip1)=sxn(1)
415	sxs(iip1)=sxs(1)
416	do i=1,iim
417	q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
418	q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
419	enddo
420	q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
421	q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
422
423	endif
424
425	ENDDO
426	endif
427
428	c bouclage en longitude
429	do iq=0,3
430	do l=1,llm
431	do j=1,jjp1
432	q(iip1,j,l,iq)=q(1,j,l,iq)
433	enddo
434	enddo
435	enddo
436
437	c PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....'
438
439	DO l = 1,llm
440	DO j = 1,jjp1
441	DO i = 1,iip1
442	IF (q(i,j,l,0).lt.0.) THEN
443	c PRINT*,'------------ BIP-----------'
444	c PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
445	c PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
446	c PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
447	c PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
448	c PRINT*,' PBL EN SORTIE DE PENTES'
449	q(i,j,l,0)=0.
450	c STOP
451	ENDIF
452	ENDDO
453	ENDDO
454	ENDDO
455
456	c PRINT*, '-------------------------------------------'
457
458	do l=1,llm
459	do j=1,jjp1
460	do i=1,iip1
461	if(q(i,j,l,0).lt.qmin)
462	, print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
463	enddo
464	enddo
465	enddo
466	RETURN
467	END
468
469
470
471
472
473
474
475
476
477
478
479

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

Download in other formats:

Original Format