! Test of albedo subroutine: (remove these lines later)
angle = 0.
do while (angle.lt.90.)
coszrs = cos(acos(-1.)/180.*angle)
call albedo(.true., coszrs, 290., asdir, aldir, asdif, aldif)
print*, 'angle=',angle, '  asdir=',asdir, '  aldir=',aldir, '  
asdif=', asdif, '  aldif=',aldif
angle = angle + 1.
end do
end

subroutine albedo(coszrs, asdir, aldir, asdif, aldif)
!-----------------------------------------------------------------------
! Computes surface albedos over ocean, ice
! and land  (the latter is added by Marat Khairoutdinov from CCM3)
! Modified version taken from SAM6.8 (Oct 2010) - MK

! Two spectral surface albedos for direct (dir) and diffuse (dif)
! incident radiation are calculated. The spectral intervals are:
!   s (shortwave)  = 0.2-0.7 micro-meters
!   l (longwave)   = 0.7-5.0 micro-meters
!
! Uses knowledge of surface type to specify albedo, 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.
!
! 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).

real, intent( in) :: coszrs   ! Cosine of solar zenith angle
real, intent(out) :: asdir, & ! Srf alb for direct rad   0.2-0.7 micro-ms
aldir, & ! Srf alb for direct rad   0.7-5.0 micro-ms
asdif, & ! Srf alb for diffuse rad  0.2-0.7 micro-ms
aldif    ! Srf alb for diffuse rad  0.7-5.0 micro-ms
real, parameter :: adif = 0.06
!-----------------------------------------------------------------------

!
!
! Ice-free ocean albedos function of solar zenith angle only, and
! independent of spectral interval:
!
if(coszrs <= 0) then
aldir = 0; asdir = 0; aldif = 0; asdif = 0
else
aldir = ( .026 / (coszrs**1.7 + .065)) + &
(.15*(coszrs - 0.10) * (coszrs - 0.50) * (coszrs - 1.00) )
asdir = aldir
aldif = adif
asdif = adif
end subroutine albedo