Sources/filtrez/inifilr.f

SUBROUTINE inifilr

  ! From filtrez/inifilr.F,v 1.1.1.1 2004/05/19 12:53:09
  ! H. Upadhyaya, O.Sharma

  !   This routine computes the eigenfunctions of the laplacien           
  !   on the stretched grid, and the filtering coefficients               

  !  We designate:                                                        
  !   eignfn   eigenfunctions of the discrete laplacien                   
  !   eigenvl  eigenvalues                                                
  !   jfiltn   indexof the last scalar line filtered in NH                
  !   jfilts   index of the first line filtered in SH                     
  !   modfrst  index of the mode from where modes are filtered            
  !   modemax  maximum number of modes ( im )                             
  !   coefil   filtering coefficients ( lamda_max*cos(rlat)/lamda )       
  !   sdd      SQRT( dx )                                                 

  !     the modes are filtered from modfrst to modemax                    

  USE dimens_m
  USE paramet_m
  USE logic
  USE comgeom
  USE serre
  USE parafilt
  USE coefils

  IMPLICIT NONE

  REAL dlonu(iim), dlatu(jjm)
  REAL rlamda(iim), eignvl(iim)

  REAL lamdamax, pi, cof
  INTEGER i, j, modemax, imx, k, kf, ii
  REAL dymin, dxmin, colat0
  REAL eignft(iim,iim), coff
  EXTERNAL inifgn

  !-----------------------------------------------------------            


  pi = 2.*asin(1.)

  DO i = 1, iim
     dlonu(i) = xprimu(i)
  END DO

  CALL inifgn(eignvl)

  PRINT *, ' EIGNVL '
  PRINT 250, eignvl
250 FORMAT (1X,5E13.6)

  ! compute eigenvalues and eigenfunctions                                


  !.................................................................      

  !  compute the filtering coefficients for scalar lines and              
  !  meridional wind v-lines                                              

  !  we filter all those latitude lines where coefil < 1                  
  !  NO FILTERING AT POLES                                                

  !  colat0 is to be used  when alpha (stretching coefficient)            
  !  is set equal to zero for the regular grid case                       

  !    .......   Calcul  de  colat0   .........                           
  !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...          


  DO  j = 1, jjm
     dlatu(j) = rlatu(j) - rlatu(j+1)
  END DO

  dxmin = dlonu(1)
  DO i = 2, iim
     dxmin = min(dxmin,dlonu(i))
  END DO
  dymin = dlatu(1)
  DO j = 2, jjm
     dymin = min(dymin,dlatu(j))
  END DO


  colat0 = min(0.5,dymin/dxmin)

  IF ( .NOT. fxyhypb .AND. ysinus) THEN
     colat0 = 0.6
     !         ...... a revoir  pour  ysinus !   .......                     
     alphax = 0.
  END IF

  PRINT 50, colat0, alphax
50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)

  IF (alphax==1.) THEN
     PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '
     STOP 1
  END IF

  lamdamax = iim/(pi*colat0*(1.-alphax))

  DO  i = 2, iim
     rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))
  END DO


  DO  j = 1, jjm
     DO  i = 1, iim
        coefilu(i,j) = 0.0
        coefilv(i,j) = 0.0
        coefilu2(i,j) = 0.0
        coefilv2(i,j) = 0.0
     end DO
  END DO


  !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....          
  !    .........................................................          

  modemax = iim

  !ccc    imx = modemax - 4 * (modemax/iim)                               

  imx = iim

  PRINT *, ' TRUNCATION AT ', imx

  DO  j = 2, jjm/2 + 1
     cof = cos(rlatu(j))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j
     END IF

     cof = cos(rlatu(jjp1-j+1))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1
     END IF
  END DO

  DO  j = 1, jjm/2
     cof = cos(rlatv(j))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j
     END IF

     cof = cos(rlatv(jjm-j+1))/colat0
     IF (cof<1.) THEN
        IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1
     END IF
  END DO


  IF (jfiltnu<=0) jfiltnu = 1
  IF (jfiltnu>jjm/2+1) THEN
     PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
     STOP 1
  END IF

  IF (jfiltsu<=0) jfiltsu = 1
  IF (jfiltsu>jjm+1) THEN
     PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
     STOP 1
  END IF

  IF (jfiltnv<=0) jfiltnv = 1
  IF (jfiltnv>jjm/2) THEN
     PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
     STOP 1
  END IF

  IF (jfiltsv<=0) jfiltsv = 1
  IF (jfiltsv>jjm) THEN
     PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
     STOP 1
  END IF

  PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
       jfiltsu


  !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....       
  !................................................................       


  DO  j = 1, jjm
     modfrstu(j) = iim
     modfrstv(j) = iim
  END DO

  DO  j = 2, jfiltnu
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatu(j))
        IF (cof<1.) GO TO 82
     end DO
     cycle
