#include "CPP_OPTIONS.h"
#include "GMREDI_OPTIONS.h"

      subroutine smooth_init3D (mythid)


      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "GAD.h"
c#include "tamc.h"
#include "FFIELDS.h"
#include "EOS.h"
#include "GMREDI.h"
#include "smooth.h"


c     input
c     bi, bj : array indices
c     k      : vertical index used for masking
      integer i,j,k, bi, bj, imin, imax, jmin, jmax
      integer itlo,ithi
      integer jtlo,jthi
      integer myThid
      character*( 80) fnamegeneric

      _RL wc01theta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL wc01salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
c parameter to restrain the Kz based on grid cells
      _RL wc01_3D_KzMax
c to rotate the diffusion
      _RL Kuxprime, Kvyprime, rotate_s2,rotate_cos,rotate_sin
      _RL rotaTmp1,rotaTmp2,rotaTmp3
      integer ii,jj,kk
# ifndef GM_EXTRA_DIAGONAL
      _RL Kuz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL Kvz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
# endif

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

      wc01_dt=1.
      wc01_T=wc01_nbt(smoothOpNbCur)*wc01_dt

            WRITE(standardMessageUnit,'(A,2I4,/,3f5.2)')
     & 'smooth 3D default parameters: ',
     & wc01_nbt(smoothOpNbCur),wc01_T,
     & wc01_3D_Lx0(smoothOpNbCur),wc01_3D_Ly0(smoothOpNbCur),
     & wc01_3D_Lz0(smoothOpNbCur)

cgf "pour rotation H: sauver nouveaux champs"
cgf "et poser une limite [deep interior -> isotropic]"
cgf "... eviter les sauts de direction"

c here fill the wc01_3D_Lz array
      if ((smooth3DtypeZ(smoothOpNbCur).NE.0).AND.
     & (smooth3DsizeZ(smoothOpNbCur).EQ.3)) then
      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3DscalesZ',smoothOpNbCur
      call mdsreadfield(fnamegeneric,64,'RL',nR,
     & wc01_3D_Lz,1,mythid)
      _EXCH_XYZ_RL( wc01_3D_Lz, mythid )
      else
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
          wc01_3D_Lz(i,j,k,bi,bj)=wc01_3D_Lz0(smoothOpNbCur)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      endif


c vertical diffusion
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
       kappaRwc01(i,j,k,bi,bj)=wc01_3D_Lz(i,j,k,bi,bj)*
     & wc01_3D_Lz(i,j,k,bi,bj)/wc01_T/2
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c begin: to restrain the Kz based on grid cells
      if (smooth3DsizeZ(smoothOpNbCur).NE.3) then
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr 
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx

       wc01_3D_KzMax=drC(k)
       wc01_3D_KzMax=wc01_3D_KzMax*wc01_3D_KzMax/wc01_T/2
       if (kappaRwc01(i,j,k,bi,bj).GT.wc01_3D_KzMax) then
       kappaRwc01(i,j,k,bi,bj)=wc01_3D_KzMax 
       endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      endif
c end: to restrain the Kz based on grid cells

      _EXCH_XYZ_RL( kappaRwc01,    myThid )

c horizontal/isopycnal operator:

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
       wc01_Kuy(i,j,k,bi,bj)=0.
       wc01_Kvx(i,j,k,bi,bj)=0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO


      if ((smooth3DtypeH(smoothOpNbCur).EQ.2).OR.
     & (smooth3DtypeH(smoothOpNbCur).EQ.3)) then

c isopycnal operator:

      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3DscalesH',smoothOpNbCur

      call mdsreadfield(fnamegeneric,64,'RL',nR,
     & wc01theta,1,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,
     & wc01salt,2,mythid)
      _EXCH_XYZ_RL( wc01theta, mythid )
      _EXCH_XYZ_RL( wc01salt, mythid )

      if (smooth3DsizeH(smoothOpNbCur).EQ.3) then
      call mdsreadfield(fnamegeneric,64,'RL',nR,
     & wc01_3D_Lx,3,mythid)
      _EXCH_XYZ_RL( wc01_3D_Lx, mythid )
      else
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
          wc01_3D_Lx(i,j,k,bi,bj)=wc01_3D_Lx0(smoothOpNbCur)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      endif

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy 
          DO i=1-OLx,sNx+OLx
