Context Navigation

pentes_ini.F @ 5139

Last change on this file since 5139 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.5 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
5
6	USE comconst_mod, ONLY: pi, dtvr
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, 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