Changeset 1601


Ignore:
Timestamp:
Dec 7, 2011, 8:05:43 AM (13 years ago)
Author:
Ehouarn Millour
Message:

Updates and upgrades of filter:

  • add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
  • adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
  • update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
  • update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM

Location:
LMDZ4/branches/LMDZ4_AR5
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4_AR5/arch/arch-PW6_VARGAS.fcm

    r1459 r1601  
    33%AR                  ar
    44%MAKE                gmake
    5 %FPP_FLAGS           -P
    6 %FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM
    7 %BASE_FFLAGS         -qautodbl=dbl4 -qxlf90=autodealloc -qextname=flush
     5%FPP_FLAGS           -P -I/usr/local/pub/FFTW/3.2/include
     6%FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_FFTW
     7%BASE_FFLAGS         -qautodbl=dbl4 -qxlf90=autodealloc -qmaxmem=-1 -qzerosize -I/usr/local/pub/FFTW/3.2/include
    88%PROD_FFLAGS         -O3
    99%DEV_FFLAGS          -O2 -qfullpath -qinitauto=7FBFFFFF -qfloat=nans -qflttrap=overflow:zerodivide:invalid:enable -qsigtrap
    … …  
    1111%MPI_FFLAGS          -I/usr/lpp/ppe.poe/include/thread64
    1212%OMP_FFLAGS          -qsmp=omp
    13 %BASE_LD             -lessl
     13%BASE_LD             -lessl -L/usr/local/pub/FFTW/3.2/lib -lfftw3
    1414%MPI_LD             
    1515%OMP_LD              -qsmp=omp

  • LMDZ4/branches/LMDZ4_AR5/libf/dyn3dpar/filtreg_p.F

    r1146 r1601  
    208208               IF( ifiltre.EQ.-2 )   THEN
    209209                  DO j = jdfil,jffil
     210#ifdef BLAS
    210211                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    211212     &                    matrinvn(1,1,j), iim,
    212213     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    213214     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     215#else
     216                     champ_fft(:,j-jdfil+1,:)
     217     &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
     218#endif
    214219                  ENDDO
    215220                 
    216221               ELSE IF ( griscal )     THEN
    217222                  DO j = jdfil,jffil
     223#ifdef BLAS
    218224                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    219225     &                    matriceun(1,1,j), iim,
    220226     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    221227     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     228#else
     229                     champ_fft(:,j-jdfil+1,:)
     230     &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
     231#endif
    222232                  ENDDO
    223233                 
    224234               ELSE
    225235                  DO j = jdfil,jffil
     236#ifdef BLAS
    226237                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    227238     &                    matricevn(1,1,j), iim,
    228239     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    229240     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     241#else
     242                     champ_fft(:,j-jdfil+1,:)
     243     &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
     244#endif
    230245                  ENDDO
    231246                 
    … …  
    236251               IF( ifiltre.EQ.-2 )   THEN
    237252                  DO j = jdfil,jffil
     253#ifdef BLAS
    238254                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    239255     &                    matrinvs(1,1,j-jfiltsu+1), iim,
    240256     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    241257     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     258#else
     259                     champ_fft(:,j-jdfil+1,:)
     260     &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
     261     &                            champ_loc(:iim,j,:))
     262#endif
    242263                  ENDDO
    243264                 
    … …  
    245266                 
    246267                  DO j = jdfil,jffil
     268#ifdef BLAS
    247269                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    248270     &                    matriceus(1,1,j-jfiltsu+1), iim,
    249271     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    250272     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     273#else
     274                     champ_fft(:,j-jdfil+1,:)
     275     &                    =matmul(matriceus(:,:,j-jfiltsu+1),
     276     &                            champ_loc(:iim,j,:))
     277#endif
    251278                  ENDDO
    252279                 
    … …  
    254281                 
    255282                  DO j = jdfil,jffil
     283#ifdef BLAS
    256284                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    257285     &                    matricevs(1,1,j-jfiltsv+1), iim,
    258286     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    259287     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     288#else
     289                     champ_fft(:,j-jdfil+1,:)
     290     &                    =matmul(matricevs(:,:,j-jfiltsv+1),
     291     &                            champ_loc(:iim,j,:))
     292#endif
    260293                  ENDDO
    261294                 

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/coefils.h

    r1146 r1601  
    11!
    2 ! $Header$
     2! $Id $
    33!
    44      COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)&
    … …  
    77     & ,coefilu2(iim,jjm),coefilv2(iim,jjm)
    88!c
    9       INTEGER jfiltnu,jfiltsu,jfiltnv,jfiltsv,modfrstu,modfrstv
     9      INTEGER jfiltnu ! index of the last lat line filtered in NH (U grid)
     10      INTEGER jfiltsu ! index of the first lat line filtered in SH (U grid)
     11      INTEGER jfiltnv ! index of the last lat line filtered in NH (V grid)
     12      INTEGER jfiltsv ! index of the first lat line filtered in SH (V grid)
     13      INTEGER modfrstu ! number of retained (ie: unfiltered) modes on U grid
     14      INTEGER modfrstv ! number of retained (ie: unfiltered) modes on V grid
    1015      REAL    sddu,sddv,unsddu,unsddv,coefilu,coefilv,eignfnu,eignfnv
    1116      REAL    coefilu2,coefilv2

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/filtreg_mod.F90

    r1279 r1601  
     1!
     2! $Id $
     3!
    14MODULE filtreg_mod
    25
    … …  
    4245    INTEGER ixmineq
    4346#endif
    44     EXTERNAL  inifgn
    4547    !
    4648    ! ------------------------------------------------------------
    … …  
    7173    CALL inifgn(eignvl)
    7274    !
    73     PRINT *,' EIGNVL '
     75    PRINT *,'inifilr: EIGNVL '
    7476    PRINT 250,eignvl
    75 250 FORMAT( 1x,5e13.6)
     77250 FORMAT( 1x,5e14.6)
    7678    !
    7779    ! compute eigenvalues and eigenfunctions
    … …  
    113115#endif
    114116    !
     117    ! For a regular grid, we want the filter to start at latitudes
     118    ! corresponding to lengths dx of the same size as dy (in terms
     119    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
     120    !  <=> latitude=60 degrees).
     121    ! Same idea for the zoomed grid: start filtering polewards as soon
     122    ! as length dx becomes of the same size as dy
    115123    !
    116124    colat0  =  MIN( 0.5, dymin/dxmin )
    … …  
    158166    imx  = iim
    159167    !
    160     PRINT *,' TRUNCATION AT ',imx
    161     !
     168    PRINT *,'inifilr: TRUNCATION AT ',imx
     169    !
     170! Ehouarn: set up some defaults
     171    jfiltnu=2 ! avoid north pole
     172    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
     173    jfiltnv=1 ! NB: no poles on the V grid
     174    jfiltsv=jjm
     175
    162176    DO j = 2, jjm/2+1
    163177       cof = COS( rlatu(j) )/ colat0
    164178       IF ( cof .LT. 1. ) THEN
    165           IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
     179          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
     180            jfiltnu= j
     181          ENDIF
    166182       ENDIF
    167183
    168184       cof = COS( rlatu(jjp1-j+1) )/ colat0
    169185       IF ( cof .LT. 1. ) THEN
    170           IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
     186          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
    171187               jfiltsu= jjp1-j+1
     188          ENDIF
    172189       ENDIF
    173190    ENDDO
    … …  
    176193       cof = COS( rlatv(j) )/ colat0
    177194       IF ( cof .LT. 1. ) THEN
    178           IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
     195          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
     196            jfiltnv= j
     197          ENDIF
    179198       ENDIF
    180199
    181200       cof = COS( rlatv(jjm-j+1) )/ colat0
    182201       IF ( cof .LT. 1. ) THEN
    183           IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
     202          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
    184203               jfiltsv= jjm-j+1
     204          ENDIF
    185205       ENDIF
    186206    ENDDO
    187207    !                                 
    188208
    189     IF ( jfiltnu.LE.0 ) jfiltnu=1
    190209    IF( jfiltnu.GT. jjm/2 +1 )  THEN
    191210       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
    … …  
    193212    ENDIF
    194213
    195     IF( jfiltsu.LE.0) jfiltsu=1
    196214    IF( jfiltsu.GT.  jjm  +1 )  THEN
    197215       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
    … …  
    199217    ENDIF
    200218
    201     IF( jfiltnv.LE.0) jfiltnv=1
    202219    IF( jfiltnv.GT. jjm/2    )  THEN
    203220       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
    … …  
    205222    ENDIF
    206223
    207     IF( jfiltsv.LE.0) jfiltsv=1
    208224    IF( jfiltsv.GT.     jjm  )  THEN
    209225       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
    … …  
    211227    ENDIF
    212228
    213     PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
     229    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
    214230         jfiltnv,jfiltsv,jfiltnu,jfiltsu
    215231
    216232    IF(first_call_inifilr) THEN
    217233       ALLOCATE(matriceun(iim,iim,jfiltnu))
    218        ALLOCATE(matriceus(iim,iim,jfiltsu))
     234       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
    219235       ALLOCATE(matricevn(iim,iim,jfiltnv))
    220        ALLOCATE(matricevs(iim,iim,jfiltsv))
     236       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
    221237       ALLOCATE( matrinvn(iim,iim,jfiltnu))
    222        ALLOCATE( matrinvs(iim,iim,jfiltsu))
     238       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
    223239       first_call_inifilr = .FALSE.
    224240    ENDIF
    … …  
    230246    !
    231247    DO j = 1,jjm
     248    !default initialization: all modes are retained (i.e. no filtering)
    232249       modfrstu( j ) = iim
    233250       modfrstv( j ) = iim
    … …  
    306323
    307324    IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
    308 
     325! Ehouarn: and what are these for??? Trying to handle a limit case
     326!          where filters extend to and meet at the equator?
    309327       IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
    310328       IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
    … …  
    334352             eignft(i,k) = eignfnv(k,i) * coff
    335353          ENDDO
    336        ENDDO
     354       ENDDO ! of DO i=1,iim
    337355#ifdef CRAY
    338356       CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
    … …  
    350368             ENDDO
    351369          ENDDO
    352        ENDDO
    353 #endif
    354 #endif
    355 
    356     ENDDO
     370       ENDDO ! of DO k = 1, iim
     371#endif
     372#endif
     373
     374    ENDDO ! of DO j = 2, jfiltnu
    357375
    358376    DO j = jfiltsu, jjm
    … …  
    364382             eignft(i,k) = eignfnv(k,i) * coff
    365383          ENDDO
    366        ENDDO
     384       ENDDO ! of DO i=1,iim
    367385#ifdef CRAY
    368386       CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
    … …  
    381399             ENDDO
    382400          ENDDO
    383        ENDDO
    384 #endif
    385 #endif
    386 
    387     ENDDO
     401       ENDDO ! of DO k = 1, iim
     402#endif
     403#endif
     404
     405    ENDDO ! of DO j = jfiltsu, jjm
    388406
    389407    !   ...................................................................
    … …  
    421439#endif
    422440
    423     ENDDO
     441    ENDDO ! of DO j = 1, jfiltnv
    424442
    425443    DO j = jfiltsv, jjm
    … …  
    452470#endif
    453471
    454     ENDDO
     472    ENDDO ! of DO j = jfiltsv, jjm
    455473
    456474    !   ...................................................................
    … …  
    488506#endif
    489507
    490     ENDDO
     508    ENDDO ! of DO j = 2, jfiltnu
    491509
    492510    DO j = jfiltsu, jjm
    … …  
    518536#endif
    519537
    520     ENDDO
     538    ENDDO ! of DO j = jfiltsu, jjm
    521539
    522540    IF (use_filtre_fft) THEN

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_fft_fftw.F90

    r986 r1601  
     1!
     2! $Id$
     3!
     4
    15MODULE mod_fft_fftw
    26
    37#ifdef FFT_FFTW
    48
    5   REAL,SAVE,ALLOCATABLE    :: Table_forward(:)
    6   REAL,SAVE,ALLOCATABLE    :: Table_backward(:)
    7   REAL,SAVE                :: scale_factor
    8   INTEGER,SAVE             :: vsize
    9   INTEGER,PARAMETER        :: inc=1
     9  REAL, SAVE                   :: scale_factor
     10  INTEGER, SAVE                :: vsize
     11  INTEGER, PARAMETER           :: inc=1
    1012 
    11   INTEGER,SAVE            :: plan_forward
    12   INTEGER,SAVE            :: plan_backward
     13  INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward
     14  INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward
    1315 
    1416CONTAINS
    1517 
    16   SUBROUTINE Init_fft(iim)
     18  SUBROUTINE Init_fft(iim,nvectmax)
    1719  IMPLICIT NONE
    18 #include <rfftw.h>
     20#include <fftw3.f>
    1921    INTEGER :: iim
    20     REAL    :: rtmp=1.
    21     COMPLEX*16 :: ctmp
    22     INTEGER :: itmp=1
    23     INTEGER :: isign=0
    24     INTEGER :: ierr
     22    INTEGER :: nvectmax
     23
     24    INTEGER :: itmp
     25
     26    INTEGER               :: rank
     27    INTEGER               :: howmany
     28    INTEGER               :: istride, idist
     29    INTEGER               :: ostride, odist
     30    INTEGER, DIMENSION(1) :: n_array, inembed, onembed
     31
     32    REAL,    DIMENSION(iim+1,nvectmax) :: dbidon
     33    COMPLEX, DIMENSION(iim/2+1,nvectmax) :: cbidon
     34
     35    vsize = iim
     36    scale_factor = 1./SQRT(1.*vsize)
     37
     38    dbidon = 0
     39    cbidon = 0
     40
     41    ALLOCATE(plan_forward(nvectmax))
     42    ALLOCATE(plan_backward(nvectmax))
    2543   
    26     vsize=iim
    27     scale_factor=1./SQRT(1.*vsize)
    28     ALLOCATE(Table_forward(2*vsize+64))
    29     ALLOCATE(Table_backward(2*vsize+64))
     44    WRITE(*,*)"!---------------------!"
     45    WRITE(*,*)"!                     !"
     46    WRITE(*,*)"! INITIALISATION FFTW !"
     47    WRITE(*,*)"!                     !"
     48    WRITE(*,*)"!---------------------!"
    3049   
    31 !    CALL DZFFTM(isign,vsize,itmp,scale_factor,rtmp,vsize+inc,ctmp,vsize/2+1,table_forward,rtmp,ierr)
    32    
    33 !    CALL ZDFFTM(isign,vsize,itmp,scale_factor,ctmp,vsize/2+1,rtmp,vsize+inc,table_backward,rtmp,ierr)
     50! On initialise tous les plans
     51    DO itmp = 1, nvectmax
     52       rank       = 1
     53       n_array(1) = iim
     54       howmany    = itmp
     55       inembed(1) = iim + 1 ; onembed(1) = iim/2 + 1
     56       istride    = 1       ; ostride    = 1
     57       idist      = iim + 1 ; odist      = iim/2 + 1
    3458
    35     CALL rfftw_f77_create_plan(plan_forward,iim,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE)
    36     CALL rfftw_f77_create_plan(plan_backward,iim,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE)
    37    
     59       CALL dfftw_plan_many_dft_r2c(plan_forward(itmp), rank, n_array, howmany, &
     60            & dbidon, inembed, istride, idist, &
     61            & cbidon, onembed, ostride, odist, &
     62            & FFTW_ESTIMATE)
     63
     64       rank       = 1
     65       n_array(1) = iim
     66       howmany    = itmp
     67       inembed(1) = iim/2 + 1 ; onembed(1) = iim + 1
     68       istride    = 1         ; ostride    = 1
     69       idist      = iim/2 + 1 ; odist      = iim + 1
     70       CALL dfftw_plan_many_dft_c2r(plan_backward(itmp), rank, n_array, howmany, &
     71            & cbidon, inembed, istride, idist, &
     72            & dbidon, onembed, ostride, odist, &
     73            & FFTW_ESTIMATE)
     74
     75    ENDDO
     76
     77    WRITE(*,*)"!-------------------------!"
     78    WRITE(*,*)"!                         !"
     79    WRITE(*,*)"! FIN INITIALISATION FFTW !"
     80    WRITE(*,*)"!                         !"
     81    WRITE(*,*)"!-------------------------!"
     82
    3883  END SUBROUTINE Init_fft
    3984 
    … …  
    4186  SUBROUTINE fft_forward(vect,TF_vect,nb_vect)
    4287    IMPLICIT NONE
    43 #include <rfftw.h>
    44     INTEGER,INTENT(IN)  :: nb_vect
    45     REAL,INTENT(IN)     :: vect(vsize+inc,nb_vect)
    46     COMPLEX*16,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect)
    47     REAL                :: work(4*vsize*nb_vect)
    48     INTEGER             :: ierr
    49     INTEGER, PARAMETER :: isign=-1
    50  
    51 !    CALL DZFFTM(isign,vsize,nb_vect,scale_factor,vect,vsize+inc,TF_vect,vsize/2+1,table_forward,work,ierr)
    52      CALL rfftwnd_f77_real_to_complex(plan_forward,nb_vect,vect, 1, vsize+inc , TF_vect, 1, vsize/2+1); 
    53    
     88#include <fftw3.f>
     89    INTEGER,INTENT(IN)     :: nb_vect
     90    REAL,INTENT(IN)        :: vect(vsize+inc,nb_vect)
     91    COMPLEX,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect)
     92
     93    CALL dfftw_execute_dft_r2c(plan_forward(nb_vect),vect,TF_vect)
     94
     95    TF_vect = scale_factor * TF_vect
     96
    5497  END SUBROUTINE fft_forward
    5598 
    5699  SUBROUTINE fft_backward(TF_vect,vect,nb_vect)
    57100    IMPLICIT NONE
    58 #include <rfftw.h>
    59     INTEGER,INTENT(IN)  :: nb_vect
    60     REAL,INTENT(OUT)    :: vect(vsize+inc,nb_vect)
    61     COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
    62     REAL                :: work(4*vsize*nb_vect)
    63     INTEGER             :: ierr
    64     INTEGER, PARAMETER :: isign=1
    65  
    66 !    CALL ZDFFTM(isign,vsize,nb_vect,scale_factor,TF_vect,vsize/2+1,vect,vsize+inc,table_backward,work,ierr)
    67     CALL rfftwnd_f77_complex_to_real(plan_forward,nb_vect,TF_vect, 1, vsize/2+1 , vect, 1, vsize+inc); 
     101#include <fftw3.f>
     102    INTEGER,INTENT(IN)     :: nb_vect
     103    REAL,INTENT(OUT)       :: vect(vsize+inc,nb_vect)
     104    COMPLEX,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
     105
     106    CALL dfftw_execute_dft_c2r(plan_backward(nb_vect),TF_vect,vect)
     107
     108    vect = scale_factor * vect
    68109
    69110  END SUBROUTINE fft_backward
    … …  
    72113 
    73114END MODULE mod_fft_fftw
    74 

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_fft_mathkeisan.F90

    r986 r1601  
    1515    INTEGER :: nb_vect_max
    1616    REAL    :: rtmp=1.
    17     COMPLEX*16 :: ctmp
     17    COMPLEX :: ctmp
    1818    INTEGER :: itmp=1
    1919    INTEGER :: isign=0
    … …  
    3737    INTEGER,INTENT(IN)  :: nb_vect
    3838    REAL,INTENT(IN)     :: vect(vsize+inc,nb_vect)
    39     COMPLEX*16,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect)
     39    COMPLEX,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect)
    4040    REAL                :: work(4*vsize*nb_vect)
    4141    INTEGER             :: ierr
    … …  
    5151    INTEGER,INTENT(IN)  :: nb_vect
    5252    REAL,INTENT(OUT)    :: vect(vsize+inc,nb_vect)
    53     COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
     53    COMPLEX,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
    5454    REAL                :: work(4*vsize*nb_vect)
    5555    INTEGER             :: ierr

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_fft_mkl.F90

    r1279 r1601  
    2424    INTEGER :: nb_vect_max
    2525    REAL    :: rtmp=1.
    26     COMPLEX*16 :: ctmp
     26    COMPLEX :: ctmp
    2727    INTEGER :: itmp=1
    2828    INTEGER :: isign=0
    … …  
    6060    INTEGER,INTENT(IN)  :: nb_vect
    6161    REAL,INTENT(IN)     :: vect((vsize+inc)*nb_vect)
    62     COMPLEX*16,INTENT(OUT) :: TF_vect((vsize/2+1)*nb_vect)
     62    COMPLEX,INTENT(OUT) :: TF_vect((vsize/2+1)*nb_vect)
    6363    REAL                :: work(4*vsize*nb_vect)
    6464    INTEGER             :: ierr
    … …  
    102102    INTEGER,INTENT(IN)  :: nb_vect
    103103    REAL,INTENT(OUT)    :: vect((vsize+inc)*nb_vect)
    104     COMPLEX*16,INTENT(IN ) :: TF_vect((vsize/2+1)*nb_vect)
     104    COMPLEX,INTENT(IN ) :: TF_vect((vsize/2+1)*nb_vect)
    105105    REAL                :: work(4*vsize*nb_vect)
    106106    INTEGER             :: ierr

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_fft_wrapper.F90

    r1279 r1601  
    1919    INTEGER,INTENT(IN)  :: nb_vect
    2020    REAL,INTENT(IN)     :: vect(vsize+inc,nb_vect)
    21     COMPLEX*16,INTENT(INOUT) :: TF_vect(vsize/2+1,nb_vect)
     21    COMPLEX,INTENT(INOUT) :: TF_vect(vsize/2+1,nb_vect)
    2222   
    2323    STOP "wrapper fft : une FFT doit etre specifiee a l'aide d'une clee CPP, sinon utiliser le filtre classique"
    … …  
    2929    INTEGER,INTENT(IN)  :: nb_vect
    3030    REAL,INTENT(INOUT)    :: vect(vsize+inc,nb_vect)
    31     COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
     31    COMPLEX,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
    3232 
    3333    STOP "wrapper fft : une FFT doit etre specifiee a l'aide d'une clee CPP, sinon utiliser le filtre classique"

  • LMDZ4/branches/LMDZ4_AR5/libf/filtrez/mod_filtre_fft.F90

    r1279 r1601  
     1!
     2! $Id$
     3!
     4
    15MODULE mod_filtre_fft
    26
    … …  
    2327    INTEGER            :: index_vp(iim)
    2428    INTEGER            :: i,j
    25    
     29    INTEGER            :: l,ll_nb
     30
    2631    index_vp(1)=1
    2732    DO i=1,iim/2
    … …  
    98103    ENDDO
    99104   
    100    
     105#ifdef FFT_FFTW
     106
     107    WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
     108    WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
     109    WRITE (*,*)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
     110    CALL Init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
     111#else   
    101112    CALL Init_FFT(iim,(jjm+1)*(llm+1))
    102        
     113#endif       
    103114   
    104115  END SUBROUTINE Init_filtre_fft
    … …  
    118129
    119130    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
    120 !    REAL               :: vect_test(iim+inc,jj_end-jj_begin+1,nbniv)
    121     COMPLEX*16         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
    122 !    COMPLEX*16         :: TF_vect_test(iim/2+1,jj_end-jj_begin+1,nbniv)
     131    COMPLEX         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
    123132    INTEGER            :: nb_vect
    124133    INTEGER :: i,j,l
    125134    INTEGER :: ll_nb
    126 !    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
    127135   
    128136    ll_nb=0
    … …  
    140148    nb_vect=(jj_end-jj_begin+1)*ll_nb
    141149
    142 !    vect_tmp=vect
    143 
    144150    CALL FFT_forward(vect,TF_vect,nb_vect)
    145 
    146 !    CALL FFT_forward(vect,TF_vect_test,nb_vect)
    147 !      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
    148 !      DO j=1,jj_end-jj_begin+1
    149 !      DO i=1,iim/2+1
    150 !         PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)
    151 !       ENDDO
    152 !      ENDDO
    153151
    154152    DO l=1,ll_nb
    … …  
    159157      ENDDO
    160158    ENDDO
    161        
     159 
    162160    CALL FFT_backward(TF_vect,vect,nb_vect)
    163 !    CALL FFT_backward(TF_vect_test,vect_test,nb_vect)
    164          
    165 !      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
    166 !      DO j=1,jj_end-jj_begin+1
    167 !         DO i=1,iim
    168 !           PRINT *,"====",i,j,"----->",vect_test(i,j,1)
    169 !         ENDDO
    170 !      ENDDO
    171 
     161     
     162     
    172163    ll_nb=0
    173164!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    … …  
    199190
    200191    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
    201     COMPLEX*16         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
     192    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
    202193    INTEGER            :: nb_vect
    203194    INTEGER :: i,j,l
    … …  
    260251    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
    261252
    262     REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
    263     COMPLEX*16         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
     253     REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
     254    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
    264255    INTEGER            :: nb_vect
    265256    INTEGER :: i,j,l
    … …  
    305296
    306297  END SUBROUTINE Filtre_inv_fft 
    307  
    308  
    309 !  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
    310 !  IMPLICIT NONE
    311 !    INTEGER,INTENT(IN)  :: nbniv
    312 !    INTEGER,INTENT(OUT) :: ll_index(nbniv)
    313 !    INTEGER,INTENT(OUT) :: ll_nb
    314 !
    315 !    INTEGER :: l,ll_begin, ll_end
    316 !   INTEGER :: omp_rank,omp_size
    317 !   INTEGER :: OMP_GET_NUM_THREADS
    318 !   INTEGER :: omp_chunk
    319 !   EXTERNAL OMP_GET_NUM_THREADS
    320 !   INTEGER :: OMP_GET_THREAD_NUM
    321 !   EXTERNAL OMP_GET_THREAD_NUM
    322 !
    323 !   
    324 !   omp_size=OMP_GET_NUM_THREADS()
    325 !   omp_rank=OMP_GET_THREAD_NUM()   
    326 !   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
    327 !   
    328 !   ll_begin=omp_rank*OMP_CHUNK+1
    329 !   ll_nb=0
    330 !   DO WHILE (ll_begin<=nbniv)
    331 !     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
    332 !     DO l=ll_begin,ll_end
    333 !       ll_nb=ll_nb+1
    334 !       ll_index(ll_nb)=l
    335 !     ENDDO
    336 !     ll_begin=ll_begin+omp_size*OMP_CHUNK
    337 !   ENDDO
    338 ! 
    339 !  END SUBROUTINE get_ll_index
    340298   
    341299END MODULE mod_filtre_fft
Note: See TracChangeset for help on using the changeset viewer.