source: LMDZ6/branches/blowing_snow/libf/filtrez/filtreg_mod.F90 @ 4933

Last change on this file since 4933 was 4440, checked in by fhourdin, 17 months ago

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