subroutine initialize_user
#define DEBUG_LEVEL 3

use mpi_parameters
use user_parameters
use dependent_variables
use dimensional_scales
use immersed_boundary
use grid_info
use counters_flags_etc

use dimensional_scales
use pde_parameters
implicit none

include '../input/problem_size.h'
include 'mpif.h'

real       :: xval,yval,zval,rho_total
real       :: near_surface_window
real       :: near_x1_window
real       :: near_x0_window,window_scale

if( user_subs_flag == 'no' ) return

#if DEBUG_LEVEL >= 1
 if( myid==0) then
  write(0,*) ' '
  write(0,*) ' '
  write(0,*) 'hello world from initialize_user'
  write(0,*) ' setting parameters for surface BCs'
 endif
#endif

 x_0 = 600.
 x_1 = 1200.            ! [m]   left ice/water interface
!!! x_2 = x_1+200.        ! [m]   right ice/water interface
 x_2 = x_1
 V_0 = 0.40           ! [m/s] max ice speed at x_1
 window_scale=5.0

allocate( surface_V_of_x(nx) )

! define surface velocity for sheared ice motion
  do i=1,nx
   xval = Grid(0)%x(i,1)*L_scale
    if(xval >= x_0) then
     near_x0_window = 1.0 - exp( -( (xval-x_0)/(window_scale) )**2  )  ! decay like gaussian w/ given scale
    else
     near_x0_window = 0.0
    endif
    if( xval >= x_1 ) then
     near_x1_window = exp( -( (xval-x_1)/(window_scale) )**2  )  ! decay like gaussian w/ given scale
    else
     near_x1_window = 1.0
    endif
     surface_V_of_x(i) = -V_0*((xval-x_0)/x_1)*near_x1_window*near_x0_window/U_scale
  enddo

!========================================================
!ENABLE VERTICALLY VARYING EDDY VISCOSITY AND DIFFUSIVITY
!
allocate( Km_of_z(locnz), Kt_of_z(locnz) )

  call set_diffusion

!========================================================

if( restart_flag == 'yes' ) return
!
  print*,'nx=',nx,'  ny=',ny,'  locnz=',locnz
! we can also specify any non-zero initial conditions
! here if desired, I've just re-set zeros here to illustrate
  do i=1,nx
!!! xlv Quelle valeur pour k ?  
!!!   xval = Grid(0)%x(i,k)*L_scale   ! dimensional x position
!!!
   do j=1,ny
    yval = Grid(0)%y(j)*L_scale    ! dimensional y position
     do k=1,locnz
!!! xlv
      xval = Grid(0)%x(i,k)*L_scale   ! dimensional x positio
!!!      
      zval = Grid(0)%z(i,k)*L_scale   ! dimensional z position

      near_surface_window = exp( -( (zval-L_scale)/(5.0) )**2  )  ! decay like gaussian w/ 5m scale
      v(i,j,k)= 0.0*surface_V_of_x(i)*near_surface_window
      u(i,j,k)=0.0/U_scale
      w(i,j,k)=0.0/U_scale

!!!      s1(i,j,k)=0.0/s1_scale
!!!      s2(i,j,k)=0.0/s2_scale
      s1(i,j,k)=0.0
      s2(i,j,k)=0.0

!     Contravariant rather than cartesian components are actually needed
!     I'm assuming a stretched but cartesian grid for this problem
!     no transformations needed for v or pd
      u(i,j,k) = u(i,j,k)/Grid(0)%x_xi(i,k)       ! u --> U
      w(i,j,k) = w(i,j,k)/Grid(0)%z_zeta(i,k)     ! w --> W

     enddo
    enddo
   enddo

end subroutine initialize_user