  PROGRAM alb_mk
! 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)
     call albedo( coszrs,  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
       endif
      end subroutine albedo


