source: LMDZ5/trunk/libf/filtrez/filtreg_mod.F90 @ 1615

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

Minor corrections in filter: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong (in extreme cases, e.g. with very little grid points combined to a zoom).
EM
File size: 14.7 KB
Line 
1!
2! $Id $
3!
4MODULE filtreg_mod
5
6  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
7  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
8
9CONTAINS
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
77250 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
13350  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
25982     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
26784     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
27887     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
28689     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
29592     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
30394     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
31297     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
32099     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    !
548334 FORMAT(1x,24i3)
549755 FORMAT(1x,6f10.3,i3)
550
551    RETURN
552  END SUBROUTINE inifilr
553
554END MODULE filtreg_mod
Note: See TracBrowser for help on using the repository browser.