source: trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90 @ 3553

Last change on this file since 3553 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
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  USE logic_mod, ONLY: fxyhypb,ysinus
17  USE serre_mod, ONLY: alphax
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
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.