82   modfrstu(j) = k

     kf = modfrstu(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatu(j))
        coefilu(k,j) = cof - 1.
        coefilu2(k,j) = cof*cof - 1.
     end DO
  END DO


  DO j = 1, jfiltnv
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatv(j))
        IF (cof<1.) GO TO 87
     end DO
     cycle
87   modfrstv(j) = k

     kf = modfrstv(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatv(j))
        coefilv(k,j) = cof - 1.
        coefilv2(k,j) = cof*cof - 1.
     end DO
  end DO

  DO  j = jfiltsu, jjm
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatu(j))
        IF (cof<1.) GO TO 92
     end DO
     cycle
92   modfrstu(j) = k

     kf = modfrstu(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatu(j))
        coefilu(k,j) = cof - 1.
        coefilu2(k,j) = cof*cof - 1.
     end DO
  end DO

  DO  j = jfiltsv, jjm
     DO  k = 2, modemax
        cof = rlamda(k)*cos(rlatv(j))
        IF (cof<1.) GO TO 97
     end DO
     cycle
97   modfrstv(j) = k

     kf = modfrstv(j)
     DO  k = kf, modemax
        cof = rlamda(k)*cos(rlatv(j))
        coefilv(k,j) = cof - 1.
        coefilv2(k,j) = cof*cof - 1.
     end DO
  END DO


  IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN

     IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv
     IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu

     PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
          jfiltsu
  END IF

  PRINT *, '   Modes premiers  v  '
  PRINT 334, modfrstv
  PRINT *, '   Modes premiers  u  '
  PRINT 334, modfrstu


  IF (nfilun<jfiltnu) THEN
     PRINT *, ' le parametre nfilun utilise pour la matrice ', &
          ' matriceun  est trop petit ! '
     PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu
     PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
     STOP 1
  END IF
  IF (nfilun>jfiltnu+2) THEN
     PRINT *, ' le parametre nfilun utilise pour la matrice ', &
          ' matriceun est trop grand ! Gachis de memoire ! '
     PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnu
     PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
  END IF
  IF (nfilus<jjm-jfiltsu+1) THEN
     PRINT *, ' le parametre nfilus utilise pour la matrice ', &
          ' matriceus  est trop petit !  '
     PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &
          jjm - jfiltsu + 1
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
     STOP 1
  END IF
  IF (nfilus>jjm-jfiltsu+3) THEN
     PRINT *, ' le parametre nfilus utilise pour la matrice ', &
          ' matriceus  est trop grand ! '
     PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &
          jjm - jfiltsu + 1
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
  END IF
  IF (nfilvn<jfiltnv) THEN
     PRINT *, ' le parametre nfilvn utilise pour la matrice ', &
          ' matricevn  est trop petit ! '
     PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
     STOP 1
  END IF
  IF (nfilvn>jfiltnv+2) THEN
     PRINT *, ' le parametre nfilvn utilise pour la matrice ', &
          ' matricevn est trop grand !  Gachis de memoire ! '
     PRINT *, 'Le changer dans parafilt.h et le mettre a  ', jfiltnv
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
  END IF
  IF (nfilvs<jjm-jfiltsv+1) THEN
     PRINT *, ' le parametre nfilvs utilise pour la matrice ', &
          ' matricevs  est trop petit !  Le changer dans parafilt.h '
     PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &
          jjm - jfiltsv + 1
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
     STOP 1
  END IF
  IF (nfilvs>jjm-jfiltsv+3) THEN
     PRINT *, ' le parametre nfilvs utilise pour la matrice ', &
          ' matricevs  est trop grand ! Gachis de memoire ! '
     PRINT *, ' Le changer dans parafilt.h et le mettre a  ', &
          jjm - jfiltsv + 1
     PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
          'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
          jfiltnv, jjm - jfiltsv + 1
  END IF

  !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes  
  !                       sur la grille scalaire                 ........ 

  DO j = 2, jfiltnu

     DO i = 1, iim
        coff = coefilu(i,j)
        IF (i<modfrstu(j)) coff = 0.
        DO k = 1, iim
           eignft(i,k) = eignfnv(k,i)*coff
        END DO
     END DO
     DO k = 1, iim
        DO i = 1, iim
           matriceun(i,k,j) = 0.0
           DO ii = 1, iim
              matriceun(i,k,j) = matriceun(i,k,j) + eignfnv(i,ii)*eignft(ii,k)
           END DO
        END DO
     END DO

  END DO

  DO j = jfiltsu, jjm

     DO i = 1, iim
        coff = coefilu(i,j)
        IF (i<modfrstu(j)) coff = 0.
        DO k = 1, iim
           eignft(i,k) = eignfnv(k,i)*coff
        END DO
     END DO
     DO k = 1, iim
        DO i = 1, iim
           matriceus(i,k,j-jfiltsu+1) = 0.0
           DO ii = 1, iim
              matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) + &
                   eignfnv(i,ii)*eignft(ii,k)
           END DO
        END DO
     END DO

  END DO

  !   ................................................................... 

  !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes  
  !                       sur la grille   de V ou de Z           ........ 
  !   ................................................................... 

  DO j = 1, jfiltnv

     DO i = 1, iim
        coff = coefilv(i,j)
        IF (i<modfrstv(j)) coff = 0.
        DO k = 1, iim
           eignft(i,k) = eignfnu(k,i)*coff
        END DO
     END DO
     DO k = 1, iim
        DO i = 1, iim
           matricevn(i,k,j) = 0.0
           DO ii = 1, iim
              matricevn(i,k,j) = matricevn(i,k,j) + eignfnu(i,ii)*eignft(ii,k)
           END DO
        END DO
     END DO

  END DO

  DO j = jfiltsv, jjm

     DO i = 1, iim
        coff = coefilv(i,j)
        IF (i<modfrstv(j)) coff = 0.
        DO k = 1, iim
           eignft(i,k) = eignfnu(k,i)*coff
        END DO
     END DO
     DO k = 1, iim
        DO i = 1, iim
           matricevs(i,k,j-jfiltsv+1) = 0.0
           DO ii = 1, iim
              matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) + &
                   eignfnu(i,ii)*eignft(ii,k)
           END DO
        END DO
     END DO

  END DO

  !   ................................................................... 

  !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes   
  !              sur la grille scalaire , pour le filtre inverse ........ 
  !   ................................................................... 

  DO j = 2, jfiltnu

     DO i = 1, iim
        coff = coefilu(i,j)/(1.+coefilu(i,j))
        IF (i<modfrstu(j)) coff = 0.
        DO k = 1, iim
           eignft(i,k) = eignfnv(k,i)*coff
        END DO
     END DO
     DO k = 1, iim
        DO i = 1, iim
           matrinvn(i,k,j) = 0.0
           DO ii = 1, iim
              matrinvn(i,k,j) = matrinvn(i,k,j) + eignfnv(i,ii)*eignft(ii,k)
           END DO
        END DO
     END DO

  END DO

  DO j = jfiltsu, jjm

     DO i = 1, iim
        coff = coefilu(i,j)/(1.+coefilu(i,j))
        IF (i<modfrstu(j)) coff = 0.
        DO k = 1, iim
           eignft(i,k) = eignfnv(k,i)*coff
        END DO
     END DO
     DO k = 1, iim
        DO i = 1, iim
           matrinvs(i,k,j-jfiltsu+1) = 0.0
           DO ii = 1, iim
              matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) + &
                   eignfnv(i,ii)*eignft(ii,k)
           END DO
        END DO
     END DO

  END DO

