Ignore:
Timestamp:
Jul 23, 2024, 7:14:34 PM (8 weeks ago)
Author:
abarral
Message:

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3dmem/filtreg_p.F90

    r5104 r5105  
    11
    22
    3       SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
    4      &     ifiltre, iaire, griscal ,iter)
    5       USE parallel_lmdz, ONLY: OMP_CHUNK
    6       USE mod_filtre_fft
    7       USE timer_filtre
    8      
    9       USE filtreg_mod
    10      
    11       IMPLICIT NONE
    12      
    13 c=======================================================================
    14 c
    15 c   Auteur: P. Le Van        07/10/97
    16 c   ------
    17 c
    18 c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
    19 c                     pour l'operateur  Filtre    .
    20 c   ------
    21 c
    22 c   Arguments:
    23 c   ----------
    24 c
    25 c     
    26 c      ibeg..iend            lattitude a filtrer
    27 c      nlat                  nombre de latitudes du champ
    28 c      nbniv                 nombre de niveaux verticaux a filtrer
    29 c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
    30 c                            en sortie : champ filtre
    31 c      ifiltre               +1  Transformee directe
    32 c                            -1  Transformee inverse
    33 c                            +2  Filtre directe
    34 c                            -2  Filtre inverse
    35 c
    36 c      iaire                 1   si champ intensif
    37 c                            2   si champ extensif (pondere par les aires)
    38 c
    39 c      iter                  1   filtre simple
    40 c
    41 c=======================================================================
    42 c
    43 c
    44 c                      Variable Intensive
    45 c                ifiltre = 1     filtre directe
    46 c                ifiltre =-1     filtre inverse
    47 c
    48 c                      Variable Extensive
    49 c                ifiltre = 2     filtre directe
    50 c                ifiltre =-2     filtre inverse
    51 c
    52 c
    53       INCLUDE "dimensions.h"
    54       INCLUDE "paramet.h"
    55       INCLUDE "coefils.h"
    56 c
    57       INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter
    58       INTEGER i,j,l,k
    59       INTEGER iim2,immjm
    60       INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
    61      
    62       REAL  champ( iip1,nlat,nbniv)
    63      
    64       LOGICAL    griscal
    65       INTEGER    hemisph, iaire
    66      
    67       REAL :: champ_fft(iip1,nlat,nbniv)
    68       REAL :: champ_in(iip1,nlat,nbniv)
    69      
    70       LOGICAL,SAVE     :: first=.TRUE.
    71 c$OMP THREADPRIVATE(first)
    72 
    73       REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
    74       INTEGER :: ll_nb, nbniv_loc
    75       REAL, SAVE :: sdd12(iim,4)
    76 c$OMP THREADPRIVATE(sdd12)
    77 
    78       INTEGER, PARAMETER :: type_sddu=1
    79       INTEGER, PARAMETER :: type_sddv=2
    80       INTEGER, PARAMETER :: type_unsddu=3
    81       INTEGER, PARAMETER :: type_unsddv=4
    82 
    83       INTEGER :: sdd1_type, sdd2_type
    84 
    85       IF (first) THEN
    86          sdd12(1:iim,type_sddu) = sddu(1:iim)
    87          sdd12(1:iim,type_sddv) = sddv(1:iim)
    88          sdd12(1:iim,type_unsddu) = unsddu(1:iim)
    89          sdd12(1:iim,type_unsddv) = unsddv(1:iim)
    90 
    91          CALL Init_timer
    92          first=.FALSE.
    93       ENDIF
    94 
    95 c$OMP MASTER     
    96       CALL start_timer
    97 c$OMP END MASTER
    98 
    99 c-------------------------------------------------------c
    100 
    101       IF(ifiltre==1.or.ifiltre==-1)
    102      & CALL abort_gcm("fitreg_p","Pas de transformee simple
    103      &dans cette version",1)
    104      
    105       IF( iter== 2 )  THEN
    106          PRINT *,' Pas d iteration du filtre dans cette version !'
    107      &        , ' Utiliser old_filtreg et repasser !'
    108          CALL abort_gcm("fitreg_p","stopped",1)
    109       ENDIF
    110 
    111       IF( ifiltre== -2 .AND..NOT.griscal )     THEN
    112          PRINT *,' Cette routine ne calcule le filtre inverse que '
    113      &        , ' sur la grille des scalaires !'
    114          CALL abort_gcm("fitreg_p","stopped",1)
    115       ENDIF
    116 
    117       IF( ifiltre/=2 .AND.ifiltre/= - 2 )  THEN
    118          PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
    119      &        , ' corriger et repasser !'
    120          CALL abort_gcm("fitreg_p","stopped",1)
    121        ENDIF
    122 c
    123 
    124       iim2   = iim * iim
    125       immjm  = iim * jjm
    126 c
    127 c
    128       IF( griscal )   THEN
    129          IF( nlat. NE. jjp1 )  THEN
    130             CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
    131          ELSE
    132 c     
    133             IF( iaire==1 )  THEN
    134                sdd1_type = type_sddv
    135                sdd2_type = type_unsddv
    136             ELSE
    137                sdd1_type = type_unsddv
    138                sdd2_type = type_sddv
    139             ENDIF
    140 c
    141             jdfil1 = 2
    142             jffil1 = jfiltnu
    143             jdfil2 = jfiltsu
    144             jffil2 = jjm
    145          ENDIF
    146       ELSE
    147          IF( nlat/=jjm )  THEN
    148             CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
    149          ELSE
    150 c
    151             IF( iaire==1 )  THEN
    152                sdd1_type = type_sddu
    153                sdd2_type = type_unsddu
    154             ELSE
    155                sdd1_type = type_unsddu
    156                sdd2_type = type_sddu
    157             ENDIF
    158 c     
    159             jdfil1 = 1
    160             jffil1 = jfiltnv
    161             jdfil2 = jfiltsv
    162             jffil2 = jjm
    163          ENDIF
    164       ENDIF
    165 c     
    166       DO hemisph = 1, 2
    167 c     
    168          IF ( hemisph==1 )  THEN
    169 cym
    170             jdfil = max(jdfil1,ibeg)
    171             jffil = min(jffil1,iend)
    172          ELSE
    173 cym
    174             jdfil = max(jdfil2,ibeg)
    175             jffil = min(jffil2,iend)
    176          ENDIF
    177 
    178 
    179 cccccccccccccccccccccccccccccccccccccccccccc
    180 c Utilisation du filtre classique
    181 cccccccccccccccccccccccccccccccccccccccccccc
    182 
    183          IF (.NOT. use_filtre_fft) THEN
    184      
    185 c     !---------------------------------!
    186 c     ! Agregation des niveau verticaux !
    187 c     ! uniquement necessaire pour une  !
    188 c     ! execution OpenMP                !
    189 c     !---------------------------------!
    190             ll_nb = 0
    191 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    192             DO l = 1, nbniv
    193                ll_nb = ll_nb+1
    194                DO j = jdfil,jffil
    195                   DO i = 1, iim
    196                      champ_loc(i,j,ll_nb) =
    197      &                    champ(i,j,l) * sdd12(i,sdd1_type)
    198                   ENDDO
    199                ENDDO
    200             ENDDO
    201 c$OMP END DO NOWAIT
    202 
    203             nbniv_loc = ll_nb
    204 
    205             IF( hemisph==1 )      THEN
    206                
    207                IF( ifiltre==-2 )   THEN
    208                   DO j = jdfil,jffil
    209 #ifdef BLAS
    210                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    211      &                    matrinvn(1,1,j), iim,
    212      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    213      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
    214 #else
    215                      champ_fft(:iim,j-jdfil+1,:)
    216      &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
    217 #endif
    218                   ENDDO
    219                  
    220                ELSE IF ( griscal )     THEN
    221                   DO j = jdfil,jffil
    222 #ifdef BLAS
    223                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    224      &                    matriceun(1,1,j), iim,
    225      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    226      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
    227 #else
    228                      champ_fft(:iim,j-jdfil+1,:)
    229      &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
    230 #endif
    231                   ENDDO
    232                  
    233                ELSE
    234                   DO j = jdfil,jffil
    235 #ifdef BLAS
    236                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    237      &                    matricevn(1,1,j), iim,
    238      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    239      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
    240 #else
    241                      champ_fft(:iim,j-jdfil+1,:)
    242      &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
    243 #endif
    244                   ENDDO
    245                  
    246                ENDIF
    247                
    248             ELSE
    249                
    250                IF( ifiltre==-2 )   THEN
    251                   DO j = jdfil,jffil
    252 #ifdef BLAS
    253                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    254      &                    matrinvs(1,1,j-jfiltsu+1), iim,
    255      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    256      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
    257 #else
    258                      champ_fft(:iim,j-jdfil+1,:)
    259      &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
    260      &                            champ_loc(:iim,j,:))
    261 #endif
    262                   ENDDO
    263                  
    264                ELSE IF ( griscal )     THEN
    265                  
    266                   DO j = jdfil,jffil
    267 #ifdef BLAS
    268                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    269      &                    matriceus(1,1,j-jfiltsu+1), iim,
    270      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    271      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
    272 #else
    273                      champ_fft(:iim,j-jdfil+1,:)
    274      &                    =matmul(matriceus(:,:,j-jfiltsu+1),
    275      &                            champ_loc(:iim,j,:))
    276 #endif
    277                   ENDDO
    278                  
    279                ELSE
    280                  
    281                   DO j = jdfil,jffil
    282 #ifdef BLAS
    283                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    284      &                    matricevs(1,1,j-jfiltsv+1), iim,
    285      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    286      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
    287 #else
    288                      champ_fft(:iim,j-jdfil+1,:)
    289      &                    =matmul(matricevs(:,:,j-jfiltsv+1),
    290      &                            champ_loc(:iim,j,:))
    291 #endif
    292                   ENDDO
    293                  
    294                ENDIF
    295                
    296             ENDIF
    297 !     c     
    298             IF( ifiltre==2 )  THEN
    299                
    300 c     !-------------------------------------!
    301 c     ! Dés-agregation des niveau verticaux !
    302 c     ! uniquement necessaire pour une      !
    303 c     ! execution OpenMP                    !
    304 c     !-------------------------------------!
    305                ll_nb = 0
    306 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    307                DO l = 1, nbniv
    308                   ll_nb = ll_nb + 1
    309                   DO j = jdfil,jffil
    310                      DO i = 1, iim
    311                         champ( i,j,l ) = (champ_loc(i,j,ll_nb)
    312      &                       + champ_fft(i,j-jdfil+1,ll_nb))
    313      &                       * sdd12(i,sdd2_type)
    314                      ENDDO
    315                   ENDDO
    316                ENDDO
    317 c$OMP END DO NOWAIT
    318                
    319             ELSE
    320                
    321                ll_nb = 0
    322 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    323                DO l = 1, nbniv_loc
    324                   ll_nb = ll_nb + 1
    325                   DO j = jdfil,jffil
    326                      DO i = 1, iim
    327                         champ( i,j,l ) = (champ_loc(i,j,ll_nb)
    328      &                       - champ_fft(i,j-jdfil+1,ll_nb))
    329      &                       * sdd12(i,sdd2_type)
    330                      ENDDO
    331                   ENDDO
    332                ENDDO
    333 c$OMP END DO NOWAIT
    334                
    335             ENDIF
    336            
    337 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    338             DO l = 1, nbniv
    339                DO j = jdfil,jffil
    340                   champ( iip1,j,l ) = champ( 1,j,l )
    341                ENDDO
    342             ENDDO
    343 c$OMP END DO NOWAIT
    344            
    345 ccccccccccccccccccccccccccccccccccccccccccccc
    346 c Utilisation du filtre FFT
    347 ccccccccccccccccccccccccccccccccccccccccccccc
    348        
    349          ELSE
    350        
    351 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    352             DO l=1,nbniv
    353                DO j=jdfil,jffil
    354                   DO  i = 1, iim
    355                      champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
    356                      champ_fft( i,j,l) = champ(i,j,l)
    357                   ENDDO
    358                ENDDO
    359             ENDDO
    360 c$OMP END DO NOWAIT
    361 
    362             IF (jdfil<=jffil) THEN
    363                IF( ifiltre == -2 )   THEN
    364                   CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
    365                ELSE IF ( griscal )     THEN
    366                   CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
    367                ELSE
    368                   CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
    369                ENDIF
    370             ENDIF
    371 
    372 
    373             IF( ifiltre== 2 )  THEN
    374 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
    375                DO l=1,nbniv
    376                   DO j=jdfil,jffil
    377                      DO  i = 1, iim
    378                         champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
    379      &                       *sdd12(i,sdd2_type)
    380                      ENDDO
    381                   ENDDO
    382                ENDDO
    383 c$OMP END DO NOWAIT       
    384             ELSE
    385        
    386 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
    387                DO l=1,nbniv
    388                   DO j=jdfil,jffil
    389                      DO  i = 1, iim
    390                         champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
    391      &                       *sdd12(i,sdd2_type)
    392                      ENDDO
    393                   ENDDO
    394                ENDDO
    395 c$OMP END DO NOWAIT         
    396             ENDIF
    397 c
    398 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    399             DO l=1,nbniv
    400                DO j=jdfil,jffil
    401 !            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
    402                   champ( iip1,j,l ) = champ( 1,j,l )
    403                ENDDO
    404             ENDDO
    405 c$OMP END DO NOWAIT             
    406          ENDIF
    407 c Fin de la zone de filtrage
    408 
    409        
    410       ENDDO
    411 
    412 !      DO j=1,nlat
    413 
    414 !          PRINT *,"check FFT ----> Delta(",j,")=",
    415 !     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
    416 !     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
    417 !      ENDDO
    418      
    419 !          PRINT *,"check FFT ----> Delta(",j,")=",
    420 !     &            sum(champ-champ_fft)/sum(champ)
    421 
    422 c
    423  1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
    424      &     filtrer, sur la grille des scalaires'/)
    425  2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
    426      &     ltrer, sur la grille de V ou de Z'/)
    427 c$OMP MASTER     
    428       CALL stop_timer
    429 c$OMP END MASTER
    430       RETURN
    431       END
     3SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv, &
     4        ifiltre, iaire, griscal ,iter)
     5  USE parallel_lmdz, ONLY: OMP_CHUNK
     6  USE mod_filtre_fft
     7  USE timer_filtre
     8
     9  USE filtreg_mod
     10
     11  IMPLICIT NONE
     12
     13  !=======================================================================
     14  !
     15  !   Auteur: P. Le Van        07/10/97
     16  !   ------
     17  !
     18  !   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
     19  !                 pour l'operateur  Filtre    .
     20  !   ------
     21  !
     22  !   Arguments:
     23  !   ----------
     24  !
     25  !
     26  !  ibeg..iend            lattitude a filtrer
     27  !  nlat                  nombre de latitudes du champ
     28  !  nbniv                 nombre de niveaux verticaux a filtrer
     29  !  champ(iip1,nblat,nbniv)  en entree : champ a filtrer
     30  !                        en sortie : champ filtre
     31  !  ifiltre               +1  Transformee directe
     32  !                        -1  Transformee inverse
     33  !                        +2  Filtre directe
     34  !                        -2  Filtre inverse
     35  !
     36  !  iaire                 1   si champ intensif
     37  !                        2   si champ extensif (pondere par les aires)
     38  !
     39  !  iter                  1   filtre simple
     40  !
     41  !=======================================================================
     42  !
     43  !
     44  !                  Variable Intensive
     45  !            ifiltre = 1     filtre directe
     46  !            ifiltre =-1     filtre inverse
     47  !
     48  !                  Variable Extensive
     49  !            ifiltre = 2     filtre directe
     50  !            ifiltre =-2     filtre inverse
     51  !
     52  !
     53  INCLUDE "dimensions.h"
     54  INCLUDE "paramet.h"
     55  INCLUDE "coefils.h"
     56  !
     57  INTEGER :: ibeg,iend,nlat,nbniv,ifiltre,iter
     58  INTEGER :: i,j,l,k
     59  INTEGER :: iim2,immjm
     60  INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
     61
     62  REAL :: champ( iip1,nlat,nbniv)
     63
     64  LOGICAL :: griscal
     65  INTEGER :: hemisph, iaire
     66
     67  REAL :: champ_fft(iip1,nlat,nbniv)
     68  REAL :: champ_in(iip1,nlat,nbniv)
     69
     70  LOGICAL,SAVE     :: first=.TRUE.
     71!$OMP THREADPRIVATE(first)
     72
     73  REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
     74  INTEGER :: ll_nb, nbniv_loc
     75  REAL, SAVE :: sdd12(iim,4)
     76!$OMP THREADPRIVATE(sdd12)
     77
     78  INTEGER, PARAMETER :: type_sddu=1
     79  INTEGER, PARAMETER :: type_sddv=2
     80  INTEGER, PARAMETER :: type_unsddu=3
     81  INTEGER, PARAMETER :: type_unsddv=4
     82
     83  INTEGER :: sdd1_type, sdd2_type
     84
     85  IF (first) THEN
     86     sdd12(1:iim,type_sddu) = sddu(1:iim)
     87     sdd12(1:iim,type_sddv) = sddv(1:iim)
     88     sdd12(1:iim,type_unsddu) = unsddu(1:iim)
     89     sdd12(1:iim,type_unsddv) = unsddv(1:iim)
     90
     91     CALL Init_timer
     92     first=.FALSE.
     93  ENDIF
     94
     95!$OMP MASTER
     96  CALL start_timer
     97!$OMP END MASTER
     98
     99  !-------------------------------------------------------c
     100
     101  IF(ifiltre==1.or.ifiltre==-1) &
     102        CALL abort_gcm("fitreg_p","Pas de transformee simple&
     103        &dans cette version",1)
     104
     105  IF( iter== 2 )  THEN
     106     PRINT *,' Pas d iteration du filtre dans cette version !'&
     107           &        , ' Utiliser old_filtreg et repasser !'
     108     CALL abort_gcm("fitreg_p","stopped",1)
     109  ENDIF
     110
     111  IF( ifiltre== -2 .AND..NOT.griscal )     THEN
     112     PRINT *,' Cette routine ne calcule le filtre inverse que ' &
     113           , ' sur la grille des scalaires !'
     114     CALL abort_gcm("fitreg_p","stopped",1)
     115  ENDIF
     116
     117  IF( ifiltre/=2 .AND.ifiltre/= - 2 )  THEN
     118     PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
     119           , ' corriger et repasser !'
     120     CALL abort_gcm("fitreg_p","stopped",1)
     121   ENDIF
     122  !
     123
     124  iim2   = iim * iim
     125  immjm  = iim * jjm
     126  !
     127  !
     128  IF( griscal )   THEN
     129     IF( nlat /= jjp1 )  THEN
     130        CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
     131     ELSE
     132  !
     133        IF( iaire==1 )  THEN
     134           sdd1_type = type_sddv
     135           sdd2_type = type_unsddv
     136        ELSE
     137           sdd1_type = type_unsddv
     138           sdd2_type = type_sddv
     139        ENDIF
     140  !
     141        jdfil1 = 2
     142        jffil1 = jfiltnu
     143        jdfil2 = jfiltsu
     144        jffil2 = jjm
     145     ENDIF
     146  ELSE
     147     IF( nlat/=jjm )  THEN
     148        CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
     149     ELSE
     150  !
     151        IF( iaire==1 )  THEN
     152           sdd1_type = type_sddu
     153           sdd2_type = type_unsddu
     154        ELSE
     155           sdd1_type = type_unsddu
     156           sdd2_type = type_sddu
     157        ENDIF
     158  !
     159        jdfil1 = 1
     160        jffil1 = jfiltnv
     161        jdfil2 = jfiltsv
     162        jffil2 = jjm
     163     ENDIF
     164  ENDIF
     165  !
     166  DO hemisph = 1, 2
     167  !
     168     IF ( hemisph==1 )  THEN
     169  !ym
     170        jdfil = max(jdfil1,ibeg)
     171        jffil = min(jffil1,iend)
     172     ELSE
     173  !ym
     174        jdfil = max(jdfil2,ibeg)
     175        jffil = min(jffil2,iend)
     176     ENDIF
     177
     178
     179  !ccccccccccccccccccccccccccccccccccccccccccc
     180  ! Utilisation du filtre classique
     181  !ccccccccccccccccccccccccccccccccccccccccccc
     182
     183     IF (.NOT. use_filtre_fft) THEN
     184
     185  ! !---------------------------------!
     186  ! ! Agregation des niveau verticaux !
     187  ! ! uniquement necessaire pour une  !
     188  ! ! execution OpenMP                !
     189  ! !---------------------------------!
     190        ll_nb = 0
     191!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     192        DO l = 1, nbniv
     193           ll_nb = ll_nb+1
     194           DO j = jdfil,jffil
     195              DO i = 1, iim
     196                 champ_loc(i,j,ll_nb) = &
     197                       champ(i,j,l) * sdd12(i,sdd1_type)
     198              ENDDO
     199           ENDDO
     200        ENDDO
     201!$OMP END DO NOWAIT
     202
     203        nbniv_loc = ll_nb
     204
     205        IF( hemisph==1 )      THEN
     206
     207           IF( ifiltre==-2 )   THEN
     208              DO j = jdfil,jffil
     209#ifdef BLAS
     210                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
     211                       matrinvn(1,1,j), iim, &
     212                       champ_loc(1,j,1), iip1*nlat, 0.0, &
     213                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
     214#else
     215                 champ_fft(:iim,j-jdfil+1,:) &
     216                       =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
     217#endif
     218              ENDDO
     219
     220           ELSE IF ( griscal )     THEN
     221              DO j = jdfil,jffil
     222#ifdef BLAS
     223                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
     224                       matriceun(1,1,j), iim, &
     225                       champ_loc(1,j,1), iip1*nlat, 0.0, &
     226                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
     227#else
     228                 champ_fft(:iim,j-jdfil+1,:) &
     229                       =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
     230#endif
     231              ENDDO
     232
     233           ELSE
     234              DO j = jdfil,jffil
     235#ifdef BLAS
     236                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
     237                       matricevn(1,1,j), iim, &
     238                       champ_loc(1,j,1), iip1*nlat, 0.0, &
     239                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
     240#else
     241                 champ_fft(:iim,j-jdfil+1,:) &
     242                       =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
     243#endif
     244              ENDDO
     245
     246           ENDIF
     247
     248        ELSE
     249
     250           IF( ifiltre==-2 )   THEN
     251              DO j = jdfil,jffil
     252#ifdef BLAS
     253                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
     254                       matrinvs(1,1,j-jfiltsu+1), iim, &
     255                       champ_loc(1,j,1), iip1*nlat, 0.0, &
     256                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
     257#else
     258                 champ_fft(:iim,j-jdfil+1,:) &
     259                       =matmul(matrinvs(:,:,j-jfiltsu+1), &
     260                       champ_loc(:iim,j,:))
     261#endif
     262              ENDDO
     263
     264           ELSE IF ( griscal )     THEN
     265
     266              DO j = jdfil,jffil
     267#ifdef BLAS
     268                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
     269                       matriceus(1,1,j-jfiltsu+1), iim, &
     270                       champ_loc(1,j,1), iip1*nlat, 0.0, &
     271                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
     272#else
     273                 champ_fft(:iim,j-jdfil+1,:) &
     274                       =matmul(matriceus(:,:,j-jfiltsu+1), &
     275                       champ_loc(:iim,j,:))
     276#endif
     277              ENDDO
     278
     279           ELSE
     280
     281              DO j = jdfil,jffil
     282#ifdef BLAS
     283                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
     284                       matricevs(1,1,j-jfiltsv+1), iim, &
     285                       champ_loc(1,j,1), iip1*nlat, 0.0, &
     286                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
     287#else
     288                 champ_fft(:iim,j-jdfil+1,:) &
     289                       =matmul(matricevs(:,:,j-jfiltsv+1), &
     290                       champ_loc(:iim,j,:))
     291#endif
     292              ENDDO
     293
     294           ENDIF
     295
     296        ENDIF
     297  ! c
     298        IF( ifiltre==2 )  THEN
     299
     300  ! !-------------------------------------!
     301  ! ! Dés-agregation des niveau verticaux !
     302  ! ! uniquement necessaire pour une      !
     303  ! ! execution OpenMP                    !
     304  ! !-------------------------------------!
     305           ll_nb = 0
     306!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     307           DO l = 1, nbniv
     308              ll_nb = ll_nb + 1
     309              DO j = jdfil,jffil
     310                 DO i = 1, iim
     311                    champ( i,j,l ) = (champ_loc(i,j,ll_nb) &
     312                          + champ_fft(i,j-jdfil+1,ll_nb)) &
     313                          * sdd12(i,sdd2_type)
     314                 ENDDO
     315              ENDDO
     316           ENDDO
     317!$OMP END DO NOWAIT
     318
     319        ELSE
     320
     321           ll_nb = 0
     322!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     323           DO l = 1, nbniv_loc
     324              ll_nb = ll_nb + 1
     325              DO j = jdfil,jffil
     326                 DO i = 1, iim
     327                    champ( i,j,l ) = (champ_loc(i,j,ll_nb) &
     328                          - champ_fft(i,j-jdfil+1,ll_nb)) &
     329                          * sdd12(i,sdd2_type)
     330                 ENDDO
     331              ENDDO
     332           ENDDO
     333!$OMP END DO NOWAIT
     334
     335        ENDIF
     336
     337!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     338        DO l = 1, nbniv
     339           DO j = jdfil,jffil
     340              champ( iip1,j,l ) = champ( 1,j,l )
     341           ENDDO
     342        ENDDO
     343!$OMP END DO NOWAIT
     344
     345  !cccccccccccccccccccccccccccccccccccccccccccc
     346  ! Utilisation du filtre FFT
     347  !cccccccccccccccccccccccccccccccccccccccccccc
     348
     349     ELSE
     350
     351!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     352        DO l=1,nbniv
     353           DO j=jdfil,jffil
     354              DO  i = 1, iim
     355                 champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
     356                 champ_fft( i,j,l) = champ(i,j,l)
     357              ENDDO
     358           ENDDO
     359        ENDDO
     360!$OMP END DO NOWAIT
     361
     362        IF (jdfil<=jffil) THEN
     363           IF( ifiltre == -2 )   THEN
     364              CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     365           ELSE IF ( griscal )     THEN
     366              CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     367           ELSE
     368              CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     369           ENDIF
     370        ENDIF
     371
     372
     373        IF( ifiltre== 2 )  THEN
     374!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     375           DO l=1,nbniv
     376              DO j=jdfil,jffil
     377                 DO  i = 1, iim
     378                    champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l)) &
     379                          *sdd12(i,sdd2_type)
     380                 ENDDO
     381              ENDDO
     382           ENDDO
     383!$OMP END DO NOWAIT     
     384        ELSE
     385
     386!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     387           DO l=1,nbniv
     388              DO j=jdfil,jffil
     389                 DO  i = 1, iim
     390                    champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l)) &
     391                          *sdd12(i,sdd2_type)
     392                 ENDDO
     393              ENDDO
     394           ENDDO
     395!$OMP END DO NOWAIT
     396        ENDIF
     397  !
     398!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     399        DO l=1,nbniv
     400           DO j=jdfil,jffil
     401         ! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
     402              champ( iip1,j,l ) = champ( 1,j,l )
     403           ENDDO
     404        ENDDO
     405!$OMP END DO NOWAIT             
     406     ENDIF
     407  ! Fin de la zone de filtrage
     408
     409
     410  ENDDO
     411
     412   ! DO j=1,nlat
     413
     414   !     PRINT *,"check FFT ----> Delta(",j,")=",
     415  ! &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
     416  ! &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
     417  !  ENDDO
     418
     419  !      PRINT *,"check FFT ----> Delta(",j,")=",
     420  ! &            sum(champ-champ_fft)/sum(champ)
     421
     422  !
     423!$OMP MASTER
     424  CALL stop_timer
     425!$OMP END MASTER
     426
     427END SUBROUTINE filtreg_p
Note: See TracChangeset for help on using the changeset viewer.