source: LMDZ6/trunk/libf/filtrez/filtreg_mod.F90 @ 5271

Last change on this file since 5271 was 5271, checked in by abarral, 28 hours ago

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