#include <misc.h>
#include <params.h>

subroutine albocean(oro     ,snowh   ,coszrs  ,asdir  ,aldir   , &
                    asdif   ,aldif   )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Compute surface albedos
!
! Method: 
! Computes surface albedos for direct/diffuse incident radiation for
! two spectral intervals:
!   s = 0.2-0.7 micro-meters
!   l = 0.7-5.0 micro-meters
!
! Albedos specified as follows:
! Ocean           Uses solar zenith angle to compute albedo for direct
!                 radiation; diffuse radiation values constant; albedo
!                 independent of spectral interval and other physical
!                 factors such as ocean surface wind speed.
! Ocean with      Surface albs specified; combined with overlying snow
!   sea ice       
!
! For more details , see Briegleb, Bruce P., 1992: Delta-Eddington
! Approximation for Solar Radiation in the NCAR Community Climate Model,
! Journal of Geophysical Research, Vol 97, D7, pp7603-7612).
! 
! Author: CCM1
! 
!-----------------------------------------------------------------------
!
! $Id: albocean.F90,v 1.1.6.1 2000/08/31 23:22:06 rosinski Exp $
! $Author: rosinski $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
  implicit none
#include <albedo.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) :: coszrs(plond)    ! Cosine solar zenith angle
  real(r8), intent(out):: asdir(plond)     ! Srf alb for direct rad   0.2-0.7 micro-ms
  real(r8), intent(out):: aldir(plond)     ! Srf alb for direct rad   0.7-5.0 micro-ms
  real(r8), intent(out):: asdif(plond)     ! Srf alb for diffuse rad  0.2-0.7 micro-ms
  real(r8), intent(out):: aldif(plond)     ! Srf alb for diffuse rad  0.7-5.0 micro-ms
!-----------------------------------------------------------------------

!---------------------------Local variables-----------------------------
  integer i                 ! Longitude index
  real(r8) frsnow           ! Horizontal fraction of snow cover
  real(r8) snwhgt           ! Physical snow height
  real(r8) rghsnw           ! Roughness for horizontal snow cover fractn
  real(r8) sasdir(plond)    ! Snow alb for direct rad  0.2-0.7 micro-ms
  real(r8) saldir(plond)    ! Snow alb for direct rad  0.7-5.0 micro-ms
  real(r8) sasdif(plond)    ! Snow alb for diffuse rad  0.2-0.7 micro-ms
  real(r8) saldif(plond)    ! Snow alb for diffuse rad  0.7-5.0 micro-ms
!-----------------------------------------------------------------------

!--------------------------Statement functions--------------------------
  logical ocean,seaice
  ocean(i) = nint(oro(i))==0
  seaice(i) = nint(oro(i))==2
!-----------------------------------------------------------------------
!
! Initialize all ocean/sea ice surface albedos to zero
!
  do i=1,plon
     if (ocean(i) .or. seaice(i)) then
        asdir(i) = 0.
        aldir(i) = 0.
        asdif(i) = 0.
        aldif(i) = 0.
     end if
  end do

  do i=1,plon
     if (seaice(i) .and. coszrs(i)>0.0) then
        asdir(i)  = sices
        aldir(i)  = sicel
        asdif(i) = asdir(i)
        aldif(i) = aldir(i)
        sasdif(i) = snws
        saldif(i) = snwl
     end if
  end do

  do i=1,plon
     if (seaice(i)) then
        if (snowh(i)>0. .and. coszrs(i)>0.) then
           if (coszrs(i)<0.5) then
!
! Zenith angle regime 1 ( coszrs < 0.5 ).
! Set direct snow albedos (limit to 0.98 max)
!
              sasdir(i) = min(0.98_r8,sasdif(i) + (1. - sasdif(i))*0.5* &
                   ((3./(1. + 4.*coszrs(i))) - 1.))
              saldir(i) = min(0.98_r8,saldif(i) + (1. - saldif(i))*0.5* &
                   ((3./(1. + 4.*coszrs(i))) - 1.))
           else
!
! Zenith angle regime 2 ( coszrs >= 0.5 )
!
              sasdir(i) = snws
              saldir(i) = snwl
           end if
!
! Compute both diffuse and direct total albedos
!
           snwhgt = 20.*snowh(i)
           rghsnw = 0.25
           frsnow = snwhgt/(rghsnw + snwhgt)
           asdir(i)  = asdir(i) *(1. - frsnow) + sasdir(i) *frsnow
           aldir(i)  = aldir(i) *(1. - frsnow) + saldir(i) *frsnow
           asdif(i) = asdif(i)*(1. - frsnow) + sasdif(i)*frsnow
           aldif(i) = aldif(i)*(1. - frsnow) + saldif(i)*frsnow
        end if
     end if
  end do
!
! Ice-free ocean albedos function of solar zenith angle only, and
! independent of spectral interval:
!
  do i=1,plon
     if (ocean(i).and.coszrs(i)>0.) then
        aldir(i)  = (.026/(coszrs(i)**1.7 + .065)) + &
             (.15*(coszrs(i) - 0.10)*(coszrs(i) - 0.50)*(coszrs(i) - 1.00))
        asdir(i)  = aldir(i)
        aldif(i) = 0.06
        asdif(i) = 0.06
     end if
  end do
!
  return
end subroutine albocean

