Context Navigation

pentes_ini.F @ 1583

Last change on this file since 1583 was 764, checked in by Laurent Fairhead, 18 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.3 KB

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

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

Download in other formats:

Original Format