#include <misc.h>
#include <params.h>
subroutine srfsice(oro     ,snowh   ,ubot    ,vbot    ,tbot    , &
                   qbot    ,thbot   ,zbot    ,pmidm1  ,srfrad  , &
                   tssub   ,qflx    ,taux    ,tauy    ,ts      , &
                   shflx   ,lhflx   ,lwup    ,tref    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Compute sea ice to atmosphere surface fluxes; then compute
! sea ice temperature change.
!
! Method: 
! Temperatures over sea-ice surfaces are specified in 'plevmx' layers of
! fixed thickness and thermal properties.  The forecast temperatures are
! determined from a backward/implicit diffusion calculation using
! linearized sensible/latent heat fluxes. The bottom ocean temperature
! is fixed at -2C, allowing heat flux exchange with underlying ocean.
! Temperature over sea ice is not allowed to exceed melting temperature.
! 
! Author: CCM1
! 
!-----------------------------------------------------------------------
!
! $Id: srfsice.F90,v 1.1.6.1 2000/08/31 23:22:07 rosinski Exp $
! $Author: rosinski $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
  implicit none
#include <compbl.h>
#include <comtsc.h>

!------------------------------Arguments--------------------------------
  real(r8), intent(in) :: oro(plond)           ! Land/ocean/seaice flag
  real(r8), intent(in) :: snowh(plond)         ! Snow depth (liquid water equivalent)
  real(r8), intent(in) :: ubot(plond)          ! Bottom level u wind
  real(r8), intent(in) :: vbot(plond)          ! Bottom level v wind
  real(r8), intent(in) :: tbot(plond)          ! Bottom level temperature
  real(r8), intent(in) :: qbot(plond)          ! Bottom level specific humidity
  real(r8), intent(in) :: thbot(plond)         ! Bottom level potential temperature
  real(r8), intent(in) :: zbot(plond)          ! Bottom level height above surface
  real(r8), intent(in) :: pmidm1(plond)        ! Bottom level pressure
  real(r8), intent(in) :: srfrad(plond)        ! Srf solar abs flux plus down longwave
  real(r8), intent(inout):: tssub(plond,plevmx)  ! Surface/sub-surface temperatures
  real(r8), intent(out):: qflx(plond,pcnst+pnats)    ! Constituent flux (kg/m2/s)
  real(r8), intent(out):: taux(plond)          ! X surface stress (N/m2)
  real(r8), intent(out):: tauy(plond)          ! Y surface stress (N/m2)
  real(r8), intent(out):: ts(plond)            ! surface temperature (K)
  real(r8), intent(out):: shflx(plond)         ! Surface sensible heat flux (J/m2/s)
  real(r8), intent(out):: lhflx(plond)         ! Surface latent   heat flux (J/m2/s)
  real(r8), intent(out):: lwup(plond)          ! surface longwave up flux (W/m2)
  real(r8), intent(out):: tref(plond)          ! 2m reference temperature
!-----------------------------------------------------------------------

!---------------------------Local variables-----------------------------
  integer indx(plond)       ! Longitude index array (land)
  integer m                 ! constituent index
  integer npts              ! Number of land points
  integer i,ii              ! Longitude indices
  integer k                 ! Sub-surface level index
  integer isrftyp(plond)    ! Integer type for surface

  real(r8) tsbsf(plond,plevmx)  ! Non-adjusted srfc/sub-srfc temperatures
  real(r8) fnt(plond)           ! Net surface flux for input conditions
  real(r8) dfntdt(plond)        ! ts partial derivative of net srf flx
  real(r8) ltheat(plond)        ! Latent heat for given sfc conditions
!-----------------------------------------------------------------------

!--------------------------Statement Functions--------------------------
  logical seaice
  seaice(i) = nint(oro(i))==2
!-----------------------------------------------------------------------
!
! Find sea ice points
!
  npts = 0
  do i=1,plon
     if (seaice(i)) then
        npts = npts + 1
        indx(npts) = i
     end if
  end do
  if (npts==0) return
!
! Set latent heat for evaporation sea ice surface
!
  do ii=1,npts
     i = indx(ii)
     isrftyp(i) = 2
     ltheat(i) = latvap + latice
  end do
!
! Compute surface fluxes, derviatives, and exchange coefficiants
!
  call flxsice(indx    ,npts    ,pmidm1  ,ubot    ,vbot    , &
               tbot    ,qbot    ,thbot   ,zbot    ,srfrad  , &
               tssub(1,1),ltheat,fnt     ,dfntdt  ,shflx   , &
               lhflx   ,taux    ,tauy    ,tref    )
!
! Initialize surface/subsurface temperatures for srftsb
!
  do k=1,plevmx
     do ii=1,npts
        i = indx(ii)
        tsbsf(i,k) = tssub(i,k)
     end do
  end do
!
! Diffusion calculation for temperature
!
  call srftsb(isrftyp ,indx    ,npts    ,fnt     ,dfntdt  , &
              snowh   ,tsbsf   )
! 
! Modification to updated surface temperatures
! Reset temperature to melting point.
!
  do ii=1,npts
     i = indx(ii)
     do k=1,plevmx
        tsbsf(i,k) = min(tsbsf(i,k),tmelt)
     end do
  end do
!
! Update surface and sub-surface temperatures
!
  do k=1,plevmx
     do ii=1,npts
        i = indx(ii)
        tssub(i,k) = tsbsf(i,k)
     end do
  end do
  do ii=1,npts
     i = indx(ii)
     ts(i) = tssub(i,1)
     lwup(i) = stebol * ts(i)**4
  end do
!
! Evaluate constituent fluxes
!
  do ii=1,npts
     i = indx(ii)
     qflx(i,1) = lhflx(i)/ltheat(i)
  end do
!
! Set non-water constituent fluxes to zero
!
  do m=2,pcnst
     do ii=1,npts
        i = indx(ii)
        qflx(i,m) = 0.
     end do
  end do
!    
  return
end subroutine srfsice

