! $Id$ MODULE lmdz_filtreg IMPLICIT NONE; PRIVATE PUBLIC matriceun, matriceus, matricevn, matricevs, matrinvn, matrinvs, & inifilr, filtreg REAL, DIMENSION(:, :, :), ALLOCATABLE :: matriceun, matriceus, matricevn REAL, DIMENSION(:, :, :), ALLOCATABLE :: matricevs, matrinvn, matrinvs CONTAINS SUBROUTINE filtreg(champ, nlat, nbniv, ifiltre, iaire, & griscal, iter) USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu !======================================================================= ! ! Auteur: P. Le Van 07/10/97 ! ------ ! ! Objet: filtre matriciel longitudinal ,avec les matrices precalculees ! pour l'operateur Filtre . ! ------ ! ! Arguments: ! ---------- ! ! nblat nombre de latitudes a filtrer ! nbniv nombre de niveaux verticaux a filtrer ! champ(iip1,nblat,nbniv) en entree : champ a filtrer ! en sortie : champ filtre ! ifiltre +1 Transformee directe ! -1 Transformee inverse ! +2 Filtre directe ! -2 Filtre inverse ! ! iaire 1 si champ intensif ! 2 si champ extensif (pondere par les aires) ! ! iter 1 filtre simple ! !======================================================================= ! ! ! Variable Intensive ! ifiltre = 1 filtre directe ! ifiltre =-1 filtre inverse ! ! Variable Extensive ! ifiltre = 2 filtre directe ! ifiltre =-2 filtre inverse ! ! INCLUDE "dimensions.h" INCLUDE "paramet.h" INTEGER :: nlat, nbniv, ifiltre, iter INTEGER :: i, j, l, k INTEGER :: iim2, immjm INTEGER :: jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil REAL :: champ(iip1, nlat, nbniv) REAL :: eignq(iim, nlat, nbniv), sdd1(iim), sdd2(iim) LOGICAL :: griscal INTEGER :: hemisph, iaire LOGICAL, SAVE :: first = .TRUE. REAL, SAVE :: sdd12(iim, 4) INTEGER, PARAMETER :: type_sddu = 1 INTEGER, PARAMETER :: type_sddv = 2 INTEGER, PARAMETER :: type_unsddu = 3 INTEGER, PARAMETER :: type_unsddv = 4 INTEGER :: sdd1_type, sdd2_type IF (iim == 1) return ! no filtre in 2D y-z IF (first) THEN sdd12(1:iim, type_sddu) = sddu(1:iim) sdd12(1:iim, type_sddv) = sddv(1:iim) sdd12(1:iim, type_unsddu) = unsddu(1:iim) sdd12(1:iim, type_unsddv) = unsddv(1:iim) first = .FALSE. ENDIF IF(ifiltre==1.OR.ifiltre==-1) & stop 'Pas de transformee simple dans cette version' IF(iter== 2) THEN PRINT *, ' Pas d iteration du filtre dans cette version !'& &, ' Utiliser old_filtreg et repasser !' STOP ENDIF IF(ifiltre== -2 .AND..NOT.griscal) THEN PRINT *, ' Cette routine ne calcule le filtre inverse que ' & , ' sur la grille des scalaires !' STOP ENDIF IF(ifiltre/=2 .AND.ifiltre/= - 2) THEN PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' & , ' corriger et repasser !' STOP ENDIF iim2 = iim * iim immjm = iim * jjm IF(griscal) THEN IF(nlat /= jjp1) THEN PRINT 1111 STOP ELSE IF(iaire==1) THEN sdd1_type = type_sddv sdd2_type = type_unsddv ELSE sdd1_type = type_unsddv sdd2_type = type_sddv ENDIF ! IF( iaire.EQ.1 ) THEN ! CALL SCOPY( iim, sddv, 1, sdd1, 1 ) ! CALL SCOPY( iim, unsddv, 1, sdd2, 1 ) ! ELSE ! CALL SCOPY( iim, unsddv, 1, sdd1, 1 ) ! CALL SCOPY( iim, sddv, 1, sdd2, 1 ) ! END IF jdfil1 = 2 jffil1 = jfiltnu jdfil2 = jfiltsu jffil2 = jjm END IF ELSE IF(nlat/=jjm) THEN PRINT 2222 STOP ELSE IF(iaire==1) THEN sdd1_type = type_sddu sdd2_type = type_unsddu ELSE sdd1_type = type_unsddu sdd2_type = type_sddu ENDIF ! IF( iaire.EQ.1 ) THEN ! CALL SCOPY( iim, sddu, 1, sdd1, 1 ) ! CALL SCOPY( iim, unsddu, 1, sdd2, 1 ) ! ELSE ! CALL SCOPY( iim, unsddu, 1, sdd1, 1 ) ! CALL SCOPY( iim, sddu, 1, sdd2, 1 ) ! END IF jdfil1 = 1 jffil1 = jfiltnv jdfil2 = jfiltsv jffil2 = jjm END IF END IF DO hemisph = 1, 2 IF (hemisph==1) THEN jdfil = jdfil1 jffil = jffil1 ELSE jdfil = jdfil2 jffil = jffil2 END IF DO l = 1, nbniv DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = champ(i, j, l) * sdd12(i, sdd1_type) ! sdd1(i) END DO END DO END DO IF(hemisph == 1) THEN IF(ifiltre == -2) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matrinvn(1,1,j), & iim, champ(1,j,1), iip1*nlat, 0.0, & eignq(1,j-jdfil+1,1), iim*nlat) #else eignq(:, j - jdfil + 1, :) & = matmul(matrinvn(:, :, j), champ(:iim, j, :)) #endif END DO ELSE IF (griscal) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matriceun(1,1,j), & iim, champ(1,j,1), iip1*nlat, 0.0, & eignq(1,j-jdfil+1,1), iim*nlat) #else eignq(:, j - jdfil + 1, :) & = matmul(matriceun(:, :, j), champ(:iim, j, :)) #endif END DO ELSE DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matricevn(1,1,j), & iim, champ(1,j,1), iip1*nlat, 0.0, & eignq(1,j-jdfil+1,1), iim*nlat) #else eignq(:, j - jdfil + 1, :) & = matmul(matricevn(:, :, j), champ(:iim, j, :)) #endif END DO ENDIF ELSE IF(ifiltre == -2) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matrinvs(1,1,j-jfiltsu+1), & iim, champ(1,j,1), iip1*nlat, 0.0, & eignq(1,j-jdfil+1,1), iim*nlat) #else eignq(:, j - jdfil + 1, :) & = matmul(matrinvs(:, :, j - jfiltsu + 1), & champ(:iim, j, :)) #endif END DO ELSE IF (griscal) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matriceus(1,1,j-jfiltsu+1), & iim, champ(1,j,1), iip1*nlat, 0.0, & eignq(1,j-jdfil+1,1), iim*nlat) #else eignq(:, j - jdfil + 1, :) & = matmul(matriceus(:, :, j - jfiltsu + 1), & champ(:iim, j, :)) #endif END DO ELSE DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matricevs(1,1,j-jfiltsv+1), & iim, champ(1,j,1), iip1*nlat, 0.0, & eignq(1,j-jdfil+1,1), iim*nlat) #else eignq(:, j - jdfil + 1, :) & = matmul(matricevs(:, :, j - jfiltsv + 1), & champ(:iim, j, :)) #endif END DO ENDIF ENDIF IF(ifiltre== 2) THEN DO l = 1, nbniv DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = & (champ(i, j, l) + eignq(i, j - jdfil + 1, l)) & * sdd12(i, sdd2_type) ! sdd2(i) END DO END DO END DO ELSE DO l = 1, nbniv DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = & (champ(i, j, l) - eignq(i, j - jdfil + 1, l)) & * sdd12(i, sdd2_type) ! sdd2(i) END DO END DO END DO ENDIF DO l = 1, nbniv DO j = jdfil, jffil champ(iip1, j, l) = champ(1, j, l) END DO END DO ENDDO 1111 FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau CHAMP a& & filtrer, sur la grille des scalaires'/) 2222 FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau CHAMP a fi& & ltrer, sur la grille de V ou de Z'/) RETURN END SUBROUTINE filtreg SUBROUTINE inifgn(dv) ! ! ... H.Upadyaya , O.Sharma ... ! USE lmdz_coefils, ONLY: sddv, sddu, unsddu, unsddv, eignfnv, eignfnu IMPLICIT NONE ! include "dimensions.h" include "paramet.h" include "comgeom.h" ! REAL :: vec(iim, iim), vec1(iim, iim) REAL :: dlonu(iim), dlonv(iim) REAL :: du(iim), dv(iim), d(iim) REAL :: pi INTEGER :: i, j, k, imm1, nrot EXTERNAL SSUM REAL :: SSUM ! imm1 = iim - 1 pi = 2. * ASIN(1.) ! DO i = 1, iim dlonu(i) = xprimu(i) dlonv(i) = xprimv(i) END DO DO i = 1, iim sddv(i) = SQRT(dlonv(i)) sddu(i) = SQRT(dlonu(i)) unsddu(i) = 1. / sddu(i) unsddv(i) = 1. / sddv(i) END DO ! DO j = 1, iim DO i = 1, iim vec(i, j) = 0. vec1(i, j) = 0. eignfnv(i, j) = 0. eignfnu(i, j) = 0. END DO END DO ! ! eignfnv(1, 1) = -1. eignfnv(iim, 1) = 1. DO i = 1, imm1 eignfnv(i + 1, i + 1) = -1. eignfnv(i, i + 1) = 1. END DO DO j = 1, iim DO i = 1, iim eignfnv(i, j) = eignfnv(i, j) / (sddu(i) * sddv(j)) END DO END DO DO j = 1, iim DO i = 1, iim eignfnu(i, j) = -eignfnv(j, i) END DO END DO ! DO j = 1, iim DO i = 1, iim vec (i, j) = 0.0 vec1(i, j) = 0.0 DO k = 1, iim vec (i, j) = vec(i, j) + eignfnu(i, k) * eignfnv(k, j) vec1(i, j) = vec1(i, j) + eignfnv(i, k) * eignfnu(k, j) ENDDO ENDDO ENDDO ! CALL jacobi(vec, iim, iim, dv, eignfnv, nrot) CALL acc(eignfnv, d, iim) CALL eigen_sort(dv, eignfnv, iim, iim) ! CALL jacobi(vec1, iim, iim, du, eignfnu, nrot) CALL acc(eignfnu, d, iim) CALL eigen_sort(du, eignfnu, iim, iim) !c ancienne version avec appels IMSL ! ! CALL MXM(eignfnu,iim,eignfnv,iim,vec,iim) ! CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim) ! CALL EVCSF(iim,vec,iim,dv,eignfnv,iim) ! CALL acc(eignfnv,d,iim) ! CALL eigen(eignfnv,dv) ! ! CALL EVCSF(iim,vec1,iim,du,eignfnu,iim) ! CALL acc(eignfnu,d,iim) ! CALL eigen(eignfnu,du) RETURN END SUBROUTINE inifgn SUBROUTINE JACOBI(A, N, NP, D, V, NROT) IMPLICIT NONE ! Arguments: INTEGER, INTENT(IN) :: N INTEGER, INTENT(IN) :: NP INTEGER, INTENT(OUT) :: NROT REAL, INTENT(INOUT) :: A(NP, NP) REAL, INTENT(OUT) :: D(NP) REAL, INTENT(OUT) :: V(NP, NP) ! local variables: INTEGER :: IP, IQ, I, J REAL :: SM, TRESH, G, H, T, THETA, C, S, TAU REAL :: B(N) REAL :: Z(N) DO IP = 1, N DO IQ = 1, N V(IP, IQ) = 0. ENDDO V(IP, IP) = 1. ENDDO DO IP = 1, N B(IP) = A(IP, IP) D(IP) = B(IP) Z(IP) = 0. ENDDO NROT = 0 DO I = 1, 50 ! 50? I suspect this should be NP ! but convergence is fast enough anyway SM = 0. DO IP = 1, N - 1 DO IQ = IP + 1, N SM = SM + ABS(A(IP, IQ)) ENDDO ENDDO IF(SM==0.)RETURN IF(I<4)THEN TRESH = 0.2 * SM / N**2 ELSE TRESH = 0. ENDIF DO IP = 1, N - 1 DO IQ = IP + 1, N G = 100. * ABS(A(IP, IQ)) IF((I>4).AND.(ABS(D(IP)) + G==ABS(D(IP))) & .AND.(ABS(D(IQ)) + G==ABS(D(IQ))))THEN A(IP, IQ) = 0. ELSE IF(ABS(A(IP, IQ))>TRESH)THEN H = D(IQ) - D(IP) IF(ABS(H) + G==ABS(H))THEN T = A(IP, IQ) / H ELSE THETA = 0.5 * H / A(IP, IQ) T = 1. / (ABS(THETA) + SQRT(1. + THETA**2)) IF(THETA<0.)T = -T ENDIF C = 1. / SQRT(1 + T**2) S = T * C TAU = S / (1. + C) H = T * A(IP, IQ) Z(IP) = Z(IP) - H Z(IQ) = Z(IQ) + H D(IP) = D(IP) - H D(IQ) = D(IQ) + H A(IP, IQ) = 0. DO J = 1, IP - 1 G = A(J, IP) H = A(J, IQ) A(J, IP) = G - S * (H + G * TAU) A(J, IQ) = H + S * (G - H * TAU) ENDDO DO J = IP + 1, IQ - 1 G = A(IP, J) H = A(J, IQ) A(IP, J) = G - S * (H + G * TAU) A(J, IQ) = H + S * (G - H * TAU) ENDDO DO J = IQ + 1, N G = A(IP, J) H = A(IQ, J) A(IP, J) = G - S * (H + G * TAU) A(IQ, J) = H + S * (G - H * TAU) ENDDO DO J = 1, N G = V(J, IP) H = V(J, IQ) V(J, IP) = G - S * (H + G * TAU) V(J, IQ) = H + S * (G - H * TAU) ENDDO NROT = NROT + 1 ENDIF ENDDO ENDDO DO IP = 1, N B(IP) = B(IP) + Z(IP) D(IP) = B(IP) Z(IP) = 0. ENDDO ENDDO ! of DO I=1,50 STOP 'Jacobi: 50 iterations should never happen' END SUBROUTINE JACOBI SUBROUTINE eigen_sort(d, v, n, np) INTEGER :: n, np REAL :: d(np), v(np, np) INTEGER :: i, j, k REAL :: p DO i = 1, n - 1 k = i p = d(i) DO j = i + 1, n IF(d(j)>=p) THEN k = j p = d(j) ENDIF ENDDO IF(k/=i) THEN d(k) = d(i) d(i) = p DO j = 1, n p = v(j, i) v(j, i) = v(j, k) v(j, k) = p ENDDO ENDIF ENDDO RETURN END SUBROUTINE eigen_sort SUBROUTINE acc(vec, d, im) IMPLICIT NONE INTEGER :: im REAL :: vec(im, im), d(im) INTEGER :: i, j REAL :: sum REAL, external :: ssum do j = 1, im do i = 1, im d(i) = vec(i, j) * vec(i, j) enddo sum = ssum(im, d, 1) sum = sqrt(sum) do i = 1, im vec(i, j) = vec(i, j) / sum enddo enddo RETURN END SUBROUTINE acc SUBROUTINE inifilr #ifdef CPP_PARA USE mod_filtre_fft, ONLY: use_filtre_fft,Init_filtre_fft USE mod_filtre_fft_loc, ONLY: Init_filtre_fft_loc=>Init_filtre_fft ! #endif USE serre_mod, ONLY: alphax USE logic_mod, ONLY: fxyhypb, ysinus USE comconst_mod, ONLY: maxlatfilter USE lmdz_coefils, ONLY: modfrstv, modfrstu, jfiltnu, jfiltnv, coefilu, coefilv, & coefilu2, coefilv2, eignfnv, eignfnu, jfiltsu, jfiltsv ! ... H. Upfiltreg_modadhyaya, O.Sharma ... ! version 3 ..... ! Correction le 28/10/97 P. Le Van . ! ------------------------------------------------------------------- include "dimensions.h" include "paramet.h" include "comgeom.h" 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 LOGICAL, SAVE :: first_call_inifilr = .TRUE. INTEGER ISMIN EXTERNAL ISMIN INTEGER iymin INTEGER ixmineq ! ------------------------------------------------------------ ! 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 !----------------------------------------------------------- IF (iim == 1) return ! No filtre in 2D y-z pi = 2. * ASIN(1.) DO i = 1, iim dlonu(i) = xprimu(i) ENDDO CALL inifgn(eignvl) PRINT *, 'inifilr: EIGNVL ' PRINT 250, eignvl 250 FORMAT(1x, 5e14.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) ENDDO 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 ! For a regular grid, we want the filter to start at latitudes ! corresponding to lengths dx of the same size as dy (in terms ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees ! <=> latitude=60 degrees). ! Same idea for the zoomed grid: start filtering polewards as soon ! as length dx becomes of the same size as dy ! if maxlatfilter >0, prescribe the colat0 value from the .def files IF (maxlatfilter < 0.) THEN colat0 = MIN(0.5, dymin / dxmin) ! colat0 = 1. IF(.NOT.fxyhypb.AND.ysinus) THEN colat0 = 0.6 ! ...... a revoir pour ysinus ! ....... alphax = 0. ENDIF ELSE colat0 = (90.0 - maxlatfilter) / 180.0 * pi ENDIF 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 ENDIF lamdamax = iim / (pi * colat0 * (1. - alphax)) ! ... Correction le 28/10/97 ( P.Le Van ) .. DO i = 2, iim rlamda(i) = lamdamax / SQRT(ABS(eignvl(i))) ENDDO 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 ENDDO ENDDO ! ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv .... ! ......................................................... modemax = iim !!!! imx = modemax - 4 * (modemax/iim) imx = iim PRINT *, 'inifilr: TRUNCATION AT ', imx ! Ehouarn: set up some defaults jfiltnu = 2 ! avoid north pole jfiltsu = jjm ! avoid south pole (which is at jjm+1) jfiltnv = 1 ! NB: no poles on the V grid jfiltsv = jjm DO j = 2, jjm / 2 + 1 cof = COS(rlatu(j)) / colat0 IF (cof < 1.) THEN IF(rlamda(imx) * COS(rlatu(j))<1.) THEN jfiltnu = j ENDIF ENDIF cof = COS(rlatu(jjp1 - j + 1)) / colat0 IF (cof < 1.) THEN IF(rlamda(imx) * COS(rlatu(jjp1 - j + 1))<1.) THEN jfiltsu = jjp1 - j + 1 ENDIF ENDIF ENDDO DO j = 1, jjm / 2 cof = COS(rlatv(j)) / colat0 IF (cof < 1.) THEN IF(rlamda(imx) * COS(rlatv(j))<1.) THEN jfiltnv = j ENDIF ENDIF cof = COS(rlatv(jjm - j + 1)) / colat0 IF (cof < 1.) THEN IF(rlamda(imx) * COS(rlatv(jjm - j + 1))<1.) THEN jfiltsv = jjm - j + 1 ENDIF ENDIF ENDDO IF(jfiltnu> jjm / 2 + 1) THEN PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu STOP ENDIF IF(jfiltsu> jjm + 1) THEN PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu STOP ENDIF IF(jfiltnv> jjm / 2) THEN PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv STOP ENDIF IF(jfiltsv> jjm) THEN PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv STOP ENDIF PRINT *, 'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ', & jfiltnv, jfiltsv, jfiltnu, jfiltsu IF(first_call_inifilr) THEN ALLOCATE(matriceun(iim, iim, jfiltnu)) ALLOCATE(matriceus(iim, iim, jjm - jfiltsu + 1)) ALLOCATE(matricevn(iim, iim, jfiltnv)) ALLOCATE(matricevs(iim, iim, jjm - jfiltsv + 1)) ALLOCATE(matrinvn(iim, iim, jfiltnu)) ALLOCATE(matrinvs(iim, iim, jjm - jfiltsu + 1)) first_call_inifilr = .FALSE. ENDIF ! ... Determination de coefilu,coefilv,n=modfrstu,modfrstv .... !................................................................ DO j = 1, jjm !default initialization: all modes are retained (i.e. no filtering) modfrstu(j) = iim modfrstv(j) = iim ENDDO DO j = 2, jfiltnu DO k = 2, modemax cof = rlamda(k) * COS(rlatu(j)) IF (cof < 1.) GOTO 82 ENDDO GOTO 84 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. ENDDO 84 CONTINUE ENDDO DO j = 1, jfiltnv DO k = 2, modemax cof = rlamda(k) * COS(rlatv(j)) IF (cof < 1.) GOTO 87 ENDDO GOTO 89 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. ENDDO 89 CONTINUE ENDDO DO j = jfiltsu, jjm DO k = 2, modemax cof = rlamda(k) * COS(rlatu(j)) IF (cof < 1.) GOTO 92 ENDDO GOTO 94 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. ENDDO 94 CONTINUE ENDDO DO j = jfiltsv, jjm DO k = 2, modemax cof = rlamda(k) * COS(rlatv(j)) IF (cof < 1.) GOTO 97 ENDDO GOTO 99 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. ENDDO 99 CONTINUE ENDDO IF(jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2)THEN ! Ehouarn: and what are these for??? Trying to handle a limit case ! where filters extend to and meet at the equator? IF(jfiltnv==jfiltsv)jfiltsv = 1 + jfiltnv IF(jfiltnu==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 ! ................................................................... ! ... 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