Ignore:
Timestamp:
Jul 23, 2024, 10:21:18 PM (2 months ago)
Author:
abarral
Message:

Turn coefils.h into lmdz_coefils.f90
Put filtreg.F90 inside lmdz_filtreg.F90
Turn mod_filtreg_p.F90 into lmdz_filtreg_p.F90
Delete obsolete parafilt.h*
(lint) remove spaces between routine name and args

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/filtrez/lmdz_filtreg.F90

    r5105 r5106  
    1 
    21! $Id$
    32
    4 MODULE filtreg_mod
    5 
    6   REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
    7   REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
     3MODULE lmdz_filtreg
     4  IMPLICIT NONE; PRIVATE
     5  PUBLIC matriceun, matriceus, matricevn, matricevs, matrinvn, matrinvs, &
     6          inifilr, filtreg
     7
     8  REAL, DIMENSION(:, :, :), ALLOCATABLE :: matriceun, matriceus, matricevn
     9  REAL, DIMENSION(:, :, :), ALLOCATABLE :: matricevs, matrinvn, matrinvs
    810
    911CONTAINS
     12
     13  SUBROUTINE filtreg(champ, nlat, nbniv, ifiltre, iaire, &
     14          griscal, iter)
     15    USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu
     16
     17    !=======================================================================
     18    !
     19    !   Auteur: P. Le Van        07/10/97
     20    !   ------
     21    !
     22    !   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
     23    !                 pour l'operateur  Filtre    .
     24    !   ------
     25    !
     26    !   Arguments:
     27    !   ----------
     28    !
     29    !  nblat                 nombre de latitudes a filtrer
     30    !  nbniv                 nombre de niveaux verticaux a filtrer
     31    !  champ(iip1,nblat,nbniv)  en entree : champ a filtrer
     32    !                        en sortie : champ filtre
     33    !  ifiltre               +1  Transformee directe
     34    !                        -1  Transformee inverse
     35    !                        +2  Filtre directe
     36    !                        -2  Filtre inverse
     37    !
     38    !  iaire                 1   si champ intensif
     39    !                        2   si champ extensif (pondere par les aires)
     40    !
     41    !  iter                  1   filtre simple
     42    !
     43    !=======================================================================
     44    !
     45    !
     46    !                  Variable Intensive
     47    !            ifiltre = 1     filtre directe
     48    !            ifiltre =-1     filtre inverse
     49    !
     50    !                  Variable Extensive
     51    !            ifiltre = 2     filtre directe
     52    !            ifiltre =-2     filtre inverse
     53    !
     54    !
     55    INCLUDE "dimensions.h"
     56    INCLUDE "paramet.h"
     57
     58    INTEGER :: nlat, nbniv, ifiltre, iter
     59    INTEGER :: i, j, l, k
     60    INTEGER :: iim2, immjm
     61    INTEGER :: jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil
     62
     63    REAL :: champ(iip1, nlat, nbniv)
     64
     65    REAL :: eignq(iim, nlat, nbniv), sdd1(iim), sdd2(iim)
     66    LOGICAL :: griscal
     67    INTEGER :: hemisph, iaire
     68
     69    LOGICAL, SAVE :: first = .TRUE.
     70
     71    REAL, SAVE :: sdd12(iim, 4)
     72
     73    INTEGER, PARAMETER :: type_sddu = 1
     74    INTEGER, PARAMETER :: type_sddv = 2
     75    INTEGER, PARAMETER :: type_unsddu = 3
     76    INTEGER, PARAMETER :: type_unsddv = 4
     77
     78    INTEGER :: sdd1_type, sdd2_type
     79
     80    if (iim == 1) return ! no filtre in 2D y-z
     81
     82    IF (first) THEN
     83      sdd12(1:iim, type_sddu) = sddu(1:iim)
     84      sdd12(1:iim, type_sddv) = sddv(1:iim)
     85      sdd12(1:iim, type_unsddu) = unsddu(1:iim)
     86      sdd12(1:iim, type_unsddv) = unsddv(1:iim)
     87
     88      first = .FALSE.
     89    ENDIF
     90
     91    IF(ifiltre==1.or.ifiltre==-1) &
     92            stop 'Pas de transformee simple dans cette version'
     93
     94    IF(iter== 2)  THEN
     95      PRINT *, ' Pas d iteration du filtre dans cette version !'&
     96              &, ' Utiliser old_filtreg et repasser !'
     97      STOP
     98    ENDIF
     99
     100    IF(ifiltre== -2 .AND..NOT.griscal)     THEN
     101      PRINT *, ' Cette routine ne calcule le filtre inverse que ' &
     102              , ' sur la grille des scalaires !'
     103      STOP
     104    ENDIF
     105
     106    IF(ifiltre/=2 .AND.ifiltre/= - 2)  THEN
     107      PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
     108              , ' corriger et repasser !'
     109      STOP
     110    ENDIF
     111
     112    iim2 = iim * iim
     113    immjm = iim * jjm
     114
     115    IF(griscal)   THEN
     116      IF(nlat /= jjp1)  THEN
     117        PRINT  1111
     118        STOP
     119      ELSE
     120
     121        IF(iaire==1)  THEN
     122          sdd1_type = type_sddv
     123          sdd2_type = type_unsddv
     124        ELSE
     125          sdd1_type = type_unsddv
     126          sdd2_type = type_sddv
     127        ENDIF
     128
     129        ! IF( iaire.EQ.1 )  THEN
     130        !    CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
     131        !    CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
     132        ! ELSE
     133        !    CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
     134        !    CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
     135        ! END IF
     136
     137        jdfil1 = 2
     138        jffil1 = jfiltnu
     139        jdfil2 = jfiltsu
     140        jffil2 = jjm
     141      END IF
     142    ELSE
     143      IF(nlat/=jjm)  THEN
     144        PRINT  2222
     145        STOP
     146      ELSE
     147
     148        IF(iaire==1)  THEN
     149          sdd1_type = type_sddu
     150          sdd2_type = type_unsddu
     151        ELSE
     152          sdd1_type = type_unsddu
     153          sdd2_type = type_sddu
     154        ENDIF
     155
     156        ! IF( iaire.EQ.1 )  THEN
     157        !    CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
     158        !    CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
     159        ! ELSE
     160        !    CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
     161        !    CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
     162        ! END IF
     163
     164        jdfil1 = 1
     165        jffil1 = jfiltnv
     166        jdfil2 = jfiltsv
     167        jffil2 = jjm
     168      END IF
     169    END IF
     170
     171    DO hemisph = 1, 2
     172
     173      IF (hemisph==1)  THEN
     174        jdfil = jdfil1
     175        jffil = jffil1
     176      ELSE
     177        jdfil = jdfil2
     178        jffil = jffil2
     179      END IF
     180
     181      DO l = 1, nbniv
     182        DO j = jdfil, jffil
     183          DO i = 1, iim
     184            champ(i, j, l) = champ(i, j, l) * sdd12(i, sdd1_type) ! sdd1(i)
     185          END DO
     186        END DO
     187      END DO
     188
     189      IF(hemisph == 1)      THEN
     190
     191        IF(ifiltre == -2)   THEN
     192
     193          DO j = jdfil, jffil
     194#ifdef BLAS
     195              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
     196                    matrinvn(1,1,j), &
     197                    iim, champ(1,j,1), iip1*nlat, 0.0, &
     198                    eignq(1,j-jdfil+1,1), iim*nlat)
     199#else
     200            eignq(:, j - jdfil + 1, :) &
     201                    = matmul(matrinvn(:, :, j), champ(:iim, j, :))
     202#endif
     203          END DO
     204
     205        ELSE IF (griscal)     THEN
     206
     207          DO j = jdfil, jffil
     208#ifdef BLAS
     209              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
     210                    matriceun(1,1,j), &
     211                    iim, champ(1,j,1), iip1*nlat, 0.0, &
     212                    eignq(1,j-jdfil+1,1), iim*nlat)
     213#else
     214            eignq(:, j - jdfil + 1, :) &
     215                    = matmul(matriceun(:, :, j), champ(:iim, j, :))
     216#endif
     217          END DO
     218
     219        ELSE
     220
     221          DO j = jdfil, jffil
     222#ifdef BLAS
     223              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
     224                    matricevn(1,1,j), &
     225                    iim, champ(1,j,1), iip1*nlat, 0.0, &
     226                    eignq(1,j-jdfil+1,1), iim*nlat)
     227#else
     228            eignq(:, j - jdfil + 1, :) &
     229                    = matmul(matricevn(:, :, j), champ(:iim, j, :))
     230#endif
     231          END DO
     232
     233        ENDIF
     234
     235      ELSE
     236
     237        IF(ifiltre == -2)   THEN
     238
     239          DO j = jdfil, jffil
     240#ifdef BLAS
     241              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
     242                    matrinvs(1,1,j-jfiltsu+1), &
     243                    iim, champ(1,j,1), iip1*nlat, 0.0, &
     244                    eignq(1,j-jdfil+1,1), iim*nlat)
     245#else
     246            eignq(:, j - jdfil + 1, :) &
     247                    = matmul(matrinvs(:, :, j - jfiltsu + 1), &
     248                    champ(:iim, j, :))
     249#endif
     250          END DO
     251
     252        ELSE IF (griscal)     THEN
     253
     254          DO j = jdfil, jffil
     255#ifdef BLAS
     256              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
     257                    matriceus(1,1,j-jfiltsu+1), &
     258                    iim, champ(1,j,1), iip1*nlat, 0.0, &
     259                    eignq(1,j-jdfil+1,1), iim*nlat)
     260#else
     261            eignq(:, j - jdfil + 1, :) &
     262                    = matmul(matriceus(:, :, j - jfiltsu + 1), &
     263                    champ(:iim, j, :))
     264#endif
     265          END DO
     266
     267        ELSE
     268
     269          DO j = jdfil, jffil
     270#ifdef BLAS
     271              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
     272                    matricevs(1,1,j-jfiltsv+1), &
     273                    iim, champ(1,j,1), iip1*nlat, 0.0, &
     274                    eignq(1,j-jdfil+1,1), iim*nlat)
     275#else
     276            eignq(:, j - jdfil + 1, :) &
     277                    = matmul(matricevs(:, :, j - jfiltsv + 1), &
     278                    champ(:iim, j, :))
     279#endif
     280          END DO
     281
     282        ENDIF
     283
     284      ENDIF
     285
     286      IF(ifiltre== 2)  THEN
     287
     288        DO l = 1, nbniv
     289          DO j = jdfil, jffil
     290            DO i = 1, iim
     291              champ(i, j, l) = &
     292                      (champ(i, j, l) + eignq(i, j - jdfil + 1, l)) &
     293                              * sdd12(i, sdd2_type) ! sdd2(i)
     294            END DO
     295          END DO
     296        END DO
     297
     298      ELSE
     299
     300        DO l = 1, nbniv
     301          DO j = jdfil, jffil
     302            DO i = 1, iim
     303              champ(i, j, l) = &
     304                      (champ(i, j, l) - eignq(i, j - jdfil + 1, l)) &
     305                              * sdd12(i, sdd2_type) ! sdd2(i)
     306            END DO
     307          END DO
     308        END DO
     309
     310      ENDIF
     311
     312      DO l = 1, nbniv
     313        DO j = jdfil, jffil
     314          champ(iip1, j, l) = champ(1, j, l)
     315        END DO
     316      END DO
     317
     318    ENDDO
     319
     320    1111   FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau  CHAMP a&
     321            &     filtrer, sur la grille des scalaires'/)
     322    2222   FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau CHAMP a fi&
     323            &     ltrer, sur la grille de V ou de Z'/)
     324    RETURN
     325  END SUBROUTINE filtreg
    10326
    11327  SUBROUTINE inifilr
    … …  
    14330  USE mod_filtre_fft_loc, ONLY: Init_filtre_fft_loc=>Init_filtre_fft    !
    15331#endif
    16   USE serre_mod, ONLY: alphax
    17   USE logic_mod, ONLY: fxyhypb, ysinus
    18   USE comconst_mod, ONLY: maxlatfilter
    19 
    20     !    ... H. Upadhyaya, O.Sharma   ...
    21 
    22     IMPLICIT NONE
     332    USE serre_mod, ONLY: alphax
     333    USE logic_mod, ONLY: fxyhypb, ysinus
     334    USE comconst_mod, ONLY: maxlatfilter
     335    USE lmdz_coefils, ONLY: modfrstv, modfrstu, jfiltnu, jfiltnv, coefilu, coefilv, &
     336            coefilu2, coefilv2, eignfnv, eignfnu, jfiltsu, jfiltsv
     337
     338    !    ... H. Upfiltreg_modadhyaya, O.Sharma   ...
    23339
    24340    !     version 3 .....
    … …  
    28344    include "dimensions.h"
    29345    include "paramet.h"
    30     !  -------------------------------------------------------------------
    31346    include "comgeom.h"
    32     include "coefils.h"
    33 
    34     REAL  dlonu(iim),dlatu(jjm)
    35     REAL  rlamda( iim ),  eignvl( iim )
    36 
    37     REAL    lamdamax,pi,cof
    38     INTEGER i,j,modemax,imx,k,kf,ii
    39     REAL dymin,dxmin,colat0
    40     REAL eignft(iim,iim), coff
     347
     348    REAL  dlonu(iim), dlatu(jjm)
     349    REAL  rlamda(iim), eignvl(iim)
     350
     351    REAL    lamdamax, pi, cof
     352    INTEGER i, j, modemax, imx, k, kf, ii
     353    REAL dymin, dxmin, colat0
     354    REAL eignft(iim, iim), coff
    41355
    42356    LOGICAL, SAVE :: first_call_inifilr = .TRUE.
    … …  
    44358    INTEGER   ISMIN
    45359    EXTERNAL  ISMIN
    46     INTEGER iymin 
     360    INTEGER iymin
    47361    INTEGER ixmineq
    48362
    … …  
    65379    !-----------------------------------------------------------
    66380
    67     if ( iim == 1 ) return ! No filtre in 2D y-z
    68 
    69     pi       = 2. * ASIN( 1. )
    70 
    71     DO i = 1,iim
    72        dlonu(i) = xprimu( i )
     381    if (iim == 1) return ! No filtre in 2D y-z
     382
     383    pi = 2. * ASIN(1.)
     384
     385    DO i = 1, iim
     386      dlonu(i) = xprimu(i)
    73387    ENDDO
    74388
    75389    CALL inifgn(eignvl)
    76390
    77     PRINT *,'inifilr: EIGNVL '
    78     PRINT 250,eignvl
    79 250 FORMAT( 1x,5e14.6)
     391    PRINT *, 'inifilr: EIGNVL '
     392    PRINT 250, eignvl
     393    250 FORMAT(1x, 5e14.6)
    80394
    81395    ! compute eigenvalues and eigenfunctions
    … …  
    96410    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
    97411
    98 
    99     DO j = 1,jjm
    100        dlatu( j ) = rlatu( j ) - rlatu( j+1 )
    101     ENDDO
    102 
    103     dxmin   =  dlonu(1)
    104     DO  i  = 2, iim
    105        dxmin = MIN( dxmin,dlonu(i) )
    106     ENDDO
    107     dymin  = dlatu(1)
    108     DO j  = 2, jjm
    109        dymin = MIN( dymin,dlatu(j) )
     412    DO j = 1, jjm
     413      dlatu(j) = rlatu(j) - rlatu(j + 1)
     414    ENDDO
     415
     416    dxmin = dlonu(1)
     417    DO  i = 2, iim
     418      dxmin = MIN(dxmin, dlonu(i))
     419    ENDDO
     420    dymin = dlatu(1)
     421    DO j = 2, jjm
     422      dymin = MIN(dymin, dlatu(j))
    110423    ENDDO
    111424
    … …  
    118431
    119432    ! if maxlatfilter >0, prescribe the colat0 value from the .def files
    120    
     433
    121434    IF (maxlatfilter < 0.) THEN
    122435
    123     colat0  =  MIN( 0.5, dymin/dxmin )
    124     ! colat0  =  1.
    125 
    126     IF( .NOT.fxyhypb.AND.ysinus )  THEN
    127        colat0 = 0.6
    128        !         ...... a revoir  pour  ysinus !   .......
    129        alphax = 0.
    130     ENDIF
     436      colat0 = MIN(0.5, dymin / dxmin)
     437      ! colat0  =  1.
     438
     439      IF(.NOT.fxyhypb.AND.ysinus)  THEN
     440        colat0 = 0.6
     441        !         ...... a revoir  pour  ysinus !   .......
     442        alphax = 0.
     443      ENDIF
    131444
    132445    ELSE
    133446
    134     colat0=(90.0-maxlatfilter)/180.0*pi       
    135 
    136     ENDIF
    137 
    138     PRINT 50, colat0,alphax
    139 50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
    140 
    141     IF(alphax==1. )  THEN
    142        PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
    143        STOP
    144     ENDIF
    145 
    146     lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
     447      colat0 = (90.0 - maxlatfilter) / 180.0 * pi
     448
     449    ENDIF
     450
     451    PRINT 50, colat0, alphax
     452    50  FORMAT(/15x, ' Inifilr colat0 alphax ', 2e16.7)
     453
     454    IF(alphax==1.)  THEN
     455      PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '
     456      STOP
     457    ENDIF
     458
     459    lamdamax = iim / (pi * colat0 * (1. - alphax))
    147460
    148461    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
    149462
    150     DO i = 2,iim
    151        rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
    152     ENDDO
    153 
    154     DO j = 1,jjm
    155        DO i = 1,iim
    156           coefilu( i,j ) = 0.0
    157           coefilv( i,j ) = 0.0
    158           coefilu2( i,j ) = 0.0
    159           coefilv2( i,j ) = 0.0
    160        ENDDO
     463    DO i = 2, iim
     464      rlamda(i) = lamdamax / SQRT(ABS(eignvl(i)))
     465    ENDDO
     466
     467    DO j = 1, jjm
     468      DO i = 1, iim
     469        coefilu(i, j) = 0.0
     470        coefilv(i, j) = 0.0
     471        coefilu2(i, j) = 0.0
     472        coefilv2(i, j) = 0.0
     473      ENDDO
    161474    ENDDO
    162475
    … …  
    166479    modemax = iim
    167480
    168 !!!!    imx = modemax - 4 * (modemax/iim)
    169 
    170     imx  = iim
    171 
    172     PRINT *,'inifilr: TRUNCATION AT ',imx
    173 
    174 ! Ehouarn: set up some defaults
    175     jfiltnu=2 ! avoid north pole
    176     jfiltsu=jjm ! avoid south pole (which is at jjm+1)
    177     jfiltnv=1 ! NB: no poles on the V grid
    178     jfiltsv=jjm
    179 
    180     DO j = 2, jjm/2+1
    181        cof = COS( rlatu(j) )/ colat0
    182        IF ( cof < 1. ) THEN
    183           IF( rlamda(imx) * COS(rlatu(j) )<1. ) THEN
    184             jfiltnu= j
    185           ENDIF
    186        ENDIF
    187 
    188        cof = COS( rlatu(jjp1-j+1) )/ colat0
    189        IF ( cof < 1. ) THEN
    190           IF( rlamda(imx) * COS(rlatu(jjp1-j+1) )<1. ) THEN
    191                jfiltsu= jjp1-j+1
    192           ENDIF
    193        ENDIF
    194     ENDDO
    195 
    196     DO j = 1, jjm/2
    197        cof = COS( rlatv(j) )/ colat0
    198        IF ( cof < 1. ) THEN
    199           IF( rlamda(imx) * COS(rlatv(j) )<1. ) THEN
    200             jfiltnv= j
    201           ENDIF
    202        ENDIF
    203 
    204        cof = COS( rlatv(jjm-j+1) )/ colat0
    205        IF ( cof < 1. ) THEN
    206           IF( rlamda(imx) * COS(rlatv(jjm-j+1) )<1. ) THEN
    207                jfiltsv= jjm-j+1
    208           ENDIF
    209        ENDIF
    210     ENDDO
    211 
    212     IF( jfiltnu> jjm/2 +1 )  THEN
    213        PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
    214        STOP
    215     ENDIF
    216 
    217     IF( jfiltsu>  jjm  +1 )  THEN
    218        PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
    219        STOP
    220     ENDIF
    221 
    222     IF( jfiltnv> jjm/2    )  THEN
    223        PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
    224        STOP
    225     ENDIF
    226 
    227     IF( jfiltsv>     jjm  )  THEN
    228        PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
    229        STOP
    230     ENDIF
    231 
    232     PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
    233          jfiltnv,jfiltsv,jfiltnu,jfiltsu
     481    !!!!    imx = modemax - 4 * (modemax/iim)
     482
     483    imx = iim
     484
     485    PRINT *, 'inifilr: TRUNCATION AT ', imx
     486
     487    ! Ehouarn: set up some defaults
     488    jfiltnu = 2 ! avoid north pole
     489    jfiltsu = jjm ! avoid south pole (which is at jjm+1)
     490    jfiltnv = 1 ! NB: no poles on the V grid
     491    jfiltsv = jjm
     492
     493    DO j = 2, jjm / 2 + 1
     494      cof = COS(rlatu(j)) / colat0
     495      IF (cof < 1.) THEN
     496        IF(rlamda(imx) * COS(rlatu(j))<1.) THEN
     497          jfiltnu = j
     498        ENDIF
     499      ENDIF
     500
     501      cof = COS(rlatu(jjp1 - j + 1)) / colat0
     502      IF (cof < 1.) THEN
     503        IF(rlamda(imx) * COS(rlatu(jjp1 - j + 1))<1.) THEN
     504          jfiltsu = jjp1 - j + 1
     505        ENDIF
     506      ENDIF
     507    ENDDO
     508
     509    DO j = 1, jjm / 2
     510      cof = COS(rlatv(j)) / colat0
     511      IF (cof < 1.) THEN
     512        IF(rlamda(imx) * COS(rlatv(j))<1.) THEN
     513          jfiltnv = j
     514        ENDIF
     515      ENDIF
     516
     517      cof = COS(rlatv(jjm - j + 1)) / colat0
     518      IF (cof < 1.) THEN
     519        IF(rlamda(imx) * COS(rlatv(jjm - j + 1))<1.) THEN
     520          jfiltsv = jjm - j + 1
     521        ENDIF
     522      ENDIF
     523    ENDDO
     524
     525    IF(jfiltnu> jjm / 2 + 1)  THEN
     526      PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
     527      STOP
     528    ENDIF
     529
     530    IF(jfiltsu>  jjm + 1)  THEN
     531      PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
     532      STOP
     533    ENDIF
     534
     535    IF(jfiltnv> jjm / 2)  THEN
     536      PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
     537      STOP
     538    ENDIF
     539
     540    IF(jfiltsv>     jjm)  THEN
     541      PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
     542      STOP
     543    ENDIF
     544
     545    PRINT *, 'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ', &
     546            jfiltnv, jfiltsv, jfiltnu, jfiltsu
    234547
    235548    IF(first_call_inifilr) THEN
    236        ALLOCATE(matriceun(iim,iim,jfiltnu))
    237        ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
    238        ALLOCATE(matricevn(iim,iim,jfiltnv))
    239        ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
    240        ALLOCATE( matrinvn(iim,iim,jfiltnu))
    241        ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
    242        first_call_inifilr = .FALSE.
     549      ALLOCATE(matriceun(iim, iim, jfiltnu))
     550      ALLOCATE(matriceus(iim, iim, jjm - jfiltsu + 1))
     551      ALLOCATE(matricevn(iim, iim, jfiltnv))
     552      ALLOCATE(matricevs(iim, iim, jjm - jfiltsv + 1))
     553      ALLOCATE(matrinvn(iim, iim, jfiltnu))
     554      ALLOCATE(matrinvs(iim, iim, jjm - jfiltsu + 1))
     555      first_call_inifilr = .FALSE.
    243556    ENDIF
    244557
    … …  
    246559    !................................................................
    247560
    248 
    249     DO j = 1,jjm
    250     !default initialization: all modes are retained (i.e. no filtering)
    251        modfrstu( j ) = iim
    252        modfrstv( j ) = iim
    253     ENDDO
    254 
    255     DO j = 2,jfiltnu
    256        DO k = 2,modemax
    257           cof = rlamda(k) * COS( rlatu(j) )
    258           IF ( cof < 1. ) GOTO 82
    259        ENDDO
    260        GOTO 84
    261 82     modfrstu( j ) = k
    262 
    263        kf = modfrstu( j )
    264        DO k = kf , modemax
    265           cof = rlamda(k) * COS( rlatu(j) )
    266           coefilu(k,j) = cof - 1.
    267           coefilu2(k,j) = cof*cof - 1.
    268        ENDDO
    269 84     CONTINUE
    270     ENDDO
    271 
    272 
    273     DO j = 1,jfiltnv
    274 
    275        DO k = 2,modemax
    276           cof = rlamda(k) * COS( rlatv(j) )
    277           IF ( cof < 1. ) GOTO 87
    278        ENDDO
    279        GOTO 89
    280 87     modfrstv( j ) = k
    281 
    282        kf = modfrstv( j )
    283        DO k = kf , modemax
    284           cof = rlamda(k) * COS( rlatv(j) )
    285           coefilv(k,j) = cof - 1.
    286           coefilv2(k,j) = cof*cof - 1.
    287        ENDDO
    288 89     CONTINUE
    289     ENDDO
    290 
    291     DO j = jfiltsu,jjm
    292        DO k = 2,modemax
    293           cof = rlamda(k) * COS( rlatu(j) )
    294           IF ( cof < 1. ) GOTO 92
    295        ENDDO
    296        GOTO 94
    297 92     modfrstu( j ) = k
    298 
    299        kf = modfrstu( j )
    300        DO k = kf , modemax
    301           cof = rlamda(k) * COS( rlatu(j) )
    302           coefilu(k,j) = cof - 1.
    303           coefilu2(k,j) = cof*cof - 1.
    304        ENDDO
    305 94     CONTINUE
    306     ENDDO
    307 
    308     DO j = jfiltsv,jjm
    309        DO k = 2,modemax
    310           cof = rlamda(k) * COS( rlatv(j) )
    311           IF ( cof < 1. ) GOTO 97
    312        ENDDO
    313        GOTO 99
    314 97     modfrstv( j ) = k
    315 
    316        kf = modfrstv( j )
    317        DO k = kf , modemax
    318           cof = rlamda(k) * COS( rlatv(j) )
    319           coefilv(k,j) = cof - 1.
    320           coefilv2(k,j) = cof*cof - 1.
    321        ENDDO
    322 99     CONTINUE
    323     ENDDO
    324 
    325     IF(jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2)THEN
    326 ! Ehouarn: and what are these for??? Trying to handle a limit case
    327 !          where filters extend to and meet at the equator?
    328        IF(jfiltnv==jfiltsv)jfiltsv=1+jfiltnv
    329        IF(jfiltnu==jfiltsu)jfiltsu=1+jfiltnu
    330 
    331        PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
    332             jfiltnv,jfiltsv,jfiltnu,jfiltsu
    333     ENDIF
    334 
    335     PRINT *,'   Modes premiers  v  '
    336     PRINT 334,modfrstv
    337     PRINT *,'   Modes premiers  u  '
    338     PRINT 334,modfrstu
     561    DO j = 1, jjm
     562      !default initialization: all modes are retained (i.e. no filtering)
     563      modfrstu(j) = iim
     564      modfrstv(j) = iim
     565    ENDDO
     566
     567    DO j = 2, jfiltnu
     568      DO k = 2, modemax
     569        cof = rlamda(k) * COS(rlatu(j))
     570        IF (cof < 1.) GOTO 82
     571      ENDDO
     572      GOTO 84
     573      82     modfrstu(j) = k
     574
     575      kf = modfrstu(j)
     576      DO k = kf, modemax
     577        cof = rlamda(k) * COS(rlatu(j))
     578        coefilu(k, j) = cof - 1.
     579        coefilu2(k, j) = cof * cof - 1.
     580      ENDDO
     581      84     CONTINUE
     582    ENDDO
     583
     584    DO j = 1, jfiltnv
     585
     586      DO k = 2, modemax
     587        cof = rlamda(k) * COS(rlatv(j))
     588        IF (cof < 1.) GOTO 87
     589      ENDDO
     590      GOTO 89
     591      87     modfrstv(j) = k
     592
     593      kf = modfrstv(j)
     594      DO k = kf, modemax
     595        cof = rlamda(k) * COS(rlatv(j))
     596        coefilv(k, j) = cof - 1.
     597        coefilv2(k, j) = cof * cof - 1.
     598      ENDDO
     599      89     CONTINUE
     600    ENDDO
     601
     602    DO j = jfiltsu, jjm
     603      DO k = 2, modemax
     604        cof = rlamda(k) * COS(rlatu(j))
     605        IF (cof < 1.) GOTO 92
     606      ENDDO
     607      GOTO 94
     608      92     modfrstu(j) = k
     609
     610      kf = modfrstu(j)
     611      DO k = kf, modemax
     612        cof = rlamda(k) * COS(rlatu(j))
     613        coefilu(k, j) = cof - 1.
     614        coefilu2(k, j) = cof * cof - 1.
     615      ENDDO
     616      94     CONTINUE
     617    ENDDO
     618
     619    DO j = jfiltsv, jjm
     620      DO k = 2, modemax
     621        cof = rlamda(k) * COS(rlatv(j))
     622        IF (cof < 1.) GOTO 97
     623      ENDDO
     624      GOTO 99
     625      97     modfrstv(j) = k
     626
     627      kf = modfrstv(j)
     628      DO k = kf, modemax
     629        cof = rlamda(k) * COS(rlatv(j))
     630        coefilv(k, j) = cof - 1.
     631        coefilv2(k, j) = cof * cof - 1.
     632      ENDDO
     633      99     CONTINUE
     634    ENDDO
     635
     636    IF(jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2)THEN
     637      ! Ehouarn: and what are these for??? Trying to handle a limit case
     638      !          where filters extend to and meet at the equator?
     639      IF(jfiltnv==jfiltsv)jfiltsv = 1 + jfiltnv
     640      IF(jfiltnu==jfiltsu)jfiltsu = 1 + jfiltnu
     641
     642      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', &
     643              jfiltnv, jfiltsv, jfiltnu, jfiltsu
     644    ENDIF
     645
     646    PRINT *, '   Modes premiers  v  '
     647    PRINT 334, modfrstv
     648    PRINT *, '   Modes premiers  u  '
     649    PRINT 334, modfrstu
    339650
    340651    !   ...................................................................
    … …  
    346657    DO j = 2, jfiltnu
    347658
    348        DO i=1,iim
    349           coff = coefilu(i,j)
    350           IF( i<modfrstu(j) ) coff = 0.
    351           DO k=1,iim
    352              eignft(i,k) = eignfnv(k,i) * coff
    353           ENDDO
    354        ENDDO ! of DO i=1,iim
     659      DO i = 1, iim
     660        coff = coefilu(i, j)
     661        IF(i<modfrstu(j)) coff = 0.
     662        DO k = 1, iim
     663          eignft(i, k) = eignfnv(k, i) * coff
     664        ENDDO
     665      ENDDO ! of DO i=1,iim
    355666
    356667#ifdef BLAS
    … …  
    358669            eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
    359670#else
    360        DO k = 1, iim
    361           DO i = 1, iim
    362              matriceun(i,k,j) = 0.0
    363              DO ii = 1, iim
    364                 matriceun(i,k,j) = matriceun(i,k,j) &
    365                      + eignfnv(i,ii)*eignft(ii,k)
    366              ENDDO
     671      DO k = 1, iim
     672        DO i = 1, iim
     673          matriceun(i, k, j) = 0.0
     674          DO ii = 1, iim
     675            matriceun(i, k, j) = matriceun(i, k, j) &
     676                    + eignfnv(i, ii) * eignft(ii, k)
    367677          ENDDO
    368        ENDDO ! of DO k = 1, iim
     678        ENDDO
     679      ENDDO ! of DO k = 1, iim
    369680#endif
    370681
    … …  
    373684    DO j = jfiltsu, jjm
    374685
    375        DO i=1,iim
    376           coff = coefilu(i,j)
    377           IF( i<modfrstu(j) ) coff = 0.
    378           DO k=1,iim
    379              eignft(i,k) = eignfnv(k,i) * coff
    380           ENDDO
    381        ENDDO ! of DO i=1,iim
     686      DO i = 1, iim
     687        coff = coefilu(i, j)
     688        IF(i<modfrstu(j)) coff = 0.
     689        DO k = 1, iim
     690          eignft(i, k) = eignfnv(k, i) * coff
     691        ENDDO
     692      ENDDO ! of DO i=1,iim
    382693#ifdef BLAS
    383694       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    … …  
    385696            matriceus(1,1,j-jfiltsu+1), iim)
    386697#else
    387        DO k = 1, iim
    388           DO i = 1, iim
    389              matriceus(i,k,j-jfiltsu+1) = 0.0
    390              DO ii = 1, iim
    391                 matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
    392                      + eignfnv(i,ii)*eignft(ii,k)
    393              ENDDO
     698      DO k = 1, iim
     699        DO i = 1, iim
     700          matriceus(i, k, j - jfiltsu + 1) = 0.0
     701          DO ii = 1, iim
     702            matriceus(i, k, j - jfiltsu + 1) = matriceus(i, k, j - jfiltsu + 1) &
     703                    + eignfnv(i, ii) * eignft(ii, k)
    394704          ENDDO
    395        ENDDO ! of DO k = 1, iim
     705        ENDDO
     706      ENDDO ! of DO k = 1, iim
    396707#endif
    397708
    … …  
    406717    DO j = 1, jfiltnv
    407718
    408        DO i = 1, iim
    409           coff = coefilv(i,j)
    410           IF( i<modfrstv(j) ) coff = 0.
    411           DO k = 1, iim
    412              eignft(i,k) = eignfnu(k,i) * coff
    413           ENDDO
    414        ENDDO
     719      DO i = 1, iim
     720        coff = coefilv(i, j)
     721        IF(i<modfrstv(j)) coff = 0.
     722        DO k = 1, iim
     723          eignft(i, k) = eignfnu(k, i) * coff
     724        ENDDO
     725      ENDDO
    415726
    416727#ifdef BLAS
    … …  
    418729            eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
    419730#else
    420        DO k = 1, iim
    421           DO i = 1, iim
    422              matricevn(i,k,j) = 0.0
    423              DO ii = 1, iim
    424                 matricevn(i,k,j) = matricevn(i,k,j) &
    425                      + eignfnu(i,ii)*eignft(ii,k)
    426              ENDDO
     731      DO k = 1, iim
     732        DO i = 1, iim
     733          matricevn(i, k, j) = 0.0
     734          DO ii = 1, iim
     735            matricevn(i, k, j) = matricevn(i, k, j) &
     736                    + eignfnu(i, ii) * eignft(ii, k)
    427737          ENDDO
    428        ENDDO
     738        ENDDO
     739      ENDDO
    429740#endif
    430741
    … …  
    433744    DO j = jfiltsv, jjm
    434745
    435        DO i = 1, iim
    436           coff = coefilv(i,j)
    437           IF( i<modfrstv(j) ) coff = 0.
    438           DO k = 1, iim
    439              eignft(i,k) = eignfnu(k,i) * coff
    440           ENDDO
    441        ENDDO
     746      DO i = 1, iim
     747        coff = coefilv(i, j)
     748        IF(i<modfrstv(j)) coff = 0.
     749        DO k = 1, iim
     750          eignft(i, k) = eignfnu(k, i) * coff
     751        ENDDO
     752      ENDDO
    442753
    443754#ifdef BLAS
    … …  
    446757            matricevs(1,1,j-jfiltsv+1), iim)
    447758#else
    448        DO k = 1, iim
    449           DO i = 1, iim
    450              matricevs(i,k,j-jfiltsv+1) = 0.0
    451              DO ii = 1, iim
    452                 matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
    453                      + eignfnu(i,ii)*eignft(ii,k)
    454              ENDDO
     759      DO k = 1, iim
     760        DO i = 1, iim
     761          matricevs(i, k, j - jfiltsv + 1) = 0.0
     762          DO ii = 1, iim
     763            matricevs(i, k, j - jfiltsv + 1) = matricevs(i, k, j - jfiltsv + 1) &
     764                    + eignfnu(i, ii) * eignft(ii, k)
    455765          ENDDO
    456        ENDDO
     766        ENDDO
     767      ENDDO
    457768#endif
    458769
    … …  
    467778    DO j = 2, jfiltnu
    468779
    469        DO i = 1,iim
    470           coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
    471           IF( i<modfrstu(j) ) coff = 0.
    472           DO k=1,iim
    473              eignft(i,k) = eignfnv(k,i) * coff
    474           ENDDO
    475        ENDDO
     780      DO i = 1, iim
     781        coff = coefilu(i, j) / (1. + coefilu(i, j))
     782        IF(i<modfrstu(j)) coff = 0.
     783        DO k = 1, iim
     784          eignft(i, k) = eignfnv(k, i) * coff
     785        ENDDO
     786      ENDDO
    476787
    477788#ifdef BLAS
    … …  
    479790            eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
    480791#else
    481        DO k = 1, iim
    482           DO i = 1, iim
    483              matrinvn(i,k,j) = 0.0
    484              DO ii = 1, iim
    485                 matrinvn(i,k,j) = matrinvn(i,k,j) &
    486                      + eignfnv(i,ii)*eignft(ii,k)
    487              ENDDO
     792      DO k = 1, iim
     793        DO i = 1, iim
     794          matrinvn(i, k, j) = 0.0
     795          DO ii = 1, iim
     796            matrinvn(i, k, j) = matrinvn(i, k, j) &
     797                    + eignfnv(i, ii) * eignft(ii, k)
    488798          ENDDO
    489        ENDDO
     799        ENDDO
     800      ENDDO
    490801#endif
    491802
    … …  
    494805    DO j = jfiltsu, jjm
    495806
    496        DO i = 1,iim
    497           coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
    498           IF( i<modfrstu(j) ) coff = 0.
    499           DO k=1,iim
    500              eignft(i,k) = eignfnv(k,i) * coff
    501           ENDDO
    502        ENDDO
     807      DO i = 1, iim
     808        coff = coefilu(i, j) / (1. + coefilu(i, j))
     809        IF(i<modfrstu(j)) coff = 0.
     810        DO k = 1, iim
     811          eignft(i, k) = eignfnv(k, i) * coff
     812        ENDDO
     813      ENDDO
    503814#ifdef BLAS
    504815       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
    505816            eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
    506817#else
    507        DO k = 1, iim
    508           DO i = 1, iim
    509              matrinvs(i,k,j-jfiltsu+1) = 0.0
    510              DO ii = 1, iim
    511                 matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
    512                      + eignfnv(i,ii)*eignft(ii,k)
    513              ENDDO
     818      DO k = 1, iim
     819        DO i = 1, iim
     820          matrinvs(i, k, j - jfiltsu + 1) = 0.0
     821          DO ii = 1, iim
     822            matrinvs(i, k, j - jfiltsu + 1) = matrinvs(i, k, j - jfiltsu + 1) &
     823                    + eignfnv(i, ii) * eignft(ii, k)
    514824          ENDDO
    515        ENDDO
     825        ENDDO
     826      ENDDO
    516827#endif
    517828
    … …  
    528839    !   ...................................................................
    529840
    530 334 FORMAT(1x,24i3)
     841    334 FORMAT(1x, 24i3)
    531842
    532843  END SUBROUTINE inifilr
    533844
    534 END MODULE filtreg_mod
     845END MODULE lmdz_filtreg
Note: See TracChangeset for help on using the changeset viewer.