#include "CPP_OPTIONS.h"

      subroutine smooth_correl3D (
     U     fld_in,smoothOpNb,mythid)


      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GRID.h"
#include "PARAMS.h"
c#include "tamc.h"
#include "smooth.h"


      _RL fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
      integer smoothOpNb
      integer nbt_in
      character*( 80) fnamegeneric
      integer i,j,k,bi,bj
      integer itlo,ithi
      integer jtlo,jthi
      integer myThid


      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)

      
      smoothOpNbCur=smoothOpNb
      nbt_in=wc01_nbt(smoothOpNbCur)/2

cgf tmp: set the diffkh to 0
cgf later: have an index for the operator number
cfg	and do not forget common block parts (nbt...)

c      write(fnamegeneric(1:80),'(1a)')
c     &    'wc01_3Doperator'
      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3Doperator',smoothOpNbCur
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kwx,
     & 1,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kwy,
     & 2,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kwz,
     & 3,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kux,
     & 4,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kvy,
     & 5,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kuz,
     & 6,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kvz,
     & 7,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kuy,
     & 8,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,wc01_Kvx,
     & 9,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,kappaRwc01,
     & 10,mythid)
      _EXCH_XYZ_RL ( wc01_Kwx, myThid )
      _EXCH_XYZ_RL ( wc01_Kwy, myThid )
      _EXCH_XYZ_RL ( wc01_Kwz, myThid )
      _EXCH_XYZ_RL ( wc01_Kux, myThid )
      _EXCH_XYZ_RL ( wc01_Kvy, myThid )
      _EXCH_XYZ_RL ( wc01_Kuz, myThid )
      _EXCH_XYZ_RL ( wc01_Kvz, myThid )
      _EXCH_XYZ_RL ( wc01_Kuy, myThid )
      _EXCH_XYZ_RL ( wc01_Kvx, myThid )
      _EXCH_XYZ_RL ( kappaRwc01, myThid )

c      write(fnamegeneric(1:80),'(1a)')
c     &    'wc01_filternorm'
      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3Dfilternorm',smoothOpNbCur
      call mdsreadfield(fnamegeneric,64,'RL',nR,normFilterwc01,
     & 1,mythid)
      _EXCH_XYZ_RL ( normFilterwc01, myThid )


      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1,sNy
          DO i=1,sNx 
c division by ~sqrt(volume):
           fld_in(i,j,k,bi,bj)=fld_in(i,j,k,bi,bj)
     & *sqrt(recip_rA(i,j,bi,bj)*recip_drF(k))
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XYZ_RL ( fld_in , myThid )

      call smooth_diff3D(fld_in,nbt_in,mythid)


       do bj = jtlo,jthi
        do bi = itlo,ithi
         DO j = 1,sNy
          DO i = 1,sNx
           DO k = 1,nR
c division by ~sqrt(var(filter)):
       fld_in(i,j,k,bi,bj)=fld_in(i,j,k,bi,bj)
     & *normFilterwc01(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO

      _EXCH_XYZ_RL ( fld_in , myThid )

      end