module biosun
  implicit none
  private
  
  public :: shortwaverad
contains

subroutine shortwaverad (i,j,xkoor,ykoor,cloudi,hum,casus,swrad)

  !-----------------------------------------------------------------------------
  ! Short wave radiation subroutine for biomod. Based on subroutine q_i
  ! 
  ! Outputs incoming radiation at the surface
  !-----------------------------------------------------------------------------

  !- modules  ------------------------------------------------------------------
  use constants,  only : zero, milli, half, one, twentyfour, zeitstu, pair,    &
                         zwopi

  implicit none

  !- arguments -----------------------------------------------------------------
  logical, intent(in)                 :: casus
  integer(4), intent(in)              :: i,j
  real(8), intent(in)                 :: xkoor(:), ykoor(:)             !(n) (m)
  real(8), intent(in)                 :: hum, cloudi 
  real(8), intent(out)                :: swrad
  
  !
  !  Formel von Zillman, siehe Maykut,
  !       Introduction to Ice in the Polar Oceans, 1985, S.80,81
  !
  real(8), parameter :: eighty = 80.0_8
  real(8), parameter :: a1 = 0.622_8, a2 = 0.378_8, a3 = 1.2_8
  real(8), parameter :: a4 = 0.045_8, a5 = 0.6_8
  real(8), parameter :: q0        = 1350.0_8
  real(8), parameter :: sindek    = 0.39870_8    ! sindek = sin (23.5_8 grad)
  real(8), parameter :: albedo0   = 0.070_8
  real(8), parameter :: albedoi   = 0.600_8
  real(8), parameter :: ddayspryr = 1.0_8/365.2422_8

  real(8)    :: phi, rlam, albedo, hd24, dekl, ea, sa, sw, sd, cd, ph
  real(8)    :: humid, cloud
 
  hd24 = zeitstu/twentyfour
  dekl = asin( sindek*sin((hd24-eighty)*ddayspryr*zwopi) )

  sd = sin(dekl)
  cd = cos(dekl)
 
  ph = (hd24-half)*zwopi
  
  !- albedo:
  if (casus) then
    albedo = albedo0
  else
    albedo = albedoi
  endif

  !- time-center humidity and cloud cover:
  humid = hum*milli
  cloud = cloudi

  !- lat/lon coords:
  phi  = ykoor(i)
  rlam = xkoor(j)

  !- sun angle:
  sa = sin(phi)*sd + cos(phi)*cd*cos(ph + rlam)
  sa = max(sa,zero)
    
  !- short wave radiation:
  ea = pair*humid/(a1 + a2*humid)
  sw = q0*sa*sa/(a3*sa + (one+sa)*ea*milli + a4) 
  swrad = sw*(one - a5*cloud)*(one - albedo)

end subroutine shortwaverad

end module biosun
