Context Navigation

pentes_ini.F @ 1300

Last change on this file since 1300 was 1300, checked in by emillour, 10 years ago
Common dynamics: Some updates to keep up with LMDZ5 Earth model evolution (up to LMDZ5 rev 1955). Main change is the introduction of a "dyn3d_common" directory to store files common to dyn3d and dyn3dpar. See file "DOC/chantiers/commit_importants.log" for detailed list of changes. These changes do not change results on test cases. EM
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, convflu,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