Context Navigation

filtreg_mod.F90 @ 5394

Last change on this file since 5394 was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM

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

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

Download in other formats:

Original Format