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

Last change on this file since 2069 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • 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.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#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    !    ... H. Upadhyaya, O.Sharma   ...
17    !
18    IMPLICIT NONE
19    !
20    !     version 3 .....
21
22    !     Correction  le 28/10/97    P. Le Van .
23    !  -------------------------------------------------------------------
24#include "dimensions.h"
25#include "paramet.h"
26    !  -------------------------------------------------------------------
27#include "comgeom.h"
28#include "coefils.h"
29#include "logic.h"
30#include "serre.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
79250 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
13550  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
26182     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
26984     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
28087     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
28889     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
29792     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
30594     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
31497     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
32299     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    !
553334 FORMAT(1x,24i3)
554755 FORMAT(1x,6f10.3,i3)
555
556    RETURN
557  END SUBROUTINE inifilr
558
559END MODULE filtreg_mod
Note: See TracBrowser for help on using the repository browser.