source: LMDZ4/branches/LMDZ4V5.0-LF/libf/filtrez/filtreg_mod.F90

Last change on this file was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

File size: 13.8 KB
Line 
1MODULE filtreg_mod
2
3  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
4  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
5
6CONTAINS
7
8  SUBROUTINE inifilr
9  USE mod_filtre_fft
10    !
11    !    ... H. Upadhyaya, O.Sharma   ...
12    !
13    IMPLICIT NONE
14    !
15    !     version 3 .....
16
17    !     Correction  le 28/10/97    P. Le Van .
18    !  -------------------------------------------------------------------
19#include "dimensions.h"
20#include "paramet.h"
21    !  -------------------------------------------------------------------
22#include "comgeom.h"
23#include "coefils.h"
24#include "logic.h"
25#include "serre.h"
26
27    REAL  dlonu(iim),dlatu(jjm)
28    REAL  rlamda( iim ),  eignvl( iim )
29    !
30
31    REAL    lamdamax,pi,cof
32    INTEGER i,j,modemax,imx,k,kf,ii
33    REAL dymin,dxmin,colat0
34    REAL eignft(iim,iim), coff
35
36    LOGICAL, SAVE :: first_call_inifilr = .TRUE.
37
38#ifdef CRAY
39    INTEGER   ISMIN
40    EXTERNAL  ISMIN
41    INTEGER iymin
42    INTEGER ixmineq
43#endif
44    EXTERNAL  inifgn
45    !
46    ! ------------------------------------------------------------
47    !   This routine computes the eigenfunctions of the laplacien
48    !   on the stretched grid, and the filtering coefficients
49    !     
50    !  We designate:
51    !   eignfn   eigenfunctions of the discrete laplacien
52    !   eigenvl  eigenvalues
53    !   jfiltn   indexof the last scalar line filtered in NH
54    !   jfilts   index of the first line filtered in SH
55    !   modfrst  index of the mode from WHERE modes are filtered
56    !   modemax  maximum number of modes ( im )
57    !   coefil   filtering coefficients ( lamda_max*COS(rlat)/lamda )
58    !   sdd      SQRT( dx )
59    !     
60    !     the modes are filtered from modfrst to modemax
61    !     
62    !-----------------------------------------------------------
63    !
64
65    pi       = 2. * ASIN( 1. )
66
67    DO i = 1,iim
68       dlonu(i) = xprimu( i )
69    ENDDO
70    !
71    CALL inifgn(eignvl)
72    !
73    PRINT *,' EIGNVL '
74    PRINT 250,eignvl
75250 FORMAT( 1x,5e13.6)
76    !
77    ! compute eigenvalues and eigenfunctions
78    !
79    !
80    !.................................................................
81    !
82    !  compute the filtering coefficients for scalar lines and
83    !  meridional wind v-lines
84    !
85    !  we filter all those latitude lines WHERE coefil < 1
86    !  NO FILTERING AT POLES
87    !
88    !  colat0 is to be used  when alpha (stretching coefficient)
89    !  is set equal to zero for the regular grid CASE
90    !
91    !    .......   Calcul  de  colat0   .........
92    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
93    !
94    !
95    DO j = 1,jjm
96       dlatu( j ) = rlatu( j ) - rlatu( j+1 )
97    ENDDO
98    !
99#ifdef CRAY
100    iymin   = ISMIN( jjm, dlatu, 1 )
101    ixmineq = ISMIN( iim, dlonu, 1 )
102    dymin   = dlatu( iymin )
103    dxmin   = dlonu( ixmineq )
104#else
105    dxmin   =  dlonu(1)
106    DO  i  = 2, iim
107       dxmin = MIN( dxmin,dlonu(i) )
108    ENDDO
109    dymin  = dlatu(1)
110    DO j  = 2, jjm
111       dymin = MIN( dymin,dlatu(j) )
112    ENDDO
113#endif
114    !
115    !
116    colat0  =  MIN( 0.5, dymin/dxmin )
117    !
118    IF( .NOT.fxyhypb.AND.ysinus )  THEN
119       colat0 = 0.6
120       !         ...... a revoir  pour  ysinus !   .......
121       alphax = 0.
122    ENDIF
123    !
124    PRINT 50, colat0,alphax
12550  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
126    !
127    IF(alphax.EQ.1. )  THEN
128       PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
129       STOP
130    ENDIF
131    !
132    lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
133
134    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
135    !
136    DO i = 2,iim
137       rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
138    ENDDO
139    !
140
141    DO j = 1,jjm
142       DO i = 1,iim
143          coefilu( i,j )  = 0.0
144          coefilv( i,j )  = 0.0
145          coefilu2( i,j ) = 0.0
146          coefilv2( i,j ) = 0.0
147       ENDDO
148    ENDDO
149
150    !
151    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
152    !    .........................................................
153    !
154    modemax = iim
155
156!!!!    imx = modemax - 4 * (modemax/iim)
157
158    imx  = iim
159    !
160    PRINT *,' TRUNCATION AT ',imx
161    !
162    DO j = 2, jjm/2+1
163       cof = COS( rlatu(j) )/ colat0
164       IF ( cof .LT. 1. ) THEN
165          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
166       ENDIF
167
168       cof = COS( rlatu(jjp1-j+1) )/ colat0
169       IF ( cof .LT. 1. ) THEN
170          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
171               jfiltsu= jjp1-j+1
172       ENDIF
173    ENDDO
174    !
175    DO j = 1, jjm/2
176       cof = COS( rlatv(j) )/ colat0
177       IF ( cof .LT. 1. ) THEN
178          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
179       ENDIF
180
181       cof = COS( rlatv(jjm-j+1) )/ colat0
182       IF ( cof .LT. 1. ) THEN
183          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
184               jfiltsv= jjm-j+1
185       ENDIF
186    ENDDO
187    !                                 
188
189    IF ( jfiltnu.LE.0 ) jfiltnu=1
190    IF( jfiltnu.GT. jjm/2 +1 )  THEN
191       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
192       STOP
193    ENDIF
194
195    IF( jfiltsu.LE.0) jfiltsu=1
196    IF( jfiltsu.GT.  jjm  +1 )  THEN
197       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
198       STOP
199    ENDIF
200
201    IF( jfiltnv.LE.0) jfiltnv=1
202    IF( jfiltnv.GT. jjm/2    )  THEN
203       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
204       STOP
205    ENDIF
206
207    IF( jfiltsv.LE.0) jfiltsv=1
208    IF( jfiltsv.GT.     jjm  )  THEN
209       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
210       STOP
211    ENDIF
212
213    PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
214         jfiltnv,jfiltsv,jfiltnu,jfiltsu
215
216    IF(first_call_inifilr) THEN
217       ALLOCATE(matriceun(iim,iim,jfiltnu))
218       ALLOCATE(matriceus(iim,iim,jfiltsu))
219       ALLOCATE(matricevn(iim,iim,jfiltnv))
220       ALLOCATE(matricevs(iim,iim,jfiltsv))
221       ALLOCATE( matrinvn(iim,iim,jfiltnu))
222       ALLOCATE( matrinvs(iim,iim,jfiltsu))
223       first_call_inifilr = .FALSE.
224    ENDIF
225
226    !                                 
227    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
228    !................................................................
229    !
230    !
231    DO j = 1,jjm
232       modfrstu( j ) = iim
233       modfrstv( j ) = iim
234    ENDDO
235    !
236    DO j = 2,jfiltnu
237       DO k = 2,modemax
238          cof = rlamda(k) * COS( rlatu(j) )
239          IF ( cof .LT. 1. ) GOTO 82
240       ENDDO
241       GOTO 84
24282     modfrstu( j ) = k
243       !
244       kf = modfrstu( j )
245       DO k = kf , modemax
246          cof = rlamda(k) * COS( rlatu(j) )
247          coefilu(k,j) = cof - 1.
248          coefilu2(k,j) = cof*cof - 1.
249       ENDDO
25084     CONTINUE
251    ENDDO
252    !                                 
253    !
254    DO j = 1,jfiltnv
255       !
256       DO k = 2,modemax
257          cof = rlamda(k) * COS( rlatv(j) )
258          IF ( cof .LT. 1. ) GOTO 87
259       ENDDO
260       GOTO 89
26187     modfrstv( j ) = k
262       !
263       kf = modfrstv( j )
264       DO k = kf , modemax
265          cof = rlamda(k) * COS( rlatv(j) )
266          coefilv(k,j) = cof - 1.
267          coefilv2(k,j) = cof*cof - 1.
268       ENDDO
26989     CONTINUE
270    ENDDO
271    !
272    DO j = jfiltsu,jjm
273       DO k = 2,modemax
274          cof = rlamda(k) * COS( rlatu(j) )
275          IF ( cof .LT. 1. ) GOTO 92
276       ENDDO
277       GOTO 94
27892     modfrstu( j ) = k
279       !
280       kf = modfrstu( j )
281       DO k = kf , modemax
282          cof = rlamda(k) * COS( rlatu(j) )
283          coefilu(k,j) = cof - 1.
284          coefilu2(k,j) = cof*cof - 1.
285       ENDDO
28694     CONTINUE
287    ENDDO
288    !                                 
289    DO j = jfiltsv,jjm
290       DO k = 2,modemax
291          cof = rlamda(k) * COS( rlatv(j) )
292          IF ( cof .LT. 1. ) GOTO 97
293       ENDDO
294       GOTO 99
29597     modfrstv( j ) = k
296       !
297       kf = modfrstv( j )
298       DO k = kf , modemax
299          cof = rlamda(k) * COS( rlatv(j) )
300          coefilv(k,j) = cof - 1.
301          coefilv2(k,j) = cof*cof - 1.
302       ENDDO
30399     CONTINUE
304    ENDDO
305    !
306
307    IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
308
309       IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
310       IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
311
312       PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
313            jfiltnv,jfiltsv,jfiltnu,jfiltsu
314    ENDIF
315
316    PRINT *,'   Modes premiers  v  '
317    PRINT 334,modfrstv
318    PRINT *,'   Modes premiers  u  '
319    PRINT 334,modfrstu
320
321    ! 
322    !   ...................................................................
323    !
324    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
325    !                       sur la grille scalaire                 ........
326    !   ...................................................................
327    !
328    DO j = 2, jfiltnu
329
330       DO i=1,iim
331          coff = coefilu(i,j)
332          IF( i.LT.modfrstu(j) ) coff = 0.
333          DO k=1,iim
334             eignft(i,k) = eignfnv(k,i) * coff
335          ENDDO
336       ENDDO
337#ifdef CRAY
338       CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
339#else
340#ifdef BLAS
341       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
342            eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
343#else
344       DO k = 1, iim
345          DO i = 1, iim
346             matriceun(i,k,j) = 0.0
347             DO ii = 1, iim
348                matriceun(i,k,j) = matriceun(i,k,j) &
349                     + eignfnv(i,ii)*eignft(ii,k)
350             ENDDO
351          ENDDO
352       ENDDO
353#endif
354#endif
355
356    ENDDO
357
358    DO j = jfiltsu, jjm
359
360       DO i=1,iim
361          coff = coefilu(i,j)
362          IF( i.LT.modfrstu(j) ) coff = 0.
363          DO k=1,iim
364             eignft(i,k) = eignfnv(k,i) * coff
365          ENDDO
366       ENDDO
367#ifdef CRAY
368       CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
369#else
370#ifdef BLAS
371       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
372            eignfnv, iim, eignft, iim, 0.0, &
373            matriceus(1,1,j-jfiltsu+1), iim)
374#else
375       DO k = 1, iim
376          DO i = 1, iim
377             matriceus(i,k,j-jfiltsu+1) = 0.0
378             DO ii = 1, iim
379                matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
380                     + eignfnv(i,ii)*eignft(ii,k)
381             ENDDO
382          ENDDO
383       ENDDO
384#endif
385#endif
386
387    ENDDO
388
389    !   ...................................................................
390    !
391    !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
392    !                       sur la grille   de V ou de Z           ........
393    !   ...................................................................
394    !
395    DO j = 1, jfiltnv
396
397       DO i = 1, iim
398          coff = coefilv(i,j)
399          IF( i.LT.modfrstv(j) ) coff = 0.
400          DO k = 1, iim
401             eignft(i,k) = eignfnu(k,i) * coff
402          ENDDO
403       ENDDO
404#ifdef CRAY
405       CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
406#else
407#ifdef BLAS
408       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
409            eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
410#else
411       DO k = 1, iim
412          DO i = 1, iim
413             matricevn(i,k,j) = 0.0
414             DO ii = 1, iim
415                matricevn(i,k,j) = matricevn(i,k,j) &
416                     + eignfnu(i,ii)*eignft(ii,k)
417             ENDDO
418          ENDDO
419       ENDDO
420#endif
421#endif
422
423    ENDDO
424
425    DO j = jfiltsv, jjm
426
427       DO i = 1, iim
428          coff = coefilv(i,j)
429          IF( i.LT.modfrstv(j) ) coff = 0.
430          DO k = 1, iim
431             eignft(i,k) = eignfnu(k,i) * coff
432          ENDDO
433       ENDDO
434#ifdef CRAY
435       CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
436#else
437#ifdef BLAS
438       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
439            eignfnu, iim, eignft, iim, 0.0, &
440            matricevs(1,1,j-jfiltsv+1), iim)
441#else
442       DO k = 1, iim
443          DO i = 1, iim
444             matricevs(i,k,j-jfiltsv+1) = 0.0
445             DO ii = 1, iim
446                matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
447                     + eignfnu(i,ii)*eignft(ii,k)
448             ENDDO
449          ENDDO
450       ENDDO
451#endif
452#endif
453
454    ENDDO
455
456    !   ...................................................................
457    !
458    !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
459    !              sur la grille scalaire , pour le filtre inverse ........
460    !   ...................................................................
461    !
462    DO j = 2, jfiltnu
463
464       DO i = 1,iim
465          coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
466          IF( i.LT.modfrstu(j) ) coff = 0.
467          DO k=1,iim
468             eignft(i,k) = eignfnv(k,i) * coff
469          ENDDO
470       ENDDO
471#ifdef CRAY
472       CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
473#else
474#ifdef BLAS
475       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
476            eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
477#else
478       DO k = 1, iim
479          DO i = 1, iim
480             matrinvn(i,k,j) = 0.0
481             DO ii = 1, iim
482                matrinvn(i,k,j) = matrinvn(i,k,j) &
483                     + eignfnv(i,ii)*eignft(ii,k)
484             ENDDO
485          ENDDO
486       ENDDO
487#endif
488#endif
489
490    ENDDO
491
492    DO j = jfiltsu, jjm
493
494       DO i = 1,iim
495          coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
496          IF( i.LT.modfrstu(j) ) coff = 0.
497          DO k=1,iim
498             eignft(i,k) = eignfnv(k,i) * coff
499          ENDDO
500       ENDDO
501#ifdef CRAY
502       CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
503#else
504#ifdef BLAS
505       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
506            eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
507#else
508       DO k = 1, iim
509          DO i = 1, iim
510             matrinvs(i,k,j-jfiltsu+1) = 0.0
511             DO ii = 1, iim
512                matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
513                     + eignfnv(i,ii)*eignft(ii,k)
514             ENDDO
515          ENDDO
516       ENDDO
517#endif
518#endif
519
520    ENDDO
521
522    IF (use_filtre_fft) THEN
523       CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,  &
524                           coefilv,modfrstv,jfiltnv,jfiltsv)
525    ENDIF
526
527    !   ...................................................................
528
529    !
530334 FORMAT(1x,24i3)
531755 FORMAT(1x,6f10.3,i3)
532
533    RETURN
534  END SUBROUTINE inifilr
535
536END MODULE filtreg_mod
Note: See TracBrowser for help on using the repository browser.