subroutine s1_bar_func(xprime,yprime,zprime,s1_bar,
     *                        s1_bar_x,s1_bar_y,s1_bar_z,g2s1b,
     *                        Lz,s1_scale)
c
c     *****input arguments xprime,yprime,zprime are DIMENSIONLESS as are all returned values
c     (xprime,yprime,zprime):     d'less spatial position
c     Lz:                         vertical size of computational domain [m]
c     s1_scale:                   characteristic range of s1 values, used to scale s1
c                                 units are dimensional, e.g. deg C for s1-->Temp.
      
c     *****output variables, these need to be specified in this routine
c             s1_bar, s1_bar_x, s1_bar_y, s1_bar_z, g2s1b
c             scalar concentration, and spatial gradients, g2s1b=grad2(s1_bar)   [dimensional units]

      implicit none
      real xprime,yprime,zprime
      real Lz,s1_scale
      real s1_bar,s1_bar_x,s1_bar_y,s1_bar_z,g2s1b

c     local variables: you may need some intermediate variables,
c                      these need to be explicitly declared here
      real x,y,z,ans,ansp,ansm,delta_z,delta_x

c     it may be convenient to know the dimensional spatial location (x,y,z) in [m], so we compute it
      x = xprime*Lz       ! use dimensional value to evaluate formulae
      y = yprime*Lz       ! use dimensional value to evaluate formulae
      z = zprime*Lz       ! use dimensional value to evaluate formulae


c****    USER CUSTOMIZED CODE SEGMENT STARTS  HERE ************************
c            scalars='TS' --> s1 = T' and s2 = S'

      delta_z = Lz/10000.    ! [m] (offset for estimating derivs)
      delta_x = 100.   ! ""
      
      ! stay underwater
      if(z >= Lz-delta_z) z=Lz-delta_z
      
      call interp2d(x,z,'T',ans) 
      call interp2d(x,z+delta_z,'T',ansp)
      call interp2d(x,z-delta_z,'T',ansm)
cc      print*,'after call interp2d on z'
      
      s1_bar = ans ! [deg C]
      s1_bar_z = (ansp-ansm)/(2.*delta_z)        ! [deg C/m]
      g2s1b = (ansp -2.*ans + ansm)/(delta_z**2) ! [deg C/m2]
      
      
      call interp2d(x+delta_x,z,'T',ansp)
      call interp2d(x-delta_x,z,'T',ansm)
cc      print*,'after call interp2d on x'

      
      s1_bar_x = (ansp-ansm)/(2.*delta_x)        ! [deg C/m]
         
         
      s1_bar_y = 0.0                             ! [deg C/m] 
      

c****    USER DEFINED CODE SEGMENT ENDS HERE ********************************




c***  LEAVE THIS SECTION AS IS, SCALING MUST BE CONSISTENT********************
c***  WITH THE REST OF THE ALGORITHM  ****************************************
      s1_bar=s1_bar/(s1_scale)
      s1_bar_x = s1_bar_x/(s1_scale/Lz)
      s1_bar_y = s1_bar_y/(s1_scale/Lz)
      s1_bar_z = s1_bar_z/(s1_scale/Lz)
      g2s1b = g2s1b/(s1_scale/Lz**2)
c*****************************************************************************
c*****************************************************************************

c      print*,'s1_bar interpoled ...'
      return
      end

c It may seem odd to have the inputs in dimensionless form
c when it is usually convenient for the user to define the
c ambient quantities in dimensional form. Also, we have a
c **do not touch** section of code that converts the user
c specified dimensional quantities back to dimensionless form.
c I chose to do it this way because this function gets called
c from so many places in the code, it would be a pain and
c a potential source of typos to do the scalings just before
c and just after each call. This is not the case with 
c user_defined_ics and user_defined_forcing. These routines
c get called in one location only and so it makes sense to
c make the scaling "invisible" to the user.