Context Navigation

pentes_ini.F @ 2306

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

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

Download in other formats:

Original Format