/[lmdze]/trunk/filtrez/inifilr.f
ViewVC logotype

Diff of /trunk/filtrez/inifilr.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/filtrez/inifilr.f revision 5 by guez, Mon Mar 3 16:32:04 2008 UTC trunk/Sources/filtrez/inifilr.f revision 143 by guez, Tue Jun 9 14:32:46 2015 UTC
# Line 1  Line 1 
1  !  module inifilr_m
 ! $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  
       IMPLICIT NONE  
 c  
 c     version 3 .....  
   
 c     Correction  le 28/10/97    P. Le Van .  
 c  -------------------------------------------------------------------  
 c  -------------------------------------------------------------------  
       include "coefils.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)  
       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)  
2    
3         RETURN    IMPLICIT NONE
4         END  
5      INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6      ! jfiltn index of the last scalar line filtered in NH
7      ! jfilts index of the first line filtered in SH
8    
9      ! North:
10      real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
11      ! (iim, iim, 2:jfiltnu)
12    
13      real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
14    
15      ! South:
16      real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
17      ! (iim, iim, jfiltsu:jjm)
18    
19      real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
20    
21    contains
22    
23      SUBROUTINE inifilr
24    
25        ! From filtrez/inifilr.F, version 1.1.1.1 2004/05/19 12:53:09
26        ! H. Upadhyaya, O. Sharma
27    
28        ! This routine computes the eigenfunctions of the laplacian on the
29        ! stretched grid, and the filtering coefficients. The modes are
30        ! filtered from modfrst to iim.
31    
32        USE dimens_m, ONLY : iim, jjm
33        USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
34        use inifgn_m, only: inifgn, eignfnu, eignfnv
35        use jumble, only: new_unit
36        use nr_util, only: pi
37    
38        ! Local:
39        REAL dlatu(jjm)
40        REAL rlamda(2: iim)
41        real eignvl(iim) ! eigenvalues
42        REAL cof
43        INTEGER i, j, k, kf, unit
44        REAL colat0
45        REAL eignft(iim, iim), coff
46    
47        ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
48        real coefilu(iim, jjm), coefilv(iim, jjm)
49        real coefilu2(iim, jjm), coefilv2(iim, jjm)
50    
51        integer modfrstu(jjm), modfrstv(jjm)
52        ! index of the mode from where modes are filtered
53    
54        !-----------------------------------------------------------
55    
56        print *, "Call sequence information: inifilr"
57    
58        CALL inifgn(eignvl)
59    
60        call new_unit(unit)
61        open(unit, file = "eignvl.txt", status = "replace", action = "write")
62        write(unit, fmt = *) EIGNVL
63        close(unit)
64    
65        ! compute eigenvalues and eigenfunctions
66        ! compute the filtering coefficients for scalar lines and
67        ! meridional wind v-lines
68        ! we filter all those latitude lines where coefil < 1
69        ! NO FILTERING AT POLES
70        ! colat0 is to be used when alpha (stretching coefficient)
71        ! is set equal to zero for the regular grid case
72    
73        ! Calcul de colat0
74        forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
75        colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
76        PRINT *, 'colat0 = ', colat0
77    
78        rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
79        coefilu = 0.
80        coefilv = 0.
81        coefilu2 = 0.
82        coefilv2 = 0.
83    
84        ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
85    
86        DO j = 2, jjm / 2 + 1
87           IF (cos(rlatu(j)) / colat0 < 1. &
88                .and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
89    
90           IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
91                .and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
92                jfiltsu = jjm - j + 2
93        END DO
94    
95        DO j = 1, jjm / 2
96           IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
97                jfiltnv = j
98    
99           IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
100                * cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
101        END DO
102    
103        IF (jfiltnu <= 0) jfiltnu = 1
104        IF (jfiltnu > jjm / 2 + 1) THEN
105           PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
106           STOP 1
107        END IF
108    
109        IF (jfiltsu <= 0) jfiltsu = 1
110        IF (jfiltsu > jjm + 1) THEN
111           PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu
112           STOP 1
113        END IF
114    
115        IF (jfiltnv <= 0) jfiltnv = 1
116        IF (jfiltnv > jjm / 2) THEN
117           PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
118           STOP 1
119        END IF
120    
121        IF (jfiltsv <= 0) jfiltsv = 1
122        IF (jfiltsv > jjm) THEN
123           PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv
124           STOP 1
125        END IF
126    
127        PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
128             jfiltsu
129    
130        ! Determination de coefilu, coefilv, n=modfrstu, modfrstv
131    
132        DO j = 1, jjm
133           modfrstu(j) = iim
134           modfrstv(j) = iim
135        END DO
136    
137        DO j = 2, jfiltnu
138           DO k = 2, iim
139              IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
140           end DO
141           if (k == iim + 1) cycle
142           modfrstu(j) = k
143    
144           kf = modfrstu(j)
145           DO k = kf, iim
146              cof = rlamda(k) * cos(rlatu(j))
147              coefilu(k, j) = cof - 1.
148              coefilu2(k, j) = cof**2 - 1.
149           end DO
150        END DO
151    
152        DO j = 1, jfiltnv
153           DO k = 2, iim
154              IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
155           end DO
156           if (k == iim + 1) cycle
157           modfrstv(j) = k
158    
159           kf = modfrstv(j)
160           DO k = kf, iim
161              cof = rlamda(k) * cos(rlatv(j))
162              coefilv(k, j) = cof - 1.
163              coefilv2(k, j) = cof**2 - 1.
164           end DO
165        end DO
166    
167        DO j = jfiltsu, jjm
168           DO k = 2, iim
169              IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
170           end DO
171           if (k == iim + 1) cycle
172           modfrstu(j) = k
173    
174           kf = modfrstu(j)
175           DO k = kf, iim
176              cof = rlamda(k) * cos(rlatu(j))
177              coefilu(k, j) = cof - 1.
178              coefilu2(k, j) = cof**2 - 1.
179           end DO
180        end DO
181    
182        DO j = jfiltsv, jjm
183           DO k = 2, iim
184              IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
185           end DO
186           if (k == iim + 1) cycle
187           modfrstv(j) = k
188    
189           kf = modfrstv(j)
190           DO k = kf, iim
191              cof = rlamda(k) * cos(rlatv(j))
192              coefilv(k, j) = cof - 1.
193              coefilv2(k, j) = cof**2 - 1.
194           end DO
195        END DO
196    
197        IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
198           IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
199           IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
200    
201           PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
202                jfiltsu
203        END IF
204    
205        PRINT *, 'Modes premiers v '
206        PRINT 334, modfrstv
207        PRINT *, 'Modes premiers u '
208        PRINT 334, modfrstu
209    
210        allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
211        allocate(matricevn(iim, iim, jfiltnv))
212        allocate(matricevs(iim, iim, jfiltsv:jjm))
213        allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
214    
215        ! Calcul de la matrice filtre 'matriceu' pour les champs situes
216        ! sur la grille scalaire
217    
218        DO j = 2, jfiltnu
219           DO i = 1, iim
220              IF (i < modfrstu(j)) then
221                 coff = 0.
222              else
223                 coff = coefilu(i, j)
224              end IF
225              eignft(i, :) = eignfnv(:, i) * coff
226           END DO
227           matriceun(:, :, j) = matmul(eignfnv, eignft)
228        END DO
229    
230        DO j = jfiltsu, jjm
231           DO i = 1, iim
232              IF (i < modfrstu(j)) then
233                 coff = 0.
234              else
235                 coff = coefilu(i, j)
236              end IF
237              eignft(i, :) = eignfnv(:, i) * coff
238           END DO
239           matriceus(:, :, j) = matmul(eignfnv, eignft)
240        END DO
241    
242        ! Calcul de la matrice filtre 'matricev' pour les champs situes
243        ! sur la grille de V ou de Z
244    
245        DO j = 1, jfiltnv
246           DO i = 1, iim
247              IF (i < modfrstv(j)) then
248                 coff = 0.
249              else
250                 coff = coefilv(i, j)
251              end IF
252              eignft(i, :) = eignfnu(:, i) * coff
253           END DO
254           matricevn(:, :, j) = matmul(eignfnu, eignft)
255        END DO
256    
257        DO j = jfiltsv, jjm
258           DO i = 1, iim
259              IF (i < modfrstv(j)) then
260                 coff = 0.
261              else
262                 coff = coefilv(i, j)
263              end IF
264              eignft(i, :) = eignfnu(:, i) * coff
265           END DO
266           matricevs(:, :, j) = matmul(eignfnu, eignft)
267        END DO
268    
269        ! Calcul de la matrice filtre 'matrinv' pour les champs situes
270        ! sur la grille scalaire , pour le filtre inverse
271    
272        DO j = 2, jfiltnu
273           DO i = 1, iim
274              IF (i < modfrstu(j)) then
275                 coff = 0.
276              else
277                 coff = coefilu(i, j) / (1. + coefilu(i, j))
278              end IF
279              eignft(i, :) = eignfnv(:, i) * coff
280           END DO
281           matrinvn(:, :, j) = matmul(eignfnv, eignft)
282        END DO
283    
284        DO j = jfiltsu, jjm
285           DO i = 1, iim
286              IF (i < modfrstu(j)) then
287                 coff = 0.
288              else
289                 coff = coefilu(i, j) / (1. + coefilu(i, j))
290              end IF
291              eignft(i, :) = eignfnv(:, i) * coff
292           END DO
293           matrinvs(:, :, j) = matmul(eignfnv, eignft)
294        END DO
295    
296    334 FORMAT (1X, 24I3)
297    
298      END SUBROUTINE inifilr
299    
300    end module inifilr_m
Legend:
Removed from v.5  
changed lines
  Added in v.143
  ViewVC Help
Powered by ViewVC 1.1.21