Ignore:
Timestamp:
Feb 3, 2009, 10:50:31 AM (16 years ago)
Author:
yann meurdesoif
Message:

Modifications Othman Bouzi : optimisation du filtre (remplacement opération matrice/vecteurs par matrice/matrice/matrice - BLAS), l'allocation des tableaux du filtre se fait maintenant dynamiquement (plus d'intervention manuelle dans parafiltre.h)

Location:
LMDZ4/branches/LMDZ4-dev/libf/filtrez
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/branches/LMDZ4-dev/libf/filtrez/coefils.h

    r524 r1086  
    22! $Header$
    33!
    4       COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)
    5      * ,unsddu(iim),unsddv(iim),coefilu(iim,jjm),coefilv(iim,jjm),
    6      * modfrstu(jjm),modfrstv(jjm),eignfnu(iim,iim),eignfnv(iim,iim)
    7      * ,coefilu2(iim,jjm),coefilv2(iim,jjm)
    8 c
     4      COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)&
     5     & ,unsddu(iim),unsddv(iim),coefilu(iim,jjm),coefilv(iim,jjm),      &
     6     & modfrstu(jjm),modfrstv(jjm),eignfnu(iim,iim),eignfnv(iim,iim)    &
     7     & ,coefilu2(iim,jjm),coefilv2(iim,jjm)
     8!c
    99      INTEGER jfiltnu,jfiltsu,jfiltnv,jfiltsv,modfrstu,modfrstv
    1010      REAL    sddu,sddv,unsddu,unsddv,coefilu,coefilv,eignfnu,eignfnv
  • LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F

    r524 r1086  
    33!
    44      SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
    5      .   griscal ,iter)
    6 
     5     &     griscal ,iter)
     6     
     7      USE filtre
     8     
    79      IMPLICIT NONE
    810c=======================================================================
     
    4648#include "dimensions.h"
    4749#include "paramet.h"
    48 #include "parafilt.h"
    4950#include "coefils.h"
    50 c
    51       INTEGER nlat,nbniv,ifiltre,iter
    52       INTEGER i,j,l,k
    53       INTEGER iim2,immjm
    54       INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
    55 
    56       REAL  champ( iip1,nlat,nbniv)
    57       REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
    58       COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
    59      ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
    60      ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
    61       REAL  eignq(iim), sdd1(iim),sdd2(iim)
     51
     52      INTEGER    nlat,nbniv,ifiltre,iter
     53      INTEGER    i,j,l,k
     54      INTEGER    iim2,immjm
     55      INTEGER    jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
     56
     57      REAL       champ( iip1,nlat,nbniv)
     58
     59      REAL       eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
    6260      LOGICAL    griscal
    6361      INTEGER    hemisph, iaire
    64 c
     62
     63      LOGICAL,SAVE     :: first=.TRUE.
     64
     65      REAL, SAVE :: sdd12(iim,4)
     66
     67      INTEGER, PARAMETER :: type_sddu=1
     68      INTEGER, PARAMETER :: type_sddv=2
     69      INTEGER, PARAMETER :: type_unsddu=3
     70      INTEGER, PARAMETER :: type_unsddv=4
     71
     72      INTEGER :: sdd1_type, sdd2_type
     73
     74      IF (first) THEN
     75         sdd12(1:iim,type_sddu) = sddu(1:iim)
     76         sdd12(1:iim,type_sddv) = sddv(1:iim)
     77         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
     78         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
     79
     80         first=.FALSE.
     81      ENDIF
    6582
    6683      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
    67      *    STOP'Pas de transformee simple dans cette version'
    68 
     84     &     STOP'Pas de transformee simple dans cette version'
     85     
    6986      IF( iter.EQ. 2 )  THEN
    70        PRINT *,' Pas d iteration du filtre dans cette version !'
    71      * , ' Utiliser old_filtreg et repasser !'
    72            STOP
    73       ENDIF
    74 
     87         PRINT *,' Pas d iteration du filtre dans cette version !'
     88     &        , ' Utiliser old_filtreg et repasser !'
     89         STOP
     90      ENDIF
     91     
    7592      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
    76        PRINT *,' Cette routine ne calcule le filtre inverse que ',
    77      * ' sur la grille des scalaires !'
    78            STOP
    79       ENDIF
    80 
     93         PRINT *,' Cette routine ne calcule le filtre inverse que '
     94     &        , ' sur la grille des scalaires !'
     95         STOP
     96      ENDIF
     97     
    8198      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
    82        PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
    83      *,' corriger et repasser !'
    84            STOP
    85       ENDIF
    86 c
    87 
     99         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
     100     &        , ' corriger et repasser !'
     101         STOP
     102      ENDIF
     103     
    88104      iim2   = iim * iim
    89105      immjm  = iim * jjm
    90 c
    91 c
     106
    92107      IF( griscal )   THEN
    93108         IF( nlat. NE. jjp1 )  THEN
    94              PRINT  1111
    95              STOP
    96          ELSE
    97 c
    98              IF( iaire.EQ.1 )  THEN
    99                 CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
    100                 CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
    101              ELSE
    102                 CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
    103                 CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
    104              END IF
    105 c
    106              jdfil1 = 2
    107              jffil1 = jfiltnu
    108              jdfil2 = jfiltsu
    109              jffil2 = jjm
    110           END IF
     109            PRINT  1111
     110            STOP
     111         ELSE
     112           
     113            IF( iaire.EQ.1 )  THEN
     114               sdd1_type = type_sddu
     115               sdd2_type = type_unsddu
     116            ELSE
     117               sdd1_type = type_unsddu
     118               sdd2_type = type_sddu
     119            ENDIF
     120
     121c            IF( iaire.EQ.1 )  THEN
     122c               CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
     123c               CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
     124c            ELSE
     125c               CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
     126c               CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
     127c            END IF
     128           
     129            jdfil1 = 2
     130            jffil1 = jfiltnu
     131            jdfil2 = jfiltsu
     132            jffil2 = jjm
     133         END IF
    111134      ELSE
    112           IF( nlat.NE.jjm )  THEN
    113              PRINT  2222
    114              STOP
    115           ELSE
    116 c
    117              IF( iaire.EQ.1 )  THEN
    118                 CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
    119                 CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
    120              ELSE
    121                 CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
    122                 CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
    123              END IF
    124 c
    125              jdfil1 = 1
    126              jffil1 = jfiltnv
    127              jdfil2 = jfiltsv
    128              jffil2 = jjm
    129           END IF
     135         IF( nlat.NE.jjm )  THEN
     136            PRINT  2222
     137            STOP
     138         ELSE
     139           
     140            IF( iaire.EQ.1 )  THEN
     141               sdd1_type = type_sddu
     142               sdd2_type = type_unsddu
     143            ELSE
     144               sdd1_type = type_unsddu
     145               sdd2_type = type_sddu
     146            ENDIF
     147
     148c            IF( iaire.EQ.1 )  THEN
     149c               CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
     150c               CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
     151c            ELSE
     152c               CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
     153c               CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
     154c            END IF
     155           
     156            jdfil1 = 1
     157            jffil1 = jfiltnv
     158            jdfil2 = jfiltsv
     159            jffil2 = jjm
     160         END IF
    130161      END IF
    131 c
    132 c
    133       DO 100  hemisph = 1, 2
    134 c
    135       IF ( hemisph.EQ.1 )  THEN
    136           jdfil = jdfil1
    137           jffil = jffil1
    138       ELSE
    139           jdfil = jdfil2
    140           jffil = jffil2
    141       END IF
    142 
    143  
    144       DO 50  l = 1, nbniv
    145       DO 30  j = jdfil,jffil
    146  
    147  
    148       DO  5  i = 1, iim
    149       champ(i,j,l) = champ(i,j,l) * sdd1(i)
    150    5  CONTINUE
    151 c
    152 
    153       IF( hemisph. EQ. 1 )      THEN
    154 
    155         IF( ifiltre. EQ. -2 )   THEN
    156 #ifdef CRAY
    157          CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq  ,
    158      *                             1, iim, iim                         )
    159 #else
    160 #ifdef BLAS
    161       CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
    162      .           champ(1,j,l), 1, 0.0, eignq, 1)
    163 #else
    164       DO k = 1, iim
    165          eignq(k) = 0.0
     162     
     163      DO hemisph = 1, 2
     164         
     165         IF ( hemisph.EQ.1 )  THEN
     166            jdfil = jdfil1
     167            jffil = jffil1
     168         ELSE
     169            jdfil = jdfil2
     170            jffil = jffil2
     171         END IF
     172         
     173         DO l = 1, nbniv
     174            DO j = jdfil,jffil
     175               DO i = 1, iim
     176                  champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
     177               END DO
     178            END DO
     179         END DO
     180         
     181         IF( hemisph. EQ. 1 )      THEN
     182           
     183            IF( ifiltre. EQ. -2 )   THEN
     184               
     185               DO j = jdfil,jffil
     186                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0,
     187     &                 matrinvn(1,1,j),
     188     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
     189     &                 eignq(1,j-jdfil+1,1), iim*nlat)
     190               END DO
     191               
     192            ELSE IF ( griscal )     THEN
     193               
     194               DO j = jdfil,jffil
     195                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0,
     196     &                 matriceun(1,1,j),
     197     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
     198     &                 eignq(1,j-jdfil+1,1), iim*nlat)
     199               END DO
     200               
     201            ELSE
     202               
     203               DO j = jdfil,jffil
     204                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0,
     205     &                 matricevn(1,1,j),
     206     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
     207     &                 eignq(1,j-jdfil+1,1), iim*nlat)
     208               END DO
     209               
     210            ENDIF
     211           
     212         ELSE
     213           
     214            IF( ifiltre. EQ. -2 )   THEN
     215               
     216               DO j = jdfil,jffil
     217                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0,
     218     &                 matrinvs(1,1,j-jfiltsu+1),
     219     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
     220     &                 eignq(1,j-jdfil+1,1), iim*nlat)
     221               END DO
     222               
     223               
     224            ELSE IF ( griscal )     THEN
     225               
     226               DO j = jdfil,jffil
     227                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0,
     228     &                 matriceus(1,1,j-jfiltsu+1),
     229     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
     230     &                 eignq(1,j-jdfil+1,1), iim*nlat)
     231               END DO
     232                             
     233            ELSE
     234               
     235               DO j = jdfil,jffil
     236                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0,
     237     &                 matricevs(1,1,j-jfiltsv+1),
     238     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
     239     &                 eignq(1,j-jdfil+1,1), iim*nlat)
     240               END DO
     241                             
     242            ENDIF
     243           
     244         ENDIF
     245         
     246         IF( ifiltre.EQ. 2 )  THEN
     247           
     248            DO l = 1, nbniv
     249               DO j = jdfil,jffil
     250                  DO i = 1, iim
     251                     champ( i,j,l ) =
     252     &                    (champ(i,j,l) + eignq(i,j-jdfil+1,l))
     253     &                    * sdd12(i,sdd2_type) ! sdd2(i)
     254                  END DO
     255               END DO
     256            END DO
     257
     258         ELSE
     259
     260            DO l = 1, nbniv
     261               DO j = jdfil,jffil
     262                  DO i = 1, iim
     263                     champ( i,j,l ) =
     264     &                    (champ(i,j,l) - eignq(i,j-jdfil+1,l))
     265     &                    * sdd12(i,sdd2_type) ! sdd2(i)
     266                  END DO
     267               END DO
     268            END DO
     269
     270         ENDIF
     271
     272         DO l = 1, nbniv
     273            DO j = jdfil,jffil
     274               champ( iip1,j,l ) = champ( 1,j,l )
     275            END DO
     276         END DO
     277
     278     
    166279      ENDDO
    167       DO k = 1, iim
    168       DO i = 1, iim
    169          eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
    170       ENDDO
    171       ENDDO
    172 #endif
    173 #endif
    174         ELSE IF ( griscal )     THEN
    175 #ifdef CRAY
    176          CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
    177      *                             1, iim, iim                         )
    178 #else
    179 #ifdef BLAS
    180       CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
    181      .           champ(1,j,l), 1, 0.0, eignq, 1)
    182 #else
    183       DO k = 1, iim
    184          eignq(k) = 0.0
    185       ENDDO
    186       DO i = 1, iim
    187       DO k = 1, iim
    188          eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
    189       ENDDO
    190       ENDDO
    191 #endif
    192 #endif
    193         ELSE
    194 #ifdef CRAY
    195          CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
    196      *                             1, iim, iim                         )
    197 #else
    198 #ifdef BLAS
    199       CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
    200      .           champ(1,j,l), 1, 0.0, eignq, 1)
    201 #else
    202       DO k = 1, iim
    203          eignq(k) = 0.0
    204       ENDDO
    205       DO i = 1, iim
    206       DO k = 1, iim
    207          eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
    208       ENDDO
    209       ENDDO
    210 #endif
    211 #endif
    212         ENDIF
    213 
    214       ELSE
    215 
    216         IF( ifiltre. EQ. -2 )   THEN
    217 #ifdef CRAY
    218          CALL MXVA( matrinvs(1,1,j-jfiltsu+1),  1, iim, champ(1,j,l),1 , 
    219      *                          eignq,  1, iim, iim                    )
    220 #else
    221 #ifdef BLAS
    222       CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
    223      .           champ(1,j,l), 1, 0.0, eignq, 1)
    224 #else
    225       DO k = 1, iim
    226          eignq(k) = 0.0
    227       ENDDO
    228       DO i = 1, iim
    229       DO k = 1, iim
    230          eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
    231       ENDDO
    232       ENDDO
    233 #endif
    234 #endif
    235         ELSE IF ( griscal )     THEN
    236 #ifdef CRAY
    237          CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
    238      *                          eignq,  1, iim, iim                    )
    239 #else
    240 #ifdef BLAS
    241       CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
    242      .           champ(1,j,l), 1, 0.0, eignq, 1)
    243 #else
    244       DO k = 1, iim
    245          eignq(k) = 0.0
    246       ENDDO
    247       DO i = 1, iim
    248       DO k = 1, iim
    249          eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
    250       ENDDO
    251       ENDDO
    252 #endif
    253 #endif
    254         ELSE
    255 #ifdef CRAY
    256          CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
    257      *                          eignq,  1, iim, iim                    )
    258 #else
    259 #ifdef BLAS
    260       CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
    261      .           champ(1,j,l), 1, 0.0, eignq, 1)
    262 #else
    263       DO k = 1, iim
    264          eignq(k) = 0.0
    265       ENDDO
    266       DO i = 1, iim
    267       DO k = 1, iim
    268          eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
    269       ENDDO
    270       ENDDO
    271 #endif
    272 #endif
    273         ENDIF
    274 
    275       ENDIF
    276 c
    277       IF( ifiltre.EQ. 2 )  THEN
    278         DO 15 i = 1, iim
    279         champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
    280   15    CONTINUE
    281       ELSE
    282         DO 16 i=1,iim
    283         champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
    284 16      CONTINUE
    285       ENDIF
    286 c
    287       champ( iip1,j,l ) = champ( 1,j,l )
    288 c
    289   30  CONTINUE
    290 c
    291   50  CONTINUE
    292 c   
    293  100  CONTINUE
    294 c
     280
    2952811111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
    296      *filtrer, sur la grille des scalaires'/)
     282     &     filtrer, sur la grille des scalaires'/)
    2972832222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
    298      *ltrer, sur la grille de V ou de Z'/)
     284     &     ltrer, sur la grille de V ou de Z'/)
    299285      RETURN
    300286      END
  • LMDZ4/branches/LMDZ4-dev/libf/filtrez/inifgn.F

    r524 r1086  
    11!
    2 ! $Header$
     2! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004-05-19 12:53:09 lmdzadmin Exp $
    33!
    44      SUBROUTINE inifgn(dv)
  • LMDZ4/branches/LMDZ4-dev/libf/filtrez/parafilt.h

    r1024 r1086  
    33!
    44        INTEGER nfilun, nfilus, nfilvn, nfilvs
    5 
    6       PARAMETER (nfilun=24, nfilus=23, nfilvn=24, nfilvs=24)
    7 
    8 c
    9 c
    10 c      Ici , on a exagere  les nombres de lignes de latitudes a filtrer .
    11 c
    12 c      La premiere fois que  le Gcm  rentrera  dans le Filtre ,
    13 c
    14 c      il indiquera  les bonnes valeurs  de  nfilun , nflius, nfilvn  et
    15 c
    16 c      nfilvs  a  mettre .  Il suffira alors de changer ces valeurs dans
    17 c
    18 c      Parameter  ci-dessus  et de relancer  le  run . 
    19 
Note: See TracChangeset for help on using the changeset viewer.