Ignore:
Timestamp:
Jul 30, 2008, 5:50:03 PM (16 years ago)
Author:
Laurent Fairhead
Message:

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/trunk/libf/dyn3dpar/filtreg_p.F

    r763 r985  
     1
     2
    13      SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
    24     .                       ifiltre, iaire, griscal ,iter)
    35      USE Parallel, only : OMP_CHUNK
     6      USE mod_filtre_fft
     7      USE timer_filtre
    48      IMPLICIT NONE
    59
     
    5963     ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
    6064     ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
    61       REAL  eignq(iim), sdd1(iim),sdd2(iim)
     65cym      REAL  eignq(iim), sdd1(iim),sdd2(iim)
     66
     67      REAL  eignq(iim)
     68      REAL :: sdd1(iim),sdd2(iim)
     69     
    6270      LOGICAL    griscal
    6371      INTEGER    hemisph, iaire
    64 c
     72     
     73      REAL :: champ_fft(iip1,nlat,nbniv)
     74      REAL :: champ_in(iip1,nlat,nbniv)
     75     
     76      REAL,SAVE,TARGET :: sddu_loc(iim)
     77      REAL,SAVE,TARGET :: sddv_loc(iim)
     78      REAL,SAVE,TARGET :: unsddu_loc(iim)
     79      REAL,SAVE,TARGET :: unsddv_loc(iim)
     80c$OMP THREADPRIVATE(sddu_loc,sddv_loc,unsddu_loc,unsddv_loc)
     81      LOGICAL,SAVE     :: first=.TRUE.
     82c$OMP THREADPRIVATE(first)
     83
     84      IF (first) THEN
     85        sddu_loc(1:iim)=sddu(1:iim)
     86        sddv_loc(1:iim)=sddv(1:iim)
     87        unsddu_loc(1:iim)=unsddu(1:iim)
     88        unsddv_loc(1:iim)=unsddv(1:iim)
     89        CALL Init_timer
     90        first=.FALSE.
     91c       PRINT *,"----> sddu_loc=",sddu_loc
     92c       PRINT *,"----> sddv_loc=",sddv_loc
     93c       PRINT *,"----> unsddu_loc=",unsddu_loc
     94c       PRINT *,"----> unsddv_loc=",unsddv_loc
     95      ENDIF
     96
     97c$OMP MASTER     
     98      CALL start_timer
     99c$OMP END MASTER
    65100
    66101      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
     
    97132c
    98133             IF( iaire.EQ.1 )  THEN
    99                 CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
    100                 CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
     134cym                CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
     135cym                CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
     136cym               sdd1=>sddv_loc
     137cym               sdd2=>unsddv_loc
     138               sdd1(1:iim)=sddv_loc(1:iim)
     139               sdd2(1:iim)=unsddv_loc(1:iim)
    101140             ELSE
    102                 CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
    103                 CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
     141cym                CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
     142cym                CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
     143               sdd1(1:iim)=unsddv_loc(1:iim)
     144               sdd2(1:iim)=sddv_loc(1:iim)
    104145             END IF
    105146c
     
    116157c
    117158             IF( iaire.EQ.1 )  THEN
    118                 CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
    119                 CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
     159cym                CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
     160cym                CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
     161cym                sdd1=>sddu_loc
     162cym                sdd2=>unsddu_loc
     163                sdd1(1:iim)=sddu_loc(1:iim)
     164                sdd2(1:iim)=unsddu_loc(1:iim)
     165
    120166             ELSE
    121                 CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
    122                 CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
     167cym                CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
     168cym                CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
     169cym               sdd1=>unsddu_loc
     170cym               sdd2=>sddu_loc
     171               sdd1(1:iim)=unsddu_loc(1:iim)
     172               sdd2(1:iim)=sddu_loc(1:iim)
    123173             END IF
    124174c
     
    129179          END IF
    130180      END IF
     181
     182c      PRINT *,"APPEL a filtreg --> sdd1=",sdd1
     183c      PRINT *,"APPEL a filtreg --> sdd2=",sdd2
     184c      PRINT *,"----> sddu_loc=",sddu_loc
     185c       PRINT *,"----> sddv_loc=",sddv_loc
     186c       PRINT *,"----> unsddu_loc=",unsddu_loc
     187c       PRINT *,"----> unsddv_loc=",unsddv_loc
     188 
    131189c
    132190c
     
    143201      END IF
    144202
     203
     204cccccccccccccccccccccccccccccccccccccccccccc
     205c Utilisation du filtre classique
     206cccccccccccccccccccccccccccccccccccccccccccc
     207
     208      IF (.NOT. use_filtre_fft) THEN
     209     
    145210c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    146211      DO 50  l = 1, nbniv
    147       DO 30  j = jdfil,jffil
     212        DO 30  j = jdfil,jffil
    148213 
    149214 
    150       DO  5  i = 1, iim
    151       champ(i,j,l) = champ(i,j,l) * sdd1(i)
    152    5  CONTINUE
    153 c
    154 
    155       IF( hemisph. EQ. 1 )      THEN
    156 
    157         IF( ifiltre. EQ. -2 )   THEN
    158 #ifdef CRAY
    159          CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq  ,
    160      *                             1, iim, iim                         )
    161 #else
    162 #ifdef BLAS
    163       CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
    164      .           champ(1,j,l), 1, 0.0, eignq, 1)
    165 #else
    166       DO k = 1, iim
    167          eignq(k) = 0.0
    168       ENDDO
    169       DO k = 1, iim
    170       DO i = 1, iim
    171          eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
    172       ENDDO
    173       ENDDO
    174 #endif
    175 #endif
    176         ELSE IF ( griscal )     THEN
    177 #ifdef CRAY
    178          CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
    179      *                             1, iim, iim                         )
    180 #else
    181 #ifdef BLAS
    182       CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
    183      .           champ(1,j,l), 1, 0.0, eignq, 1)
    184 #else
    185       DO k = 1, iim
    186          eignq(k) = 0.0
    187       ENDDO
    188       DO i = 1, iim
    189       DO k = 1, iim
    190          eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
    191       ENDDO
    192       ENDDO
    193 #endif
    194 #endif
    195         ELSE
    196 #ifdef CRAY
    197          CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
    198      *                             1, iim, iim                         )
    199 #else
    200 #ifdef BLAS
    201       CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
    202      .           champ(1,j,l), 1, 0.0, eignq, 1)
    203 #else
    204       DO k = 1, iim
    205          eignq(k) = 0.0
    206       ENDDO
    207       DO i = 1, iim
    208       DO k = 1, iim
    209          eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
    210       ENDDO
    211       ENDDO
    212 #endif
    213 #endif
    214         ENDIF
    215 
    216       ELSE
    217 
    218         IF( ifiltre. EQ. -2 )   THEN
    219 #ifdef CRAY
    220          CALL MXVA( matrinvs(1,1,j-jfiltsu+1),  1, iim, champ(1,j,l),1 , 
    221      *                          eignq,  1, iim, iim                    )
    222 #else
    223 #ifdef BLAS
    224       CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
    225      .           champ(1,j,l), 1, 0.0, eignq, 1)
    226 #else
    227       DO k = 1, iim
    228          eignq(k) = 0.0
    229       ENDDO
    230       DO i = 1, iim
    231       DO k = 1, iim
    232          eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
    233       ENDDO
    234       ENDDO
    235 #endif
    236 #endif
    237         ELSE IF ( griscal )     THEN
    238 #ifdef CRAY
    239          CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
    240      *                          eignq,  1, iim, iim                    )
    241 #else
    242 #ifdef BLAS
    243       CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
    244      .           champ(1,j,l), 1, 0.0, eignq, 1)
    245 #else
    246       DO k = 1, iim
    247          eignq(k) = 0.0
    248       ENDDO
    249       DO i = 1, iim
    250       DO k = 1, iim
    251          eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
    252       ENDDO
    253       ENDDO
    254 #endif
    255 #endif
    256         ELSE
    257 #ifdef CRAY
    258          CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
    259      *                          eignq,  1, iim, iim                    )
    260 #else
    261 #ifdef BLAS
    262       CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
    263      .           champ(1,j,l), 1, 0.0, eignq, 1)
    264 #else
    265       DO k = 1, iim
    266          eignq(k) = 0.0
    267       ENDDO
    268       DO i = 1, iim
    269       DO k = 1, iim
    270          eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
    271       ENDDO
    272       ENDDO
    273 #endif
    274 #endif
    275         ENDIF
    276 
    277       ENDIF
    278 c
    279       IF( ifiltre.EQ. 2 )  THEN
    280         DO 15 i = 1, iim
    281         champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
    282   15    CONTINUE
    283       ELSE
    284         DO 16 i=1,iim
    285         champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
    286 16      CONTINUE
    287       ENDIF
    288 c
    289       champ( iip1,j,l ) = champ( 1,j,l )
    290 c
    291   30  CONTINUE
     215          DO  5  i = 1, iim
     216            champ(i,j,l) = champ(i,j,l) * sdd1(i)
     217   5      CONTINUE
     218c
     219
     220          IF( hemisph. EQ. 1 )      THEN
     221
     222            IF( ifiltre. EQ. -2 )   THEN
     223
     224
     225              CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
     226     .                     champ(1,j,l), 1, 0.0, eignq, 1)
     227
     228
     229            ELSE IF ( griscal )     THEN
     230
     231              CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
     232     .                    champ(1,j,l), 1, 0.0, eignq, 1)
     233
     234            ELSE
     235
     236              CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
     237     .                   champ(1,j,l), 1, 0.0, eignq, 1)
     238            ENDIF
     239
     240          ELSE
     241
     242            IF( ifiltre. EQ. -2 )   THEN
     243     
     244              CALL SGEMV("N",iim,iim,1.0, matrinvs(1,1,j-jfiltsu+1),iim,
     245     .                   champ(1,j,l), 1, 0.0, eignq, 1)
     246     
     247            ELSE IF ( griscal )     THEN
     248     
     249              CALL SGEMV("N",iim,iim,1.0,matriceus(1,1,j-jfiltsu+1),iim,
     250     .                   champ(1,j,l), 1, 0.0, eignq, 1)
     251            ELSE
     252         
     253              CALL SGEMV("N",iim,iim,1.0,matricevs(1,1,j-jfiltsv+1),iim,
     254     .                    champ(1,j,l), 1, 0.0, eignq, 1)
     255            ENDIF
     256
     257          ENDIF
     258
     259
     260c
     261          IF( ifiltre.EQ. 2 )  THEN
     262         
     263            DO 15 i = 1, iim
     264              champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
     265  15        CONTINUE
     266         
     267          ELSE
     268       
     269            DO 16 i=1,iim
     270               champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
     27116          CONTINUE
     272         
     273          ENDIF
     274c
     275          champ( iip1,j,l ) = champ( 1,j,l )
     276c
     277  30    CONTINUE
    292278c
    293279  50  CONTINUE
    294280c$OMP END DO NOWAIT
    295 c   
     281
     282ccccccccccccccccccccccccccccccccccccccccccccc
     283c Utilisation du filtre FFT
     284ccccccccccccccccccccccccccccccccccccccccccccc
     285       
     286       ELSE
     287       
     288c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     289          DO l=1,nbniv
     290            DO j=jdfil,jffil
     291              DO  i = 1, iim
     292                champ( i,j,l)= champ(i,j,l)*sdd1(i)
     293                champ_fft( i,j,l) = champ(i,j,l)
     294              ENDDO
     295            ENDDO
     296          ENDDO
     297c$OMP END DO NOWAIT
     298
     299      IF (jdfil<=jffil) THEN
     300        IF( ifiltre. EQ. -2 )   THEN
     301          CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     302        ELSE IF ( griscal )     THEN
     303          CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     304        ELSE
     305          CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
     306        ENDIF
     307      ENDIF
     308
     309
     310        IF( ifiltre.EQ. 2 )  THEN
     311c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
     312          DO l=1,nbniv
     313            DO j=jdfil,jffil
     314              DO  i = 1, iim
     315                champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
     316     .                             *sdd2(i)
     317              ENDDO
     318            ENDDO
     319          ENDDO
     320c$OMP END DO NOWAIT       
     321        ELSE
     322       
     323c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
     324          DO l=1,nbniv
     325            DO j=jdfil,jffil
     326              DO  i = 1, iim
     327                champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
     328     .                            *sdd2(i)
     329              ENDDO
     330            ENDDO
     331          ENDDO
     332c$OMP END DO NOWAIT         
     333        ENDIF
     334c
     335c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     336        DO l=1,nbniv
     337          DO j=jdfil,jffil
     338!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
     339            champ( iip1,j,l ) = champ( 1,j,l )
     340          ENDDO
     341        ENDDO
     342c$OMP END DO NOWAIT             
     343      ENDIF
     344c Fin de la zone de filtrage
     345
     346       
    296347 100  CONTINUE
     348
     349!      DO j=1,nlat
     350!     
     351!          PRINT *,"check FFT ----> Delta(",j,")=",
     352!     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
     353!     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
     354!      ENDDO
     355     
     356!          PRINT *,"check FFT ----> Delta(",j,")=",
     357!     &            sum(champ-champ_fft)/sum(champ)
     358!     
     359     
    297360c
    2983611111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
     
    3003632222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
    301364     *ltrer, sur la grille de V ou de Z'/)
     365c$OMP MASTER     
     366      CALL stop_timer
     367c$OMP END MASTER
    302368      RETURN
    303369      END
Note: See TracChangeset for help on using the changeset viewer.