# Copyright (C) 2011, 2012 David Maxwell
# 
# This file is part of PISM.
# 
# PISM is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
# 
# PISM is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
# 
# You should have received a copy of the GNU General Public License
# along with PISM; if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

import PISM

def computeSIASurfaceVelocities(modeldata,siasolver=PISM.SIAFD):
  """Generates surface horizontal velocities corresponding to solving the SIA with zero
  basal sliding.
  
  :param `modeldata`: :class:`PISM.model.ModelData` containing variables and model physics
  :param `siasolver`: specific class used for solving the SIA
  """
  md = modeldata
  grid = md.grid
  sia = siasolver(md.grid, md.enthalpyconverter, md.config)
  sia.init(md.vecs.asPISMVars())

  zero_sliding = PISM.IceModelVec2V()
  zero_sliding.create(grid, 'basal_velocity', False )
  zero_sliding.set(0.)

  sia.update(zero_sliding,False)
  (u,v) = sia.get_horizontal_3d_velocity()

  vel_sia = PISM.model.create2dVelocityVec(grid,name="_sia",stencil_width=1)
  tmp = PISM.IceModelVec2S()
  tmp.create(grid,'tmp',False)
  u.getSurfaceValues(tmp,md.vecs.thickness)
  vel_sia.set_component(0,tmp)
  v.getSurfaceValues(tmp,md.vecs.thickness)
  vel_sia.set_component(1,tmp)
  
  return vel_sia