#include "CPP_OPTIONS.h"

      subroutine smooth_ensemble3D (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)
      _RL fld_mean(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
      _RL fld_meansq(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
      _RL fld_norm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
      character*( 80) fnamegeneric
      integer i,j,k,bi,bj
      integer itlo,ithi
      integer jtlo,jthi
      integer myThid
      _RL port_rand, port_rand_norm
      integer iii,ii,nbRand

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

      

C initialize the random number generator:
      fld_in(1,1,1,1,1)=port_rand(1)
      nbRand=1000

C main loop:
      DO ii=1,nbRand
            WRITE(standardMessageUnit,'(A,I4)')
     & 'smooth_ensemble3D ii: ',ii

C define a perturbation:
      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           fld_in(i,j,k,bi,bj)=0.
           if (maskC(i,j,k,bi,bj).NE.0) then
           fld_in(i,j,k,bi,bj)=port_rand_norm()
           endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_RL ( fld_in, myThid )

C the smoothing itself:
      call smooth_correl3D(fld_in,1,mythid)

      write(fnamegeneric(1:80),'(1a,i3.3)')
     &    'wc01_3Densemble.',ii
      call mdswritefield(fnamegeneric,32,.false.,'RL',
     & nR,fld_in,1,1,mythid)

      ENDDO


      end