SUBROUTINE inifilr c c ... H. Upadhyaya, O.Sharma ... c IMPLICIT NONE c c version 3 ..... c Correction le 28/10/97 P. Le Van . c ------------------------------------------------------------------- #include "dimensions.h" #include "paramet.h" #include "parafilt.h" c ------------------------------------------------------------------- #include "comgeom.h" #include "coefils.h" #include "logic.h" #include "serre.h" 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) #ifdef CRAY INTEGER ISMIN EXTERNAL ISMIN INTEGER iymin INTEGER ixmineq #endif 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 #ifdef CRAY iymin = ISMIN( jjm, dlatu, 1 ) ixmineq = ISMIN( iim, dlonu, 1 ) dymin = dlatu( iymin ) dxmin = dlonu( ixmineq ) #else 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 #endif c c PRINT *, 'TB15 : dymin, dxmin =',dymin,dxmin 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 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 c Forcing filters: TB22 this was added in some version jfiltnu=2 jfiltsu=2 DO 75 j = 2, jjm/2+1 cof = COS( rlatu(j) )/ colat0 !PRINT *,' j,cof,expr= ',j,cof,rlamda(imx)*COS(rlatu(j)) 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 PRINT *,'TB22 jfiltsu ' ,jfiltsu IF( jfiltnu.LE.0 .OR. jfiltnu.GT. jjm/2 +1 ) THEN PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu STOP ENDIF IF( jfiltsu.LE.0 .OR. jfiltsu.GT. jjm +1 ) THEN PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu STOP ENDIF IF( jfiltnv.LE.0 .OR. jfiltnv.GT. jjm/2 ) THEN PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv STOP ENDIF IF( jfiltsv.LE.0 .OR. jfiltsv.GT. jjm ) THEN PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv STOP 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 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 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 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 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 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 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 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 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 #ifdef CRAY CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim ) #else #ifdef BLAS CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, $ eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim) #else 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 #endif #endif 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 #ifdef CRAY CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim) #else #ifdef BLAS CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, $ eignfnv, iim, eignft, iim, 0.0, $ matriceus(1,1,j-jfiltsu+1), iim) #else 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 #endif #endif 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 #ifdef CRAY CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim ) #else #ifdef BLAS CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, $ eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim) #else 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 #endif #endif 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 #ifdef CRAY CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim) #else #ifdef BLAS CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, $ eignfnu, iim, eignft, iim, 0.0, $ matricevs(1,1,j-jfiltsv+1), iim) #else 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 #endif #endif 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 #ifdef CRAY CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim ) #else #ifdef BLAS CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, $ eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim) #else 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 #endif #endif 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 #ifdef CRAY CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim) #else #ifdef BLAS CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, $ eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim) #else 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 #endif #endif ENDDO c ................................................................... c 334 FORMAT(1x,24i3) 755 FORMAT(1x,6f10.3,i3) RETURN END