c here wc01_3D_Lx contains K divided by Kgmredi(=1000) 
      wc01_3D_Lx(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*
     & wc01_3D_Lx(i,j,k,bi,bj)/wc01_T/2 /1000
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XYZ_RL( wc01_3D_Lx, mythid )

        iMin = 1-OLx
        iMax = sNx+OLx
        jMin = 1-OLy
        jMax = sNy+OLy

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=Nr,1,-1
            CALL FIND_RHO_2D(
     I        iMin, iMax, jMin, jMax, k,
     I        wc01theta(1-OLx,1-OLy,k,bi,bj),
     I        wc01salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoK(1-OLx,1-OLy,bi,bj),
     I        k, bi, bj, myThid ) 
            IF (k.GT.1) THEN
             CALL FIND_RHO_2D(
     I        iMin, iMax, jMin, jMax, k,
     I        wc01theta(1-OLx,1-OLy,k-1,bi,bj),
     I        wc01salt(1-OLx,1-OLy,k-1,bi,bj),
     O        rhoKm1(1-OLx,1-OLy,bi,bj),
     I        k-1, bi, bj, myThid ) 
            ELSE
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
              rhoKm1(i,j,bi,bj)=rhoK(i,j,bi,bj)
          ENDDO
         ENDDO
            ENDIF
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
              rhoKp1(i,j,bi,bj)=rhoK(i,j,bi,bj)
          ENDDO
         ENDDO
cgf rk: GRAD_SIGMA ne calcule la derivee verticale au point w, entre Km1 et K
            CALL GRAD_SIGMA(
     &             bi, bj, iMin, iMax, jMin, jMax, k,
     &             rhoK(1-OLx,1-OLy,bi,bj), 
     &             rhoKm1(1-OLx,1-OLy,bi,bj), 
     &             rhoKp1(1-OLx,1-OLy,bi,bj),
     &             sigmaX(1-OLx,1-OLy,1,bi,bj), 
     &             sigmaY(1-OLx,1-OLy,1,bi,bj), 
     &             sigmaR(1-OLx,1-OLy,1,bi,bj),
     I             myThid )
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XYZ_RL( sigmaX, myThid )
      _EXCH_XYZ_RL( sigmaY, myThid )
      _EXCH_XYZ_RL( sigmaR, myThid )

      DO bj=jtlo,jthi
       DO bi=itlo,ithi

          CALL GMREDI_CALC_TENSOR(
     I             bi, bj, iMin, iMax, jMin, jMax,
     I             sigmaX(1-OLx,1-OLy,1,bi,bj), 
     &             sigmaY(1-OLx,1-OLy,1,bi,bj), 
     &             sigmaR(1-OLx,1-OLy,1,bi,bj),
     I             myThid )

        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
      if (smooth3DtypeH(smoothOpNbCur).EQ.2) then
       Kwx(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kwx(i,j,k,bi,bj)
       Kwy(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kwy(i,j,k,bi,bj)
# ifdef GM_EXTRA_DIAGONAL
       Kuz(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kuz(i,j,k,bi,bj)
       Kvz(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kvz(i,j,k,bi,bj)
# else
       Kuz(i,j,k,bi,bj)=0.
       Kvz(i,j,k,bi,bj)=0.
# endif
      else
       Kwx(i,j,k,bi,bj)=2*wc01_3D_Lx(i,j,k,bi,bj)*Kwx(i,j,k,bi,bj)
       Kwy(i,j,k,bi,bj)=2*wc01_3D_Lx(i,j,k,bi,bj)*Kwy(i,j,k,bi,bj)
       Kuz(i,j,k,bi,bj)=0.
       Kvz(i,j,k,bi,bj)=0.
      endif 
       Kwz(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kwz(i,j,k,bi,bj)
       Kux(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kux(i,j,k,bi,bj)
       Kvy(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*Kvy(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO

c begin: to restrain the Kz based on grid cells
        DO k=1,Nr 
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx

       wc01_3D_KzMax=drC(k)
       wc01_3D_KzMax=wc01_3D_KzMax*wc01_3D_KzMax/wc01_T/2
         
       if (Kwz(i,j,k,bi,bj).GT.wc01_3D_KzMax) then
       Kwx(i,j,k,bi,bj)=Kwx(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       Kwy(i,j,k,bi,bj)=Kwy(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       Kuz(i,j,k,bi,bj)=Kuz(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       Kvz(i,j,k,bi,bj)=Kvz(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       Kux(i,j,k,bi,bj)=Kux(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       Kvy(i,j,k,bi,bj)=Kvy(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       Kwz(i,j,k,bi,bj)=Kwz(i,j,k,bi,bj)
     & *wc01_3D_KzMax/Kwz(i,j,k,bi,bj)
       endif
          ENDDO
         ENDDO
        ENDDO
c end: to restrain the Kz based on grid cells


          CALL GMREDI_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRwc01(1-OLx,1-OLy,1,bi,bj),
     I        myThid)

       ENDDO
      ENDDO

      else

c hoizontal operator:


      if ((smooth3DtypeH(smoothOpNbCur).NE.0).AND.
     & (smooth3DsizeH(smoothOpNbCur).EQ.3)) then
      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3DscalesH',smoothOpNbCur
      call mdsreadfield(fnamegeneric,64,'RL',nR,
     & wc01_3D_Lx,1,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',nR,
     & wc01_3D_Ly,2,mythid)
      _EXCH_XYZ_RL( wc01_3D_Lx, mythid )
      _EXCH_XYZ_RL( wc01_3D_Ly, mythid )
      else
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
          wc01_3D_Lx(i,j,k,bi,bj)=wc01_3D_Lx0(smoothOpNbCur)
          wc01_3D_Ly(i,j,k,bi,bj)=wc01_3D_Ly0(smoothOpNbCur)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      endif

      if (smooth3DtypeH(smoothOpNbCur).NE.4) then  
c along model axes
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
      Kwx(i,j,k,bi,bj)=0.
      Kwy(i,j,k,bi,bj)=0.
      Kwz(i,j,k,bi,bj)=0.
      Kux(i,j,k,bi,bj)=wc01_3D_Lx(i,j,k,bi,bj)*
     & wc01_3D_Lx(i,j,k,bi,bj)/wc01_T/2
      Kvy(i,j,k,bi,bj)=wc01_3D_Ly(i,j,k,bi,bj)*
     & wc01_3D_Ly(i,j,k,bi,bj)/wc01_T/2
      Kuz(i,j,k,bi,bj)=0.
      Kvz(i,j,k,bi,bj)=0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      else

c along rotated axes

      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3DscalesH',smoothOpNbCur
      if (smooth3DsizeH(smoothOpNbCur).EQ.3) then
      	call mdsreadfield(fnamegeneric,64,'RL',nR,
     & 	wc01theta,3,mythid)
      else
      	call mdsreadfield(fnamegeneric,64,'RL',nR,
     & 	wc01theta,1,mythid)
      endif
      _EXCH_XYZ_RL( wc01theta, mythid )

        iMin = 1-OLx
        iMax = sNx+OLx
        jMin = 1-OLy
        jMax = sNy+OLy 

      write(fnamegeneric(1:80),'(1a)') 'wc01_3Dtest'
      
c compute the gradients from the "direction" field
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=Nr,1,-1
            CALL GRAD_SIGMA(
     &             bi, bj, iMin, iMax, jMin, jMax, k,
     &             wc01theta(1-OLx,1-OLy,k,bi,bj),
     &             wc01theta(1-OLx,1-OLy,k,bi,bj),
     &             wc01theta(1-OLx,1-OLy,k,bi,bj),
     &             sigmaX(1-OLx,1-OLy,1,bi,bj),
     &             sigmaY(1-OLx,1-OLy,1,bi,bj),
     &             sigmaR(1-OLx,1-OLy,1,bi,bj),
     I             myThid )
        ENDDO
       ENDDO 
      ENDDO
      _EXCH_XYZ_RL( sigmaX, myThid )
      _EXCH_XYZ_RL( sigmaY, myThid )
      _EXCH_XYZ_RL( sigmaR, myThid )

      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,sigmaX,1,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,sigmaY,2,1,mythid)

c available for the following computation:
c rhok,rhokm1,rhokp1,wc01salt,sigmar

c compute the associated cos and sin
c	rk1: Kwx is cos // Kwy is sin
c	rk2: 2 is used as a missing value
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
      rotate_s2=sigmaX(i,j,k,bi,bj)*sigmaX(i,j,k,bi,bj)
     & +sigmaY(i,j,k,bi,bj)*sigmaY(i,j,k,bi,bj)
      if ((rotate_s2.GT.0.).AND.(_maskS(i,j,k,bi,bj).NE.0.)
     & .AND.(_maskW(i,j,k,bi,bj).NE.0.) ) then
      Kwx(i,j,k,bi,bj)=sigmaY(i,j,k,bi,bj)/sqrt(rotate_s2)
      Kwy(i,j,k,bi,bj)=-sigmaX(i,j,k,bi,bj)/sqrt(rotate_s2)
      else
      Kwx(i,j,k,bi,bj)=2.
      Kwy(i,j,k,bi,bj)=2.
      endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,kwx,3,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,kwy,4,1,mythid)

c compute a saturation coefficient: where the angle is changing (heterogeneous)
c	we will stay isotropic
c	rk1: Kwz is the angle // wc01salt is the saturation coeff
c 	rk2: the computation uses atan to compute the angle, and has 
c	to be done twice because atan is not continuous at pi/2 
      DO kk=1,2

c initialization
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
       Kwz(i,j,k,bi,bj)=999.
      if ((Kwx(i,j,k,bi,bj).NE.2.).AND.
     & (Kwx(i,j,k,bi,bj).NE.0.)) then 
       		Kwz(i,j,k,bi,bj)=atan(Kwy(i,j,k,bi,bj)
     & 		/Kwx(i,j,k,bi,bj))
      elseif (Kwx(i,j,k,bi,bj).NE.2.) then
       		Kwz(i,j,k,bi,bj)=sign(pi/2.,Kwy(i,j,k,bi,bj))  
      endif
      if (kk.EQ.1) then
      		wc01salt(i,j,k,bi,bj)=999. 
      endif
c 	rk: it is important that the missing value is a (large) positive value
      if ((kk.EQ.2).AND.(Kwz(i,j,k,bi,bj).LT.0.)) then
      		Kwz(i,j,k,bi,bj)=Kwz(i,j,k,bi,bj)+pi 
      endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c      _EXCH_XYZ_RL( Kwz,    myThid )
c      _EXCH_XYZ_RL( wc01salt,    myThid )

c the computation/update of the saturation coeff
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1,sNy
          DO i=1,sNx
       if (Kwz(i,j,k,bi,bj).NE.999.) then
       	rotaTmp1=0.
       	rotaTmp2=0.
       	rotaTmp3=0.
       	do ii=-1,1
       	do jj=-1,1
       		if (Kwz(i+ii,j+jj,k,bi,bj).NE.999.) then
       		rotaTmp1=rotaTmp1+Kwz(i+ii,j+jj,k,bi,bj)
       		rotaTmp2=rotaTmp2+Kwz(i+ii,j+jj,k,bi,bj)*Kwz(i+ii,j+jj,k,bi,bj)
       		rotaTmp3=rotaTmp3+1.
       		endif
       	enddo
       	enddo
        rotaTmp1=rotaTmp1/rotaTmp3
        rotaTmp2=rotaTmp2/rotaTmp3
        rotaTmp3=rotaTmp2-rotaTmp1*rotaTmp1
        wc01salt(i,j,k,bi,bj)=min(wc01salt(i,j,k,bi,bj),rotaTmp3)
      		if (kk.EQ.2) then
        rotaTmp3= (1 _d +00 - wc01salt(i,j,k,bi,bj)/(pi/2/6))  
        wc01salt(i,j,k,bi,bj)=max(0 _d +00 , rotaTmp3)
      		endif
       endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_RL( wc01salt,    myThid )
      ENDDO !      DO kk=1,2

      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,wc01salt,5,1,mythid)

c finally, compute the diffusion operator 
c       rk: I will need to double-check the limit case (boundary)
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
      if (Kwz(i,j,k,bi,bj).NE.999.) then

      	rotaTmp1=wc01_3D_Lx(i,j,k,bi,bj)*
     & 	wc01_3D_Lx(i,j,k,bi,bj)/wc01_T/2
      	rotaTmp2=wc01_3D_Ly(i,j,k,bi,bj)*
     & 	wc01_3D_Ly(i,j,k,bi,bj)/wc01_T/2

      	Kuxprime=rotaTmp1
      	Kvyprime=rotaTmp2*wc01salt(i,j,k,bi,bj)
     & 	+ rotaTmp1*(1.-wc01salt(i,j,k,bi,bj))

      	Kux(i,j,k,bi,bj)=Kwx(i,j,k,bi,bj)*Kwx(i,j,k,bi,bj)*Kuxprime
     & 	+Kwy(i,j,k,bi,bj)*Kwy(i,j,k,bi,bj)*Kvyprime
      	wc01_Kuy(i,j,k,bi,bj)=Kwx(i,j,k,bi,bj)*Kwy(i,j,k,bi,bj)
     & 	*(-Kuxprime+Kvyprime)
      	Kvy(i,j,k,bi,bj)=Kwy(i,j,k,bi,bj)*Kwy(i,j,k,bi,bj)*Kuxprime
     & 	+Kwx(i,j,k,bi,bj)*Kwx(i,j,k,bi,bj)*Kvyprime
      	wc01_Kvx(i,j,k,bi,bj)=Kwx(i,j,k,bi,bj)*Kwy(i,j,k,bi,bj)
     & 	*(-Kuxprime+Kvyprime)

      else

        rotaTmp1=wc01_3D_Lx(i,j,k,bi,bj)*
     &  wc01_3D_Lx(i,j,k,bi,bj)/wc01_T/2
      Kux(i,j,k,bi,bj)=rotaTmp1
      wc01_Kuy(i,j,k,bi,bj)=0.
      Kvy(i,j,k,bi,bj)=rotaTmp1
      wc01_Kvx(i,j,k,bi,bj)=0.      

      endif

      Kwx(i,j,k,bi,bj)=0.
      Kwy(i,j,k,bi,bj)=0.
      Kwz(i,j,k,bi,bj)=0.
      Kuz(i,j,k,bi,bj)=0.
      Kvz(i,j,k,bi,bj)=0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c OLD VERSION
c      Kuxprime=wc01_3D_Lx(i,j,k,bi,bj)*
c     & wc01_3D_Lx(i,j,k,bi,bj)/wc01_T/2
c      Kvyprime=wc01_3D_Ly(i,j,k,bi,bj)*
c     & wc01_3D_Ly(i,j,k,bi,bj)/wc01_T/2

c      rotate_cos=0.7071
c      rotate_sin=0.7071
cc      rotate_cos=0.
cc      rotate_sin=1.

c      Kux(i,j,k,bi,bj)=rotate_cos*Kuxprime
c     & -rotate_sin*Kvyprime
c      Kvy(i,j,k,bi,bj)=rotate_sin*Kuxprime
c     & +rotate_cos*Kvyprime

c      DO bj=jtlo,jthi
c       DO bi=itlo,ithi
c        DO k=1,Nr
c         DO j=1-OLy,sNy+OLy
c          DO i=1-OLx,sNx+OLx
c      Kux(i,j,k,bi,bj)=rotate_cos*rotate_cos*Kuxprime
c     & +rotate_sin*rotate_sin*Kvyprime
c      wc01_Kuy(i,j,k,bi,bj)=rotate_cos*rotate_sin
c     & *(-Kuxprime+Kvyprime)
c      Kvy(i,j,k,bi,bj)=rotate_sin*rotate_sin*Kuxprime
c     & +rotate_cos*rotate_cos*Kvyprime
c      wc01_Kvx(i,j,k,bi,bj)=rotate_cos*rotate_sin
c     & *(-Kuxprime+Kvyprime)
c      Kwx(i,j,k,bi,bj)=0.
c      Kwy(i,j,k,bi,bj)=0.
c      Kwz(i,j,k,bi,bj)=0.
c      Kuz(i,j,k,bi,bj)=0.
c      Kvz(i,j,k,bi,bj)=0.
c          ENDDO
c         ENDDO
c        ENDDO
c       ENDDO
c      ENDDO

      endif

      endif


c finalize the set sup: 

      _EXCH_XYZ_RL( kappaRwc01,    myThid )

      _EXCH_XYZ_RL( Kwx,    myThid )
      _EXCH_XYZ_RL( Kwy,    myThid )
      _EXCH_XYZ_RL( Kwz,    myThid )
      _EXCH_XYZ_RL( Kux,    myThid )
      _EXCH_XYZ_RL( Kvy,    myThid )
      _EXCH_XYZ_RL( Kuz,    myThid )
      _EXCH_XYZ_RL( Kvz,    myThid )

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
       wc01_Kwx(i,j,k,bi,bj)=Kwx(i,j,k,bi,bj)
       wc01_Kwy(i,j,k,bi,bj)=Kwy(i,j,k,bi,bj)
       wc01_Kwz(i,j,k,bi,bj)=Kwz(i,j,k,bi,bj)
       wc01_Kux(i,j,k,bi,bj)=Kux(i,j,k,bi,bj)
       wc01_Kvy(i,j,k,bi,bj)=Kvy(i,j,k,bi,bj)
       wc01_Kuz(i,j,k,bi,bj)=Kuz(i,j,k,bi,bj)
       wc01_Kvz(i,j,k,bi,bj)=Kvz(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c write the diffusion operator into file:

      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3Doperator',smoothOpNbCur

      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kwx,1,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kwy,2,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kwz,3,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kux,4,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kvy,5,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kuz,6,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,Kvz,7,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,wc01_Kuy,8,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,wc01_Kvx,9,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & nR,kappaRwc01,10,1,mythid)


      END