--- trunk/libf/filtrez/inifilr.f 2010/03/05 16:43:45 25 +++ trunk/libf/filtrez/inifilr.f90 2012/01/30 12:54:02 57 @@ -1,522 +1,416 @@ -! -! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifilr.F,v 1.1.1.1 2004/05/19 12:53:09 lmdzadmin Exp $ -! - SUBROUTINE inifilr -c -c ... H. Upadhyaya, O.Sharma ... -c - use dimens_m - use paramet_m - use logic - use comgeom - use serre - use parafilt - use coefils - IMPLICIT NONE -c -c version 3 ..... - -c Correction le 28/10/97 P. Le Van . -c ------------------------------------------------------------------- -c ------------------------------------------------------------------- - - REAL dlonu(iim),dlatu(jjm) - REAL rlamda( iim ), eignvl( iim ) -c - - REAL lamdamax,pi,cof - INTEGER i,j,modemax,imx,k,kf,ii - REAL dymin,dxmin,colat0 - REAL eignft(iim,iim), coff - REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs - COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus) - , , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs) - , , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus) - EXTERNAL inifgn -c -c ------------------------------------------------------------ -c This routine computes the eigenfunctions of the laplacien -c on the stretched grid, and the filtering coefficients -c -c We designate: -c eignfn eigenfunctions of the discrete laplacien -c eigenvl eigenvalues -c jfiltn indexof the last scalar line filtered in NH -c jfilts index of the first line filtered in SH -c modfrst index of the mode from where modes are filtered -c modemax maximum number of modes ( im ) -c coefil filtering coefficients ( lamda_max*cos(rlat)/lamda ) -c sdd SQRT( dx ) -c -c the modes are filtered from modfrst to modemax -c -c----------------------------------------------------------- -c - - pi = 2. * ASIN( 1. ) - - DO i = 1,iim - dlonu(i) = xprimu( i ) - ENDDO -c - CALL inifgn(eignvl) -c - print *,' EIGNVL ' - PRINT 250,eignvl -250 FORMAT( 1x,5e13.6) -c -c compute eigenvalues and eigenfunctions -c -c -c................................................................. -c -c compute the filtering coefficients for scalar lines and -c meridional wind v-lines -c -c we filter all those latitude lines where coefil < 1 -c NO FILTERING AT POLES -c -c colat0 is to be used when alpha (stretching coefficient) -c is set equal to zero for the regular grid case -c -c ....... Calcul de colat0 ......... -c ..... colat0 = minimum de ( 0.5, min dy/ min dx ) ... -c -c - DO 45 j = 1,jjm - dlatu( j ) = rlatu( j ) - rlatu( j+1 ) - 45 CONTINUE -c - dxmin = dlonu(1) - DO i = 2, iim - dxmin = MIN( dxmin,dlonu(i) ) - ENDDO - dymin = dlatu(1) - DO j = 2, jjm - dymin = MIN( dymin,dlatu(j) ) - ENDDO -c -c - colat0 = MIN( 0.5, dymin/dxmin ) -c - IF( .NOT.fxyhypb.AND.ysinus ) THEN - colat0 = 0.6 -c ...... a revoir pour ysinus ! ....... - alphax = 0. - ENDIF -c - PRINT 50, colat0,alphax - 50 FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7) -c - IF(alphax.EQ.1. ) THEN - PRINT *,' Inifilr alphax doit etre < a 1. Corriger ' - STOP 1 - ENDIF -c - lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) ) - -cc ... Correction le 28/10/97 ( P.Le Van ) .. -c - DO 71 i = 2,iim - rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) ) - 71 CONTINUE -c - - DO 72 j = 1,jjm - DO 73 i = 1,iim - coefilu( i,j ) = 0.0 - coefilv( i,j ) = 0.0 - coefilu2( i,j ) = 0.0 - coefilv2( i,j ) = 0.0 - 73 CONTINUE - 72 CONTINUE - -c -c ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv .... -c ......................................................... -c - modemax = iim - -cccc imx = modemax - 4 * (modemax/iim) - - imx = iim -c - PRINT *,' TRUNCATION AT ',imx -c - DO 75 j = 2, jjm/2+1 - cof = COS( rlatu(j) )/ colat0 - IF ( cof .LT. 1. ) THEN - IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j - ENDIF - - cof = COS( rlatu(jjp1-j+1) )/ colat0 - IF ( cof .LT. 1. ) THEN - IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) - $ jfiltsu= jjp1-j+1 - ENDIF - 75 CONTINUE -c - DO 76 j = 1, jjm/2 - cof = COS( rlatv(j) )/ colat0 - IF ( cof .LT. 1. ) THEN - IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j - ENDIF - - cof = COS( rlatv(jjm-j+1) )/ colat0 - IF ( cof .LT. 1. ) THEN - IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) - $ jfiltsv= jjm-j+1 - ENDIF - 76 CONTINUE -c - - if ( jfiltnu.LE.0 ) jfiltnu=1 - IF( jfiltnu.GT. jjm/2 +1 ) THEN - PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu - STOP 1 - ENDIF - - IF( jfiltsu.LE.0) jfiltsu=1 - IF( jfiltsu.GT. jjm +1 ) THEN - PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu - STOP 1 - ENDIF - - IF( jfiltnv.LE.0) jfiltnv=1 - IF( jfiltnv.GT. jjm/2 ) THEN - PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv - STOP 1 - ENDIF - - IF( jfiltsv.LE.0) jfiltsv=1 - IF( jfiltsv.GT. jjm ) THEN - PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv - STOP 1 - ENDIF - - PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , - * jfiltnv,jfiltsv,jfiltnu,jfiltsu - -c -c ... Determination de coefilu,coefilv,n=modfrstu,modfrstv .... -c................................................................ -c -c - DO 77 j = 1,jjm - modfrstu( j ) = iim - modfrstv( j ) = iim - 77 CONTINUE -c - DO 84 j = 2,jfiltnu - DO 81 k = 2,modemax - cof = rlamda(k) * COS( rlatu(j) ) - IF ( cof .LT. 1. ) GOTO 82 - 81 CONTINUE - GOTO 84 - 82 modfrstu( j ) = k -c - kf = modfrstu( j ) - DO 83 k = kf , modemax - cof = rlamda(k) * COS( rlatu(j) ) - coefilu(k,j) = cof - 1. - coefilu2(k,j) = cof*cof - 1. - 83 CONTINUE - 84 CONTINUE -c -c - DO 89 j = 1,jfiltnv -c - DO 86 k = 2,modemax - cof = rlamda(k) * COS( rlatv(j) ) - IF ( cof .LT. 1. ) GOTO 87 - 86 CONTINUE - GOTO 89 - 87 modfrstv( j ) = k -c - kf = modfrstv( j ) - DO 88 k = kf , modemax - cof = rlamda(k) * COS( rlatv(j) ) - coefilv(k,j) = cof - 1. - coefilv2(k,j) = cof*cof - 1. - 88 CONTINUE -c - 89 CONTINUE -c - DO 94 j = jfiltsu,jjm - DO 91 k = 2,modemax - cof = rlamda(k) * COS( rlatu(j) ) - IF ( cof .LT. 1. ) GOTO 92 - 91 CONTINUE - GOTO 94 - 92 modfrstu( j ) = k -c - kf = modfrstu( j ) - DO 93 k = kf , modemax - cof = rlamda(k) * COS( rlatu(j) ) - coefilu(k,j) = cof - 1. - coefilu2(k,j) = cof*cof - 1. - 93 CONTINUE - 94 CONTINUE -c - DO 99 j = jfiltsv,jjm - DO 96 k = 2,modemax - cof = rlamda(k) * COS( rlatv(j) ) - IF ( cof .LT. 1. ) GOTO 97 - 96 CONTINUE - GOTO 99 - 97 modfrstv( j ) = k -c - kf = modfrstv( j ) - DO 98 k = kf , modemax - cof = rlamda(k) * COS( rlatv(j) ) - coefilv(k,j) = cof - 1. - coefilv2(k,j) = cof*cof - 1. - 98 CONTINUE - 99 CONTINUE -c - - IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN - - IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv - IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu - - PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , - * jfiltnv,jfiltsv,jfiltnu,jfiltsu - ENDIF - - PRINT *,' Modes premiers v ' - PRINT 334,modfrstv - PRINT *,' Modes premiers u ' - PRINT 334,modfrstu - - - IF( nfilun.LT. 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 - ENDIF - IF( nfilun.GT. 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 -c STOP 1 - ENDIF - IF( nfilus.LT. 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 - ENDIF - IF( nfilus.GT. 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 -c STOP 1 - ENDIF - IF( nfilvn.LT. 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 - ENDIF - IF( nfilvn.GT. 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 -c STOP 1 - ENDIF - IF( nfilvs.LT. 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 - ENDIF - IF( nfilvs.GT. 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 -c STOP 1 - ENDIF - -c -c ................................................................... -c -c ... Calcul de la matrice filtre 'matriceu' pour les champs situes -c sur la grille scalaire ........ -c ................................................................... -c - DO j = 2, jfiltnu - - DO i=1,iim - coff = coefilu(i,j) - IF( i.LT.modfrstu(j) ) coff = 0. - DO k=1,iim - eignft(i,k) = eignfnv(k,i) * coff - ENDDO - ENDDO - 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) - ENDDO - ENDDO - ENDDO - - ENDDO - - DO j = jfiltsu, jjm - - DO i=1,iim - coff = coefilu(i,j) - IF( i.LT.modfrstu(j) ) coff = 0. - DO k=1,iim - eignft(i,k) = eignfnv(k,i) * coff - ENDDO - ENDDO - 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) - ENDDO - ENDDO - ENDDO - - ENDDO - -c ................................................................... -c -c ... Calcul de la matrice filtre 'matricev' pour les champs situes -c sur la grille de V ou de Z ........ -c ................................................................... -c - DO j = 1, jfiltnv - - DO i = 1, iim - coff = coefilv(i,j) - IF( i.LT.modfrstv(j) ) coff = 0. - DO k = 1, iim - eignft(i,k) = eignfnu(k,i) * coff - ENDDO - ENDDO - 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) - ENDDO - ENDDO - ENDDO - - ENDDO - - DO j = jfiltsv, jjm - - DO i = 1, iim - coff = coefilv(i,j) - IF( i.LT.modfrstv(j) ) coff = 0. - DO k = 1, iim - eignft(i,k) = eignfnu(k,i) * coff - ENDDO - ENDDO - 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) - ENDDO - ENDDO - ENDDO - - ENDDO - -c ................................................................... -c -c ... Calcul de la matrice filtre 'matrinv' pour les champs situes -c sur la grille scalaire , pour le filtre inverse ........ -c ................................................................... -c - DO j = 2, jfiltnu - - DO i = 1,iim - coff = coefilu(i,j)/ ( 1. + coefilu(i,j) ) - IF( i.LT.modfrstu(j) ) coff = 0. - DO k=1,iim - eignft(i,k) = eignfnv(k,i) * coff - ENDDO - ENDDO - 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) - ENDDO - ENDDO - ENDDO - - ENDDO - - DO j = jfiltsu, jjm - - DO i = 1,iim - coff = coefilu(i,j) / ( 1. + coefilu(i,j) ) - IF( i.LT.modfrstu(j) ) coff = 0. - DO k=1,iim - eignft(i,k) = eignfnv(k,i) * coff - ENDDO - ENDDO - 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) - ENDDO - ENDDO - ENDDO - - ENDDO - -c ................................................................... - -c -334 FORMAT(1x,24i3) -755 FORMAT(1x,6f10.3,i3) +module inifilr_m - RETURN - END + use dimens_m, only: iim + + IMPLICIT NONE + + INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv + INTEGER, PARAMETER:: nfilun=3, nfilus=2, nfilvn=2, nfilvs=2 + + real matriceun(iim,iim,nfilun), matriceus(iim,iim,nfilus) + real matricevn(iim,iim,nfilvn), matricevs(iim,iim,nfilvs) + real matrinvn(iim,iim,nfilun), matrinvs(iim,iim,nfilus) + + private iim, nfilun, nfilus, nfilvn, nfilvs + +contains + + SUBROUTINE inifilr + + ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09 + ! H. Upadhyaya, O. Sharma + + ! This routine computes the eigenfunctions of the laplacian on the + ! stretched grid, and the filtering coefficients. + ! We designate: + ! eignfn eigenfunctions of the discrete laplacian + ! eigenvl eigenvalues + ! jfiltn index of 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, ONLY : iim, jjm + use conf_gcm_m, ONLY : fxyhypb, ysinus + USE comgeom, ONLY : rlatu, rlatv, xprimu + use nr_util, only: pi + USE serre, ONLY : alphax + USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, & + eignfnv, modfrstu, modfrstv + + ! Local: + REAL dlonu(iim), dlatu(jjm) + REAL rlamda(2: iim), eignvl(iim) + + REAL lamdamax, cof + INTEGER i, j, modemax, imx, k, kf + REAL dymin, dxmin, colat0 + REAL eignft(iim, iim), coff + EXTERNAL inifgn + + !----------------------------------------------------------- + + print *, "Call sequence information: inifilr" + + DO i = 1, iim + dlonu(i) = xprimu(i) + END DO + + CALL inifgn(eignvl) + + PRINT *, 'EIGNVL ' + PRINT "(1X, 5E13.6)", eignvl + + ! 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 + + 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 + ! À revoir pour ysinus + alphax = 0. + END IF + + PRINT *, 'colat0 = ', colat0 + PRINT *, 'alphax = ', alphax + + IF (alphax == 1.) THEN + PRINT *, 'alphax doit etre < a 1. Corriger ' + STOP 1 + END IF + + lamdamax = iim / (pi * colat0 * (1. - alphax)) + rlamda = lamdamax / sqrt(abs(eignvl(2: iim))) + + DO j = 1, jjm + DO i = 1, iim + coefilu(i, j) = 0. + coefilv(i, j) = 0. + coefilu2(i, j) = 0. + coefilv2(i, j) = 0. + end DO + END DO + + ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv + + modemax = iim + imx = iim + + PRINT *, 'TRUNCATION AT ', imx + + DO j = 2, jjm / 2 + 1 + IF (cos(rlatu(j)) / colat0 < 1. & + .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j + + IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. & + .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) & + jfiltsu = jjm - j + 2 + 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.) exit + end DO + if (k == modemax + 1) cycle + 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.) exit + end DO + if (k == modemax + 1) cycle + 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.) exit + end DO + if (k == modemax + 1) cycle + 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.) exit + end DO + if (k == modemax + 1) cycle + 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 + IF (i < modfrstu(j)) then + coff = 0. + else + coff = coefilu(i, j) + end IF + eignft(i, :) = eignfnv(:, i)*coff + END DO + matriceun(:, :, j) = matmul(eignfnv, eignft) + END DO + + DO j = jfiltsu, jjm + DO i = 1, iim + IF (i < modfrstu(j)) then + coff = 0. + else + coff = coefilu(i, j) + end IF + eignft(i, :) = eignfnv(:, i) * coff + END DO + matriceus(:, :, j - jfiltsu + 1) = matmul(eignfnv, eignft) + 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 + IF (i < modfrstv(j)) then + coff = 0. + else + coff = coefilv(i, j) + end IF + eignft(i, :) = eignfnu(:, i)*coff + END DO + matricevn(:, :, j) = matmul(eignfnu, eignft) + END DO + + DO j = jfiltsv, jjm + DO i = 1, iim + IF (i < modfrstv(j)) then + coff = 0. + else + coff = coefilv(i, j) + end IF + eignft(i, :) = eignfnu(:, i)*coff + END DO + matricevs(:, :, j-jfiltsv+1) = matmul(eignfnu, eignft) + 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 + IF (i < modfrstu(j)) then + coff = 0. + else + coff = coefilu(i, j)/(1.+coefilu(i, j)) + end IF + eignft(i, :) = eignfnv(:, i)*coff + END DO + matrinvn(:, :, j) = matmul(eignfnv, eignft) + END DO + + DO j = jfiltsu, jjm + DO i = 1, iim + IF (i < modfrstu(j)) then + coff = 0. + else + coff = coefilu(i, j)/(1.+coefilu(i, j)) + end IF + eignft(i, :) = eignfnv(:, i)*coff + END DO + matrinvs(:, :, j-jfiltsu+1) = matmul(eignfnv, eignft) + END DO + +334 FORMAT (1X, 24I3) + + END SUBROUTINE inifilr + +end module inifilr_m