source: LMDZ5/branches/testing/libf/filtrez/filtreg_mod.F90 @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706
File size: 15.0 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, ONLY : use_filtre_fft,Init_filtre_fft
13  USE mod_filtre_fft_loc, ONLY : Init_filtre_fft_loc=>Init_filtre_fft    !
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       CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu,  &
544                           coefilv,modfrstv,jfiltnv,jfiltsv)
545    ENDIF
546
547    !   ...................................................................
548
549    !
550334 FORMAT(1x,24i3)
551755 FORMAT(1x,6f10.3,i3)
552
553    RETURN
554  END SUBROUTINE inifilr
555
556END MODULE filtreg_mod
Note: See TracBrowser for help on using the repository browser.