Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (22 hours ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File:
1 moved

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/dyn3dmem/mod_filtreg_p.F90

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