LMDZ
filtreg_mod.F90
Go to the documentation of this file.
1 !
2 ! $Id$
3 !
4 MODULE filtreg_mod
5 
6  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
7  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
8 
9 CONTAINS
10 
11  SUBROUTINE inifilr
12 #ifdef CPP_PARA
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
79 250 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
135 50 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 ' , &
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
261 82 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
269 84 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
280 87 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
288 89 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
297 92 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
305 94 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
314 97 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
322 99 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?
331 
332  print *,'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
546  CALL init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu, &
548  ENDIF
549 #endif
550  ! ...................................................................
551 
552  !
553 334 FORMAT(1x,24i3)
554 755 FORMAT(1x,6f10.3,i3)
555 
556  RETURN
557  END SUBROUTINE inifilr
558 
559 END MODULE filtreg_mod
!$Header!c!c!c include serre h!c REAL alphax
Definition: serre.h:8
subroutine init_filtre_fft(coeffu, modfrstu, jfiltnu, jfiltsu, coeffv, modfrstv, jfiltnv, jfiltsv)
!$Id!COMMON coefils jfiltsu
Definition: coefils.h:4
real, dimension(:,:,:), allocatable matriceus
Definition: filtreg_mod.F90:6
logical, save use_filtre_fft
real, dimension(:,:,:), allocatable matrinvs
Definition: filtreg_mod.F90:7
!$Id!COMMON coefils eignfnv(iim, iim)&&
!$Id!COMMON coefils coefilu(iim, jjm)
!$Id!COMMON coefils eignfnu(iim, iim)
!$Id!COMMON coefils jfiltsv
Definition: coefils.h:4
subroutine inifilr
Definition: filtreg_mod.F90:12
!$Header!CDK comgeom COMMON comgeom rlatu
Definition: comgeom.h:25
!$Header!CDK comgeom2 COMMON comgeom unsaire xprimu
Definition: comgeom2.h:26
!$Id!COMMON coefils jfiltnu
Definition: coefils.h:4
subroutine inifgn(dv)
Definition: inifgn.F:5
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
!$Id!COMMON coefils && modfrstu(jjm)
!$Id!COMMON coefils coefilv(iim, jjm)
!$Id!COMMON coefils modfrstv(jjm)
!$Id fxyhypb
Definition: logic.h:10
!$Header jjp1
Definition: paramet.h:14
!$Header!CDK comgeom COMMON comgeom rlatv
Definition: comgeom.h:25
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
real, dimension(:,:,:), allocatable matrinvn
Definition: filtreg_mod.F90:7
!$Id!COMMON coefils jfiltnv
Definition: coefils.h:4
!$Id!COMMON coefils coefilu2(iim, jjm)
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
real, dimension(:,:,:), allocatable matriceun
Definition: filtreg_mod.F90:6
real, dimension(:,:,:), allocatable matricevn
Definition: filtreg_mod.F90:6
real, dimension(:,:,:), allocatable matricevs
Definition: filtreg_mod.F90:7