Context Navigation

pentes_ini.f90 @ 5106

Last change on this file since 5106 was 5106, checked in by abarral, 4 months ago
Turn coefils.h into lmdz_coefils.f90 Put filtreg.F90 inside lmdz_filtreg.F90 Turn mod_filtreg_p.F90 into lmdz_filtreg_p.F90 Delete obsolete parafilt.h* (lint) remove spaces between routine name and args
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: 12.2 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	!=======================================================================
11	! Adaptation LMDZ: A.Armengaud (LGGE)
12	! ----------------
13	!
14	! ********************************************************************
15	! Transport des traceurs par la methode des pentes
16	! ********************************************************************
17	! Reference possible : Russel. G.L., Lerner J.A.:
18	! A new Finite-Differencing Scheme for Traceur Transport
19	! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
20	! ********************************************************************
21	! q,w,masse,pbaru et pbarv
22	! sont des arguments d'entree pour le s-pg ....
23	!
24	!=======================================================================
25
26
27	include "dimensions.h"
28	include "paramet.h"
29	include "comgeom2.h"
30
31	! Arguments:
32	! ----------
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	! Local:
39	! ------
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	! 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	!
57	REAL :: SSUM
58	integer :: ismax,ismin,lati,latf
59	EXTERNAL SSUM, ismin,ismax
60	logical :: first
61	save first
62	! fin modif
63
64	! EXTERNAL masskg
65	EXTERNAL advx
66	EXTERNAL advy
67	EXTERNAL advz
68
69	! modif Fred 24 03 96
70	data first/.TRUE./
71
72	limit = .TRUE.
73	pente_max=2
74	! if (mode.eq.1.or.mode.eq.3) then
75	! if (mode.eq.1) then
76	if (mode>=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	! Fin modif Fred
114
115	! *** q contient les qqtes de traceur avant l'advection
116
117	! *** Affectation des tableaux S a partir de Q
118	! *** 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	! PRINT*,'----- S0 just before conversion -------'
132	! PRINT*,'S0(16,12,1)=',s0(16,12,1)
133	! PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
134
135	! *** 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	! *** On converti les champs S en atome (resp. kg)
146	! *** Les routines d'advection traitent les champs
147	! *** a advecter si ces derniers sont en atome (resp. kg)
148	! *** 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	! ss0 = 0.
162	! DO l = 1,llm
163	! DO j = 1,jjp1
164	! DO i = 1,iim
165	! ss0 = ss0 + s0 ( i,j,l )
166	! ENDDO
167	! ENDDO
168	! ENDDO
169	! PRINT*, 'valeur tot s0 avant advection=',ss0
170
171	! *** Appel des subroutines d'advection en X, en Y et en Z
172	! *** Advection avec "time-splitting"
173
174	!-----------------------------------------------------------
175	! PRINT*,'----- S0 just before ADVX -------'
176	! PRINT*,'S0(16,12,1)=',s0(16,12,1)
177
178	!-----------------------------------------------------------
179	! do l=1,llm
180	! do j=1,jjp1
181	! do i=1,iip1
182	! zq=s0(i,j,l)/sm(i,j,l)
183	! if(zq.lt.qmin)
184	! , PRINT*,'avant advx1, s0(',i,',',j,',',l,')=',zq
185	! enddo
186	! enddo
187	! enddo
188	!CC
189	if(mode==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	! 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==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	! CALL minmaxq(zq,1.e33,-1.e33,'avant advx ')
261	CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
262	! CALL minmaxq(zq,1.e33,-1.e33,'avant advy ')
263	if (mode==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	! 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==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	! CALL minmaxq(zq,1.e33,-1.e33,'avant advx ')
308	if (mode==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	! CALL minmaxq(zq,1.e33,-1.e33,'apres advx ')
321	! do l=1,llm
322	! do j=1,jjp1
323	! do i=1,iip1
324	! zq=s0(i,j,l)/sm(i,j,l)
325	! if(zq.lt.qmin)
326	! , PRINT*,'apres advx2, s0(',i,',',j,',',l,')=',zq
327	! enddo
328	! enddo
329	! enddo
330	! *** On repasse les S dans la variable q directement 14/10/94
331	! On revient a des rapports de melange en divisant par la masse
332
333	! 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	! Traitements specifiques au pole
347
348	if(mode>=1) then
349	DO l=1,llm
350	! 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==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	(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	(sinlon(i)dys1+coslon(i)dys2)
380	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) &
381	-q(i,jjp1,llm+1-l,2)
382	enddo
383	endif
384	if(mode==1) then
385	! 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	(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	(sinlon(i)dys1+coslon(i)dys2)/2.
404	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) &
405	-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	! 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	! 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)<0.) THEN
443	! PRINT*,'------------ BIP-----------'
444	! PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
445	! PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
446	! PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
447	! PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
448	! PRINT*,' PBL EN SORTIE DE PENTES'
449	q(i,j,l,0)=0.
450	! STOP
451	ENDIF
452	ENDDO
453	ENDDO
454	ENDDO
455
456	! PRINT*, '-------------------------------------------'
457
458	do l=1,llm
459	do j=1,jjp1
460	do i=1,iip1
461	if(q(i,j,l,0)<qmin) &
462	PRINT*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
463	enddo
464	enddo
465	enddo
466	RETURN
467	END SUBROUTINE pentes_ini
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