334 FORMAT (1X,24I3)

END SUBROUTINE inifilr
1	guez	32	SUBROUTINE inifilr
2	guez	3
3	guez	32	! From filtrez/inifilr.F,v 1.1.1.1 2004/05/19 12:53:09
4			! H. Upadhyaya, O.Sharma
5	guez	3
6	guez	32	! This routine computes the eigenfunctions of the laplacien
7			! on the stretched grid, and the filtering coefficients
8	guez	3
9	guez	32	! We designate:
10			! eignfn eigenfunctions of the discrete laplacien
11			! eigenvl eigenvalues
12			! jfiltn indexof the last scalar line filtered in NH
13			! jfilts index of the first line filtered in SH
14			! modfrst index of the mode from where modes are filtered
15			! modemax maximum number of modes ( im )
16			! coefil filtering coefficients ( lamda_max*cos(rlat)/lamda )
17			! sdd SQRT( dx )
18	guez	3
19	guez	32	! the modes are filtered from modfrst to modemax
20	guez	3
21	guez	32	USE dimens_m
22			USE paramet_m
23			USE logic
24			USE comgeom
25			USE serre
26			USE parafilt
27			USE coefils
28	guez	3
29	guez	32	IMPLICIT NONE
30	guez	3
31	guez	32	REAL dlonu(iim), dlatu(jjm)
32			REAL rlamda(iim), eignvl(iim)
33	guez	3
34	guez	32	REAL lamdamax, pi, cof
35			INTEGER i, j, modemax, imx, k, kf, ii
36			REAL dymin, dxmin, colat0
37			REAL eignft(iim,iim), coff
38			EXTERNAL inifgn
39	guez	3
40	guez	32	!-----------------------------------------------------------
41	guez	3
42
43	guez	32	pi = 2.*asin(1.)
44	guez	3
45	guez	32	DO i = 1, iim
46			dlonu(i) = xprimu(i)
47			END DO
48	guez	3
49	guez	32	CALL inifgn(eignvl)
50	guez	3
51	guez	32	PRINT *, ' EIGNVL '
52			PRINT 250, eignvl
53			250 FORMAT (1X,5E13.6)
54	guez	3
55	guez	32	! compute eigenvalues and eigenfunctions
56	guez	3
57
58	guez	32	!.................................................................
59	guez	3
60	guez	32	! compute the filtering coefficients for scalar lines and
61			! meridional wind v-lines
62	guez	3
63	guez	32	! we filter all those latitude lines where coefil < 1
64			! NO FILTERING AT POLES
65	guez	3
66	guez	32	! colat0 is to be used when alpha (stretching coefficient)
67			! is set equal to zero for the regular grid case
68	guez	3
69	guez	32	! ....... Calcul de colat0 .........
70			! ..... colat0 = minimum de ( 0.5, min dy/ min dx ) ...
71	guez	3
72
73	guez	32	DO j = 1, jjm
74			dlatu(j) = rlatu(j) - rlatu(j+1)
75			END DO
76	guez	3
77	guez	32	dxmin = dlonu(1)
78			DO i = 2, iim
79			dxmin = min(dxmin,dlonu(i))
80			END DO
81			dymin = dlatu(1)
82			DO j = 2, jjm
83			dymin = min(dymin,dlatu(j))
84			END DO
85	guez	3
86
87	guez	32	colat0 = min(0.5,dymin/dxmin)
88	guez	3
89	guez	32	IF ( .NOT. fxyhypb .AND. ysinus) THEN
90			colat0 = 0.6
91			! ...... a revoir pour ysinus ! .......
92			alphax = 0.
93			END IF
94	guez	3
95	guez	32	PRINT 50, colat0, alphax
96			50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)
97	guez	3
98	guez	32	IF (alphax==1.) THEN
99			PRINT *, ' Inifilr alphax doit etre < a 1. Corriger '
100			STOP 1
101			END IF
102	guez	3
103	guez	32	lamdamax = iim/(picolat0(1.-alphax))
104	guez	3
105	guez	32	DO i = 2, iim
106			rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))
107			END DO
108	guez	3
109
110	guez	32	DO j = 1, jjm
111			DO i = 1, iim
112			coefilu(i,j) = 0.0
113			coefilv(i,j) = 0.0
114			coefilu2(i,j) = 0.0
115			coefilv2(i,j) = 0.0
116			end DO
117			END DO
118	guez	3
119
120	guez	32	! ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
121			! .........................................................
122	guez	3
123	guez	32	modemax = iim
124	guez	3
125	guez	32	!ccc imx = modemax - 4 * (modemax/iim)
126	guez	3
127	guez	32	imx = iim
128	guez	3
129	guez	32	PRINT *, ' TRUNCATION AT ', imx
130	guez	3
131	guez	32	DO j = 2, jjm/2 + 1
132			cof = cos(rlatu(j))/colat0
133			IF (cof<1.) THEN
134			IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j
135			END IF
136	guez	3
137	guez	32	cof = cos(rlatu(jjp1-j+1))/colat0
138			IF (cof<1.) THEN
139			IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1
140			END IF
141			END DO
142	guez	3
143	guez	32	DO j = 1, jjm/2
144			cof = cos(rlatv(j))/colat0
145			IF (cof<1.) THEN
146			IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j
147			END IF
148	guez	3
149	guez	32	cof = cos(rlatv(jjm-j+1))/colat0
150			IF (cof<1.) THEN
151			IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1
152			END IF
153			END DO
154	guez	3
155	guez	32
156			IF (jfiltnu<=0) jfiltnu = 1
157			IF (jfiltnu>jjm/2+1) THEN
158			PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
159			STOP 1
160			END IF
161
162			IF (jfiltsu<=0) jfiltsu = 1
163			IF (jfiltsu>jjm+1) THEN
164			PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
165			STOP 1
166			END IF
167
168			IF (jfiltnv<=0) jfiltnv = 1
169			IF (jfiltnv>jjm/2) THEN
170			PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
171			STOP 1
172			END IF
173
174			IF (jfiltsv<=0) jfiltsv = 1
175			IF (jfiltsv>jjm) THEN
176			PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
177			STOP 1
178			END IF
179
180			PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
181			jfiltsu
182
183
184			! ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
185			!................................................................
186
187
188			DO j = 1, jjm
189			modfrstu(j) = iim
190			modfrstv(j) = iim
191			END DO
192
193			DO j = 2, jfiltnu
194			DO k = 2, modemax
195			cof = rlamda(k)*cos(rlatu(j))
196			IF (cof<1.) GO TO 82
197			end DO
198			cycle
199			82 modfrstu(j) = k
200
201			kf = modfrstu(j)
202			DO k = kf, modemax
203			cof = rlamda(k)*cos(rlatu(j))
204			coefilu(k,j) = cof - 1.
205			coefilu2(k,j) = cof*cof - 1.
206			end DO
207			END DO
208
209
210			DO j = 1, jfiltnv
211			DO k = 2, modemax
212			cof = rlamda(k)*cos(rlatv(j))
213			IF (cof<1.) GO TO 87
214			end DO
215			cycle
216			87 modfrstv(j) = k
217
218			kf = modfrstv(j)
219			DO k = kf, modemax
220			cof = rlamda(k)*cos(rlatv(j))
221			coefilv(k,j) = cof - 1.
222			coefilv2(k,j) = cof*cof - 1.
223			end DO
224			end DO
225
226			DO j = jfiltsu, jjm
227			DO k = 2, modemax
228			cof = rlamda(k)*cos(rlatu(j))
229			IF (cof<1.) GO TO 92
230			end DO
231			cycle
232			92 modfrstu(j) = k
233
234			kf = modfrstu(j)
235			DO k = kf, modemax
236			cof = rlamda(k)*cos(rlatu(j))
237			coefilu(k,j) = cof - 1.
238			coefilu2(k,j) = cof*cof - 1.
239			end DO
240			end DO
241
242			DO j = jfiltsv, jjm
243			DO k = 2, modemax
244			cof = rlamda(k)*cos(rlatv(j))
245			IF (cof<1.) GO TO 97
246			end DO
247			cycle
248			97 modfrstv(j) = k
249
250			kf = modfrstv(j)
251			DO k = kf, modemax
252			cof = rlamda(k)*cos(rlatv(j))
253			coefilv(k,j) = cof - 1.
254			coefilv2(k,j) = cof*cof - 1.
255			end DO
256			END DO
257
258
259			IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
260
261			IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv
262			IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu
263
264			PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
265			jfiltsu
266			END IF
267
268			PRINT *, ' Modes premiers v '
269			PRINT 334, modfrstv
270			PRINT *, ' Modes premiers u '
271			PRINT 334, modfrstu
272
273
274			IF (nfilun<jfiltnu) THEN
275			PRINT *, ' le parametre nfilun utilise pour la matrice ', &
276			' matriceun est trop petit ! '
277			PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu
278			PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &
279			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
280			jfiltnv, jjm - jfiltsv + 1
281			STOP 1
282			END IF
283			IF (nfilun>jfiltnu+2) THEN
284			PRINT *, ' le parametre nfilun utilise pour la matrice ', &
285			' matriceun est trop grand ! Gachis de memoire ! '
286			PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu
287			PRINT *, 'Pour information, nfilun,nfilus,nfilvn,nfilvs ', &
288			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
289			jfiltnv, jjm - jfiltsv + 1
290			END IF
291			IF (nfilus<jjm-jfiltsu+1) THEN
292			PRINT *, ' le parametre nfilus utilise pour la matrice ', &
293			' matriceus est trop petit ! '
294			PRINT *, ' Le changer dans parafilt.h et le mettre a ', &
295			jjm - jfiltsu + 1
296			PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
297			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
298			jfiltnv, jjm - jfiltsv + 1
299			STOP 1
300			END IF
301			IF (nfilus>jjm-jfiltsu+3) THEN
302			PRINT *, ' le parametre nfilus utilise pour la matrice ', &
303			' matriceus est trop grand ! '
304			PRINT *, ' Le changer dans parafilt.h et le mettre a ', &
305			jjm - jfiltsu + 1
306			PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
307			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
308			jfiltnv, jjm - jfiltsv + 1
309			END IF
310			IF (nfilvn<jfiltnv) THEN
311			PRINT *, ' le parametre nfilvn utilise pour la matrice ', &
312			' matricevn est trop petit ! '
313			PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv
314			PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
315			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
316			jfiltnv, jjm - jfiltsv + 1
317			STOP 1
318			END IF
319			IF (nfilvn>jfiltnv+2) THEN
320			PRINT *, ' le parametre nfilvn utilise pour la matrice ', &
321			' matricevn est trop grand ! Gachis de memoire ! '
322			PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv
323			PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
324			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
325			jfiltnv, jjm - jfiltsv + 1
326			END IF
327			IF (nfilvs<jjm-jfiltsv+1) THEN
328			PRINT *, ' le parametre nfilvs utilise pour la matrice ', &
329			' matricevs est trop petit ! Le changer dans parafilt.h '
330			PRINT *, ' Le changer dans parafilt.h et le mettre a ', &
331			jjm - jfiltsv + 1
332			PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
333			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
334			jfiltnv, jjm - jfiltsv + 1
335			STOP 1
336			END IF
337			IF (nfilvs>jjm-jfiltsv+3) THEN
338			PRINT *, ' le parametre nfilvs utilise pour la matrice ', &
339			' matricevs est trop grand ! Gachis de memoire ! '
340			PRINT *, ' Le changer dans parafilt.h et le mettre a ', &
341			jjm - jfiltsv + 1
342			PRINT *, ' Pour information , nfilun,nfilus,nfilvn,nfilvs ', &
343			'doivent etre egaux successivement a ', jfiltnu, jjm - jfiltsu + 1, &
344			jfiltnv, jjm - jfiltsv + 1
345			END IF
346
347			! ... Calcul de la matrice filtre 'matriceu' pour les champs situes
348			! sur la grille scalaire ........
349
350			DO j = 2, jfiltnu
351
352			DO i = 1, iim
353			coff = coefilu(i,j)
354			IF (i<modfrstu(j)) coff = 0.
355			DO k = 1, iim
356			eignft(i,k) = eignfnv(k,i)*coff
357			END DO
358			END DO
359			DO k = 1, iim
360			DO i = 1, iim
361			matriceun(i,k,j) = 0.0
362			DO ii = 1, iim
363			matriceun(i,k,j) = matriceun(i,k,j) + eignfnv(i,ii)*eignft(ii,k)
364			END DO
365			END DO
366			END DO
367
368			END DO
369
370			DO j = jfiltsu, jjm
371
372			DO i = 1, iim
373			coff = coefilu(i,j)
374			IF (i<modfrstu(j)) coff = 0.
375			DO k = 1, iim
376			eignft(i,k) = eignfnv(k,i)*coff
377			END DO
378			END DO
379			DO k = 1, iim
380			DO i = 1, iim
381			matriceus(i,k,j-jfiltsu+1) = 0.0
382			DO ii = 1, iim
383			matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) + &
384			eignfnv(i,ii)*eignft(ii,k)
385			END DO
386			END DO
387			END DO
388
389			END DO
390
391			! ...................................................................
392
393			! ... Calcul de la matrice filtre 'matricev' pour les champs situes
394			! sur la grille de V ou de Z ........
395			! ...................................................................
396
397			DO j = 1, jfiltnv
398
399			DO i = 1, iim
400			coff = coefilv(i,j)
401			IF (i<modfrstv(j)) coff = 0.
402			DO k = 1, iim
403			eignft(i,k) = eignfnu(k,i)*coff
404			END DO
405			END DO
406			DO k = 1, iim
407			DO i = 1, iim
408			matricevn(i,k,j) = 0.0
409			DO ii = 1, iim
410			matricevn(i,k,j) = matricevn(i,k,j) + eignfnu(i,ii)*eignft(ii,k)
411			END DO
412			END DO
413			END DO
414
415			END DO
416
417			DO j = jfiltsv, jjm
418
419			DO i = 1, iim
420			coff = coefilv(i,j)
421			IF (i<modfrstv(j)) coff = 0.
422			DO k = 1, iim
423			eignft(i,k) = eignfnu(k,i)*coff
424			END DO
425			END DO
426			DO k = 1, iim
427			DO i = 1, iim
428			matricevs(i,k,j-jfiltsv+1) = 0.0
429			DO ii = 1, iim
430			matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) + &
431			eignfnu(i,ii)*eignft(ii,k)
432			END DO
433			END DO
434			END DO
435
436			END DO
437
438			! ...................................................................
439
440			! ... Calcul de la matrice filtre 'matrinv' pour les champs situes
441			! sur la grille scalaire , pour le filtre inverse ........
442			! ...................................................................
443
444			DO j = 2, jfiltnu
445
446			DO i = 1, iim
447			coff = coefilu(i,j)/(1.+coefilu(i,j))
448			IF (i<modfrstu(j)) coff = 0.
449			DO k = 1, iim
450			eignft(i,k) = eignfnv(k,i)*coff
451			END DO
452			END DO
453			DO k = 1, iim
454			DO i = 1, iim
455			matrinvn(i,k,j) = 0.0
456			DO ii = 1, iim
457			matrinvn(i,k,j) = matrinvn(i,k,j) + eignfnv(i,ii)*eignft(ii,k)
458			END DO
459			END DO
460			END DO
461
462			END DO
463
464			DO j = jfiltsu, jjm
465
466			DO i = 1, iim
467			coff = coefilu(i,j)/(1.+coefilu(i,j))
468			IF (i<modfrstu(j)) coff = 0.
469			DO k = 1, iim
470			eignft(i,k) = eignfnv(k,i)*coff
471			END DO
472			END DO
473			DO k = 1, iim
474			DO i = 1, iim
475			matrinvs(i,k,j-jfiltsu+1) = 0.0
476			DO ii = 1, iim
477			matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) + &
478			eignfnv(i,ii)*eignft(ii,k)
479			END DO
480			END DO
481			END DO
482
483			END DO
484
485			334 FORMAT (1X,24I3)
486
487			END SUBROUTINE inifilr