! $Id$ MODULE lmdz_filtreg USE lmdz_paramet 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_dimensions, ONLY: iim, jjm, llm, ndm 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 ! 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 USE lmdz_ssum_scopy, ONLY: ssum USE lmdz_comgeom USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE 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 ! 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) USE lmdz_ssum_scopy, ONLY: ssum IMPLICIT NONE INTEGER :: im REAL :: vec(im, im), d(im) INTEGER :: i, j REAL :: sum 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 lmdz_filtre_fft, ONLY: use_filtre_fft,Init_filtre_fft USE lmdz_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 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_comgeom USE lmdz_paramet ! ... H. Upfiltreg_modadhyaya, O.Sharma ... ! version 3 ..... ! Correction le 28/10/97 P. Le Van . ! ------------------------------------------------------------------- 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