#include "CPP_OPTIONS.h"

      subroutine smooth_init2D ( mythid )

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

      integer i,j,k, bi, bj
      integer itlo,ithi
      integer jtlo,jthi
      integer myThid
      character*( 80) fnamegeneric

c wc01_norm : 
c   1) in the 2D isotropic case without boundaries, wc01_norm is known
c   analytically ~ sqrt(2*3.1415*wc01_L*wc01_L)
c   2) in general case: need to compute the filter variance

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


      wc01_2D_dt=1.
      wc01_2D_T=wc01_2D_nbt(smoothOpNbCur)*wc01_2D_dt

      if ((smooth2Dtype(smoothOpNbCur).NE.0).AND.
     & (smooth2Dsize(smoothOpNbCur).EQ.2)) then
      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_2Dscales',smoothOpNbCur
      call mdsreadfield(fnamegeneric,64,'RL',1,
     & wc01_2D_Lx,1,mythid)
      call mdsreadfield(fnamegeneric,64,'RL',1,
     & wc01_2D_Ly,2,mythid)
      _EXCH_XY_R8 ( wc01_2D_Lx, myThid )
      _EXCH_XY_R8 ( wc01_2D_Ly, myThid )
      else
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
          wc01_2D_Lx(i,j,bi,bj)=wc01_2D_Lx0(smoothOpNbCur)
          wc01_2D_Ly(i,j,bi,bj)=wc01_2D_Ly0(smoothOpNbCur)
          ENDDO
         ENDDO
       ENDDO
      ENDDO
      endif

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
          diffKhwc01_2D_x(i,j,bi,bj)=wc01_2D_Lx(i,j,bi,bj)*
     & wc01_2D_Lx(i,j,bi,bj)/wc01_2D_T/2
          diffKhwc01_2D_y(i,j,bi,bj)=wc01_2D_Ly(i,j,bi,bj)*
     & wc01_2D_Ly(i,j,bi,bj)/wc01_2D_T/2
          ENDDO
         ENDDO
       ENDDO
      ENDDO

      _EXCH_XY_R8 ( diffKhwc01_2D_x , myThid )
      _EXCH_XY_R8 ( diffKhwc01_2D_y , myThid )

            WRITE(standardMessageUnit,'(A,2I4,/,2f5.2)')
     & 'smooth 2D default parameters: ',
     & wc01_2D_nbt(smoothOpNbCur),wc01_2D_T,
     & wc01_2D_Lx0(smoothOpNbCur),wc01_2D_Ly0(smoothOpNbCur)

cgf write the diffusion operator in files
cgf ... that will eventually be indexed
c      write(fnamegeneric(1:80),'(1a)')
c     &    'wc01_2Doperator'
      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_2Doperator',smoothOpNbCur
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & 1,diffKhwc01_2D_x,1,1,mythid)
      call mdswritefield(fnamegeneric,64,.false.,'RL',
     & 1,diffKhwc01_2D_y,2,1,mythid)

      end