Context Navigation

filtreg_mod.F90 @ 5304

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

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

Download in other formats:

Original Format