Index: trunk/LMDZ.COMMON/libf/filtrez/filtreg.F =================================================================== --- trunk/LMDZ.COMMON/libf/filtrez/filtreg.F (revision 492) +++ trunk/LMDZ.COMMON/libf/filtrez/filtreg.F (revision 979) @@ -185,5 +185,5 @@ DO j = jdfil,jffil #ifdef BLAS - CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, + CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matrinvn(1,1,j), & iim, champ(1,j,1), iip1*nlat, 0.0, @@ -199,5 +199,5 @@ DO j = jdfil,jffil #ifdef BLAS - CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, + CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matriceun(1,1,j), & iim, champ(1,j,1), iip1*nlat, 0.0, @@ -213,5 +213,5 @@ DO j = jdfil,jffil #ifdef BLAS - CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, + CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matricevn(1,1,j), & iim, champ(1,j,1), iip1*nlat, 0.0, @@ -231,5 +231,5 @@ DO j = jdfil,jffil #ifdef BLAS - CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, + CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matrinvs(1,1,j-jfiltsu+1), & iim, champ(1,j,1), iip1*nlat, 0.0, @@ -247,5 +247,5 @@ DO j = jdfil,jffil #ifdef BLAS - CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, + CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matriceus(1,1,j-jfiltsu+1), & iim, champ(1,j,1), iip1*nlat, 0.0, @@ -262,5 +262,5 @@ DO j = jdfil,jffil #ifdef BLAS - CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, + CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & matricevs(1,1,j-jfiltsv+1), & iim, champ(1,j,1), iip1*nlat, 0.0, Index: trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90 (revision 492) +++ trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90 (revision 979) @@ -10,6 +10,6 @@ SUBROUTINE inifilr - USE mod_filtre_fft - ! + USE mod_filtre_fft, ONLY : use_filtre_fft,Init_filtre_fft + USE mod_filtre_fft_loc, ONLY : Init_filtre_fft_loc=>Init_filtre_fft ! ! ... H. Upadhyaya, O.Sharma ... ! @@ -541,4 +541,6 @@ CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu, & coefilv,modfrstv,jfiltnv,jfiltsv) + CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu, & + coefilv,modfrstv,jfiltnv,jfiltsv) ENDIF Index: trunk/LMDZ.COMMON/libf/filtrez/mod_fft.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/filtrez/mod_fft.F90 (revision 492) +++ trunk/LMDZ.COMMON/libf/filtrez/mod_fft.F90 (revision 979) @@ -3,11 +3,15 @@ #ifdef FFT_MATHKEISAN USE mod_fft_mathkeisan -#elif FFT_FFTW +#else +#ifdef FFT_FFTW USE mod_fft_fftw -#elif FFT_MKL +#else +#ifdef FFT_MKL USE mod_fft_mkl #else USE mod_fft_wrapper #endif +#endif +#endif END MODULE mod_fft Index: trunk/LMDZ.COMMON/libf/filtrez/mod_filtre_fft_loc.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/filtrez/mod_filtre_fft_loc.F90 (revision 979) +++ trunk/LMDZ.COMMON/libf/filtrez/mod_filtre_fft_loc.F90 (revision 979) @@ -0,0 +1,343 @@ +MODULE mod_filtre_fft_loc + + LOGICAL,SAVE :: use_filtre_fft + REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:) + REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:) + REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:) + +CONTAINS + + SUBROUTINE Init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv) + USE mod_fft + IMPLICIT NONE + include 'dimensions.h' + REAL, INTENT(IN) :: coeffu(iim,jjm) + INTEGER,INTENT(IN) :: modfrstu(jjm) + INTEGER,INTENT(IN) :: jfiltnu + INTEGER,INTENT(IN) :: jfiltsu + REAL, INTENT(IN) :: coeffv(iim,jjm) + INTEGER,INTENT(IN) :: modfrstv(jjm) + INTEGER,INTENT(IN) :: jfiltnv + INTEGER,INTENT(IN) :: jfiltsv + + INTEGER :: index_vp(iim) + INTEGER :: i,j + + index_vp(1)=1 + DO i=1,iim/2 + index_vp(i+1)=i*2 + ENDDO + + DO i=1,iim/2-1 + index_vp(iim/2+i+1)=iim-2*i+1 + ENDDO + + ALLOCATE(Filtre_u(iim,jjm)) + ALLOCATE(Filtre_v(iim,jjm)) + ALLOCATE(Filtre_inv(iim,jjm)) + + + DO j=2,jfiltnu + DO i=1,iim + IF (index_vp(i) < modfrstu(j)) THEN + Filtre_u(i,j)=0 + ELSE + Filtre_u(i,j)=coeffu(index_vp(i),j) + ENDIF + ENDDO + ENDDO + + DO j=jfiltsu,jjm + DO i=1,iim + IF (index_vp(i) < modfrstu(j)) THEN + Filtre_u(i,j)=0 + ELSE + Filtre_u(i,j)=coeffu(index_vp(i),j) + ENDIF + ENDDO + ENDDO + + DO j=1,jfiltnv + DO i=1,iim + IF (index_vp(i) < modfrstv(j)) THEN + Filtre_v(i,j)=0 + ELSE + Filtre_v(i,j)=coeffv(index_vp(i),j) + ENDIF + ENDDO + ENDDO + + DO j=jfiltsv,jjm + DO i=1,iim + IF (index_vp(i) < modfrstv(j)) THEN + Filtre_v(i,j)=0 + ELSE + Filtre_v(i,j)=coeffv(index_vp(i),j) + ENDIF + ENDDO + ENDDO + + DO j=2,jfiltnu + DO i=1,iim + IF (index_vp(i) < modfrstu(j)) THEN + Filtre_inv(i,j)=0 + ELSE + Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j)) + ENDIF + ENDDO + ENDDO + + DO j=jfiltsu,jjm + DO i=1,iim + IF (index_vp(i) < modfrstu(j)) THEN + Filtre_inv(i,j)=0 + ELSE + Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j)) + ENDIF + ENDDO + ENDDO + + +! CALL Init_FFT(iim,(jjm+1)*(llm+1)) + + + END SUBROUTINE Init_filtre_fft + + SUBROUTINE Filtre_u_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv) + USE mod_fft +#ifdef CPP_PARA + USE parallel,ONLY : OMP_CHUNK +#endif + IMPLICIT NONE + include 'dimensions.h' + INTEGER,INTENT(IN) :: jjb + INTEGER,INTENT(IN) :: jje + INTEGER,INTENT(IN) :: jj_begin + INTEGER,INTENT(IN) :: jj_end + INTEGER,INTENT(IN) :: nbniv + REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv) + + REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) + COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) + INTEGER :: nb_vect + INTEGER :: i,j,l + INTEGER :: ll_nb +! REAL :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv) + + ll_nb=0 +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,nbniv + ll_nb=ll_nb+1 + DO j=1,jj_end-jj_begin+1 + DO i=1,iim+1 + vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l) + ENDDO + ENDDO + ENDDO +!$OMP END DO NOWAIT + + nb_vect=(jj_end-jj_begin+1)*ll_nb + +! vect_tmp=vect + + CALL FFT_forward(vect,TF_vect,nb_vect) + +! CALL FFT_forward(vect,TF_vect_test,nb_vect) +! PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx" +! DO j=1,jj_end-jj_begin+1 +! DO i=1,iim/2+1 +! PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1) +! ENDDO +! ENDDO + + DO l=1,ll_nb + DO j=1,jj_end-jj_begin+1 + DO i=1,iim/2+1 + TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1) + ENDDO + ENDDO + ENDDO + + CALL FFT_backward(TF_vect,vect,nb_vect) +! CALL FFT_backward(TF_vect_test,vect_test,nb_vect) + +! PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx" +! DO j=1,jj_end-jj_begin+1 +! DO i=1,iim +! PRINT *,"====",i,j,"----->",vect_test(i,j,1) +! ENDDO +! ENDDO + + ll_nb=0 +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,nbniv + ll_nb=ll_nb+1 + DO j=1,jj_end-jj_begin+1 + DO i=1,iim+1 + vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb) + ENDDO + ENDDO + ENDDO +!$OMP END DO NOWAIT + + END SUBROUTINE Filtre_u_fft + + + SUBROUTINE Filtre_v_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv) + USE mod_fft +#ifdef CPP_PARA + USE parallel,ONLY : OMP_CHUNK +#endif + IMPLICIT NONE + INCLUDE 'dimensions.h' + INTEGER,INTENT(IN) :: jjb + INTEGER,INTENT(IN) :: jje + INTEGER,INTENT(IN) :: jj_begin + INTEGER,INTENT(IN) :: jj_end + INTEGER,INTENT(IN) :: nbniv + REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv) + + REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) + COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) + INTEGER :: nb_vect + INTEGER :: i,j,l + INTEGER :: ll_nb + + ll_nb=0 +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,nbniv + ll_nb=ll_nb+1 + DO j=1,jj_end-jj_begin+1 + DO i=1,iim+1 + vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l) + ENDDO + ENDDO + ENDDO +!$OMP END DO NOWAIT + + + nb_vect=(jj_end-jj_begin+1)*ll_nb + + CALL FFT_forward(vect,TF_vect,nb_vect) + + DO l=1,ll_nb + DO j=1,jj_end-jj_begin+1 + DO i=1,iim/2+1 + TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1) + ENDDO + ENDDO + ENDDO + + CALL FFT_backward(TF_vect,vect,nb_vect) + + + ll_nb=0 +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,nbniv + ll_nb=ll_nb+1 + DO j=1,jj_end-jj_begin+1 + DO i=1,iim+1 + vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb) + ENDDO + ENDDO + ENDDO +!$OMP END DO NOWAIT + + END SUBROUTINE Filtre_v_fft + + + SUBROUTINE Filtre_inv_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv) + USE mod_fft +#ifdef CPP_PARA + USE parallel,ONLY : OMP_CHUNK +#endif + IMPLICIT NONE + INCLUDE 'dimensions.h' + INTEGER,INTENT(IN) :: jjb + INTEGER,INTENT(IN) :: jje + INTEGER,INTENT(IN) :: jj_begin + INTEGER,INTENT(IN) :: jj_end + INTEGER,INTENT(IN) :: nbniv + REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv) + + REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv) + COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv) + INTEGER :: nb_vect + INTEGER :: i,j,l + INTEGER :: ll_nb + + ll_nb=0 +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,nbniv + ll_nb=ll_nb+1 + DO j=1,jj_end-jj_begin+1 + DO i=1,iim+1 + vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l) + ENDDO + ENDDO + ENDDO +!$OMP END DO NOWAIT + + nb_vect=(jj_end-jj_begin+1)*ll_nb + + CALL FFT_forward(vect,TF_vect,nb_vect) + + DO l=1,ll_nb + DO j=1,jj_end-jj_begin+1 + DO i=1,iim/2+1 + TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1) + ENDDO + ENDDO + ENDDO + + CALL FFT_backward(TF_vect,vect,nb_vect) + + ll_nb=0 +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,nbniv + ll_nb=ll_nb+1 + DO j=1,jj_end-jj_begin+1 + DO i=1,iim+1 + vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb) + ENDDO + ENDDO + ENDDO +!$OMP END DO NOWAIT + + END SUBROUTINE Filtre_inv_fft + + +! SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb) +! IMPLICIT NONE +! INTEGER,INTENT(IN) :: nbniv +! INTEGER,INTENT(OUT) :: ll_index(nbniv) +! INTEGER,INTENT(OUT) :: ll_nb +! +! INTEGER :: l,ll_begin, ll_end +! INTEGER :: omp_rank,omp_size +! INTEGER :: OMP_GET_NUM_THREADS +! INTEGER :: omp_chunk +! EXTERNAL OMP_GET_NUM_THREADS +! INTEGER :: OMP_GET_THREAD_NUM +! EXTERNAL OMP_GET_THREAD_NUM +! +! +! omp_size=OMP_GET_NUM_THREADS() +! omp_rank=OMP_GET_THREAD_NUM() +! omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size)) +! +! ll_begin=omp_rank*OMP_CHUNK+1 +! ll_nb=0 +! DO WHILE (ll_begin<=nbniv) +! ll_end=min(ll_begin+OMP_CHUNK-1,nbniv) +! DO l=ll_begin,ll_end +! ll_nb=ll_nb+1 +! ll_index(ll_nb)=l +! ENDDO +! ll_begin=ll_begin+omp_size*OMP_CHUNK +! ENDDO +! +! END SUBROUTINE get_ll_index + +END MODULE mod_filtre_fft_loc +