Context Navigation

filtreg_mod.F90 @ 4672

Last change on this file since 4672 was 4519, checked in by evignon, 19 months ago
modif for polar runs: control of the maximum latitude at which the polar filter is active (new key in .def files) + dependency of ratqs upon resolution (cell area)
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
File size: 15.4 KB

Line
1	!
2	! $Id$
3	!
4	MODULE filtreg_mod
5
6	REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
7	REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
8
9	CONTAINS
10
11	SUBROUTINE inifilr
12	#ifdef CPP_PARA
13	USE mod_filtre_fft, ONLY : use_filtre_fft,Init_filtre_fft
14	USE mod_filtre_fft_loc, ONLY : Init_filtre_fft_loc=>Init_filtre_fft !
15	#endif
16	USE serre_mod, ONLY: alphax
17	USE logic_mod, ONLY: fxyhypb, ysinus
18	USE comconst_mod, ONLY: maxlatfilter
19
20	! ... H. Upadhyaya, O.Sharma ...
21	!
22	IMPLICIT NONE
23	!
24	! version 3 .....
25
26	! Correction le 28/10/97 P. Le Van .
27	! -------------------------------------------------------------------
28	include "dimensions.h"
29	include "paramet.h"
30	! -------------------------------------------------------------------
31	include "comgeom.h"
32	include "coefils.h"
33
34	REAL dlonu(iim),dlatu(jjm)
35	REAL rlamda( iim ), eignvl( iim )
36	!
37
38	REAL lamdamax,pi,cof
39	INTEGER i,j,modemax,imx,k,kf,ii
40	REAL dymin,dxmin,colat0
41	REAL eignft(iim,iim), coff
42
43	LOGICAL, SAVE :: first_call_inifilr = .TRUE.
44
45	#ifdef CRAY
46	INTEGER ISMIN
47	EXTERNAL ISMIN
48	INTEGER iymin
49	INTEGER ixmineq
50	#endif
51	!
52	! ------------------------------------------------------------
53	! This routine computes the eigenfunctions of the laplacien
54	! on the stretched grid, and the filtering coefficients
55	!
56	! We designate:
57	! eignfn eigenfunctions of the discrete laplacien
58	! eigenvl eigenvalues
59	! jfiltn indexof the last scalar line filtered in NH
60	! jfilts index of the first line filtered in SH
61	! modfrst index of the mode from WHERE modes are filtered
62	! modemax maximum number of modes ( im )
63	! coefil filtering coefficients ( lamda_max*COS(rlat)/lamda )
64	! sdd SQRT( dx )
65	!
66	! the modes are filtered from modfrst to modemax
67	!
68	!-----------------------------------------------------------
69	!
70	if ( iim == 1 ) return ! No filtre in 2D y-z
71
72	pi = 2. * ASIN( 1. )
73
74	DO i = 1,iim
75	dlonu(i) = xprimu( i )
76	ENDDO
77	!
78	CALL inifgn(eignvl)
79	!
80	PRINT *,'inifilr: EIGNVL '
81	PRINT 250,eignvl
82	250 FORMAT( 1x,5e14.6)
83	!
84	! compute eigenvalues and eigenfunctions
85	!
86	!
87	!.................................................................
88	!
89	! compute the filtering coefficients for scalar lines and
90	! meridional wind v-lines
91	!
92	! we filter all those latitude lines WHERE coefil < 1
93	! NO FILTERING AT POLES
94	!
95	! colat0 is to be used when alpha (stretching coefficient)
96	! is set equal to zero for the regular grid CASE
97	!
98	! ....... Calcul de colat0 .........
99	! ..... colat0 = minimum de ( 0.5, min dy/ min dx ) ...
100	!
101	!
102	DO j = 1,jjm
103	dlatu( j ) = rlatu( j ) - rlatu( j+1 )
104	ENDDO
105	!
106	#ifdef CRAY
107	iymin = ISMIN( jjm, dlatu, 1 )
108	ixmineq = ISMIN( iim, dlonu, 1 )
109	dymin = dlatu( iymin )
110	dxmin = dlonu( ixmineq )
111	#else
112	dxmin = dlonu(1)
113	DO i = 2, iim
114	dxmin = MIN( dxmin,dlonu(i) )
115	ENDDO
116	dymin = dlatu(1)
117	DO j = 2, jjm
118	dymin = MIN( dymin,dlatu(j) )
119	ENDDO
120	#endif
121	!
122	! For a regular grid, we want the filter to start at latitudes
123	! corresponding to lengths dx of the same size as dy (in terms
124	! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
125	! <=> latitude=60 degrees).
126	! Same idea for the zoomed grid: start filtering polewards as soon
127	! as length dx becomes of the same size as dy
128	!
129	! if maxlatfilter >0, prescribe the colat0 value from the .def files
130
131	IF (maxlatfilter .LT. 0.) THEN
132
133	colat0 = MIN( 0.5, dymin/dxmin )
134	! colat0 = 1.
135	!
136	IF( .NOT.fxyhypb.AND.ysinus ) THEN
137	colat0 = 0.6
138	! ...... a revoir pour ysinus ! .......
139	alphax = 0.
140	ENDIF
141
142	ELSE
143
144	colat0=(90.0-maxlatfilter)/180.0*pi
145
146	ENDIF
147
148
149
150	!
151	PRINT 50, colat0,alphax
152	50 FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
153	!
154	IF(alphax.EQ.1. ) THEN
155	PRINT *,' Inifilr alphax doit etre < a 1. Corriger '
156	STOP
157	ENDIF
158	!
159	lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
160
161	! ... Correction le 28/10/97 ( P.Le Van ) ..
162	!
163	DO i = 2,iim
164	rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
165	ENDDO
166	!
167
168	DO j = 1,jjm
169	DO i = 1,iim
170	coefilu( i,j ) = 0.0
171	coefilv( i,j ) = 0.0
172	coefilu2( i,j ) = 0.0
173	coefilv2( i,j ) = 0.0
174	ENDDO
175	ENDDO
176
177	!
178	! ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
179	! .........................................................
180	!
181	modemax = iim
182
183	!!!! imx = modemax - 4 * (modemax/iim)
184
185	imx = iim
186	!
187	PRINT *,'inifilr: TRUNCATION AT ',imx
188	!
189	! Ehouarn: set up some defaults
190	jfiltnu=2 ! avoid north pole
191	jfiltsu=jjm ! avoid south pole (which is at jjm+1)
192	jfiltnv=1 ! NB: no poles on the V grid
193	jfiltsv=jjm
194
195	DO j = 2, jjm/2+1
196	cof = COS( rlatu(j) )/ colat0
197	IF ( cof .LT. 1. ) THEN
198	IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
199	jfiltnu= j
200	ENDIF
201	ENDIF
202
203	cof = COS( rlatu(jjp1-j+1) )/ colat0
204	IF ( cof .LT. 1. ) THEN
205	IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
206	jfiltsu= jjp1-j+1
207	ENDIF
208	ENDIF
209	ENDDO
210	!
211	DO j = 1, jjm/2
212	cof = COS( rlatv(j) )/ colat0
213	IF ( cof .LT. 1. ) THEN
214	IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
215	jfiltnv= j
216	ENDIF
217	ENDIF
218
219	cof = COS( rlatv(jjm-j+1) )/ colat0
220	IF ( cof .LT. 1. ) THEN
221	IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
222	jfiltsv= jjm-j+1
223	ENDIF
224	ENDIF
225	ENDDO
226	!
227
228	IF( jfiltnu.GT. jjm/2 +1 ) THEN
229	PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
230	STOP
231	ENDIF
232
233	IF( jfiltsu.GT. jjm +1 ) THEN
234	PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
235	STOP
236	ENDIF
237
238	IF( jfiltnv.GT. jjm/2 ) THEN
239	PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
240	STOP
241	ENDIF
242
243	IF( jfiltsv.GT. jjm ) THEN
244	PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
245	STOP
246	ENDIF
247
248	PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
249	jfiltnv,jfiltsv,jfiltnu,jfiltsu
250
251	IF(first_call_inifilr) THEN
252	ALLOCATE(matriceun(iim,iim,jfiltnu))
253	ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
254	ALLOCATE(matricevn(iim,iim,jfiltnv))
255	ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
256	ALLOCATE( matrinvn(iim,iim,jfiltnu))
257	ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
258	first_call_inifilr = .FALSE.
259	ENDIF
260
261	!
262	! ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
263	!................................................................
264	!
265	!
266	DO j = 1,jjm
267	!default initialization: all modes are retained (i.e. no filtering)
268	modfrstu( j ) = iim
269	modfrstv( j ) = iim
270	ENDDO
271	!
272	DO j = 2,jfiltnu
273	DO k = 2,modemax
274	cof = rlamda(k) * COS( rlatu(j) )
275	IF ( cof .LT. 1. ) GOTO 82
276	ENDDO
277	GOTO 84
278	82 modfrstu( j ) = k
279	!
280	kf = modfrstu( j )
281	DO k = kf , modemax
282	cof = rlamda(k) * COS( rlatu(j) )
283	coefilu(k,j) = cof - 1.
284	coefilu2(k,j) = cof*cof - 1.
285	ENDDO
286	84 CONTINUE
287	ENDDO
288	!
289	!
290	DO j = 1,jfiltnv
291	!
292	DO k = 2,modemax
293	cof = rlamda(k) * COS( rlatv(j) )
294	IF ( cof .LT. 1. ) GOTO 87
295	ENDDO
296	GOTO 89
297	87 modfrstv( j ) = k
298	!
299	kf = modfrstv( j )
300	DO k = kf , modemax
301	cof = rlamda(k) * COS( rlatv(j) )
302	coefilv(k,j) = cof - 1.
303	coefilv2(k,j) = cof*cof - 1.
304	ENDDO
305	89 CONTINUE
306	ENDDO
307	!
308	DO j = jfiltsu,jjm
309	DO k = 2,modemax
310	cof = rlamda(k) * COS( rlatu(j) )
311	IF ( cof .LT. 1. ) GOTO 92
312	ENDDO
313	GOTO 94
314	92 modfrstu( j ) = k
315	!
316	kf = modfrstu( j )
317	DO k = kf , modemax
318	cof = rlamda(k) * COS( rlatu(j) )
319	coefilu(k,j) = cof - 1.
320	coefilu2(k,j) = cof*cof - 1.
321	ENDDO
322	94 CONTINUE
323	ENDDO
324	!
325	DO j = jfiltsv,jjm
326	DO k = 2,modemax
327	cof = rlamda(k) * COS( rlatv(j) )
328	IF ( cof .LT. 1. ) GOTO 97
329	ENDDO
330	GOTO 99
331	97 modfrstv( j ) = k
332	!
333	kf = modfrstv( j )
334	DO k = kf , modemax
335	cof = rlamda(k) * COS( rlatv(j) )
336	coefilv(k,j) = cof - 1.
337	coefilv2(k,j) = cof*cof - 1.
338	ENDDO
339	99 CONTINUE
340	ENDDO
341	!
342
343	IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
344	! Ehouarn: and what are these for??? Trying to handle a limit case
345	! where filters extend to and meet at the equator?
346	IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
347	IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
348
349	PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
350	jfiltnv,jfiltsv,jfiltnu,jfiltsu
351	ENDIF
352
353	PRINT *,' Modes premiers v '
354	PRINT 334,modfrstv
355	PRINT *,' Modes premiers u '
356	PRINT 334,modfrstu
357
358	!
359	! ...................................................................
360	!
361	! ... Calcul de la matrice filtre 'matriceu' pour les champs situes
362	! sur la grille scalaire ........
363	! ...................................................................
364	!
365	DO j = 2, jfiltnu
366
367	DO i=1,iim
368	coff = coefilu(i,j)
369	IF( i.LT.modfrstu(j) ) coff = 0.
370	DO k=1,iim
371	eignft(i,k) = eignfnv(k,i) * coff
372	ENDDO
373	ENDDO ! of DO i=1,iim
374	#ifdef CRAY
375	CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
376	#else
377	#ifdef BLAS
378	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
379	eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
380	#else
381	DO k = 1, iim
382	DO i = 1, iim
383	matriceun(i,k,j) = 0.0
384	DO ii = 1, iim
385	matriceun(i,k,j) = matriceun(i,k,j) &
386	+ eignfnv(i,ii)*eignft(ii,k)
387	ENDDO
388	ENDDO
389	ENDDO ! of DO k = 1, iim
390	#endif
391	#endif
392
393	ENDDO ! of DO j = 2, jfiltnu
394
395	DO j = jfiltsu, jjm
396
397	DO i=1,iim
398	coff = coefilu(i,j)
399	IF( i.LT.modfrstu(j) ) coff = 0.
400	DO k=1,iim
401	eignft(i,k) = eignfnv(k,i) * coff
402	ENDDO
403	ENDDO ! of DO i=1,iim
404	#ifdef CRAY
405	CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
406	#else
407	#ifdef BLAS
408	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
409	eignfnv, iim, eignft, iim, 0.0, &
410	matriceus(1,1,j-jfiltsu+1), iim)
411	#else
412	DO k = 1, iim
413	DO i = 1, iim
414	matriceus(i,k,j-jfiltsu+1) = 0.0
415	DO ii = 1, iim
416	matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
417	+ eignfnv(i,ii)*eignft(ii,k)
418	ENDDO
419	ENDDO
420	ENDDO ! of DO k = 1, iim
421	#endif
422	#endif
423
424	ENDDO ! of DO j = jfiltsu, jjm
425
426	! ...................................................................
427	!
428	! ... Calcul de la matrice filtre 'matricev' pour les champs situes
429	! sur la grille de V ou de Z ........
430	! ...................................................................
431	!
432	DO j = 1, jfiltnv
433
434	DO i = 1, iim
435	coff = coefilv(i,j)
436	IF( i.LT.modfrstv(j) ) coff = 0.
437	DO k = 1, iim
438	eignft(i,k) = eignfnu(k,i) * coff
439	ENDDO
440	ENDDO
441	#ifdef CRAY
442	CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
443	#else
444	#ifdef BLAS
445	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
446	eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
447	#else
448	DO k = 1, iim
449	DO i = 1, iim
450	matricevn(i,k,j) = 0.0
451	DO ii = 1, iim
452	matricevn(i,k,j) = matricevn(i,k,j) &
453	+ eignfnu(i,ii)*eignft(ii,k)
454	ENDDO
455	ENDDO
456	ENDDO
457	#endif
458	#endif
459
460	ENDDO ! of DO j = 1, jfiltnv
461
462	DO j = jfiltsv, jjm
463
464	DO i = 1, iim
465	coff = coefilv(i,j)
466	IF( i.LT.modfrstv(j) ) coff = 0.
467	DO k = 1, iim
468	eignft(i,k) = eignfnu(k,i) * coff
469	ENDDO
470	ENDDO
471	#ifdef CRAY
472	CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
473	#else
474	#ifdef BLAS
475	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
476	eignfnu, iim, eignft, iim, 0.0, &
477	matricevs(1,1,j-jfiltsv+1), iim)
478	#else
479	DO k = 1, iim
480	DO i = 1, iim
481	matricevs(i,k,j-jfiltsv+1) = 0.0
482	DO ii = 1, iim
483	matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
484	+ eignfnu(i,ii)*eignft(ii,k)
485	ENDDO
486	ENDDO
487	ENDDO
488	#endif
489	#endif
490
491	ENDDO ! of DO j = jfiltsv, jjm
492
493	! ...................................................................
494	!
495	! ... Calcul de la matrice filtre 'matrinv' pour les champs situes
496	! sur la grille scalaire , pour le filtre inverse ........
497	! ...................................................................
498	!
499	DO j = 2, jfiltnu
500
501	DO i = 1,iim
502	coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
503	IF( i.LT.modfrstu(j) ) coff = 0.
504	DO k=1,iim
505	eignft(i,k) = eignfnv(k,i) * coff
506	ENDDO
507	ENDDO
508	#ifdef CRAY
509	CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
510	#else
511	#ifdef BLAS
512	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
513	eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
514	#else
515	DO k = 1, iim
516	DO i = 1, iim
517	matrinvn(i,k,j) = 0.0
518	DO ii = 1, iim
519	matrinvn(i,k,j) = matrinvn(i,k,j) &
520	+ eignfnv(i,ii)*eignft(ii,k)
521	ENDDO
522	ENDDO
523	ENDDO
524	#endif
525	#endif
526
527	ENDDO ! of DO j = 2, jfiltnu
528
529	DO j = jfiltsu, jjm
530
531	DO i = 1,iim
532	coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
533	IF( i.LT.modfrstu(j) ) coff = 0.
534	DO k=1,iim
535	eignft(i,k) = eignfnv(k,i) * coff
536	ENDDO
537	ENDDO
538	#ifdef CRAY
539	CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
540	#else
541	#ifdef BLAS
542	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
543	eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
544	#else
545	DO k = 1, iim
546	DO i = 1, iim
547	matrinvs(i,k,j-jfiltsu+1) = 0.0
548	DO ii = 1, iim
549	matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
550	+ eignfnv(i,ii)*eignft(ii,k)
551	ENDDO
552	ENDDO
553	ENDDO
554	#endif
555	#endif
556
557	ENDDO ! of DO j = jfiltsu, jjm
558
559	#ifdef CPP_PARA
560	IF (use_filtre_fft) THEN
561	CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu, &
562	coefilv,modfrstv,jfiltnv,jfiltsv)
563	CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu, &
564	coefilv,modfrstv,jfiltnv,jfiltsv)
565	ENDIF
566	#endif
567	! ...................................................................
568
569	!
570	334 FORMAT(1x,24i3)
571	755 FORMAT(1x,6f10.3,i3)
572
573	RETURN
574	END SUBROUTINE inifilr
575
576	END MODULE filtreg_mod

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

Download in other formats:

Original Format