program test

 integer :: year, month, day    ! month (1-12)
 integer :: md(12)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
 integer, external :: day_of_year
 integer :: yday
 real    :: hour, elevation, azimuth, fact, lat, lon
 real    :: elevation3, azimuth3
 integer :: i, j

 fact = 3.1415/180.0

 year = 2008; month = 7; day = 1
 lat = (35+10.0/60.0)*fact; lon = (136+57.9/60.0)*fact    ! Nagoya (JMA)

! do month=1, 12

  yday = day_of_year(year, month, day)
  write(*,*) 'yday=', yday

  do i=10, 12
  hour = i
  call  solang3( yday, hour, lat, lon, elevation, azimuth )

   write(*,*) year, month, day, hour,  elevation/fact, azimuth/fact
  end do
!  end do

end





integer function day_of_year( year, month, day )
implicit none
! calculate day of year ( Jan.1 = 1 ... )
 integer :: year, month, day    ! month (1-12)
 integer :: md(12)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
 integer :: i, result

 if( mod(year,4)==0 ) then
   if( mod(year, 100)==0 ) then
     if( mod(year, 400)==0 ) then
        md(2)=29
     else
        md(2)=28
     end if
   else
       md(2)=29
   end if
 else
       md(2)=28
 end if
  
 result = 0
 do i=1, month-1
  result = result+md(i)
 end do

 result = result + day

 day_of_year = result

end function


 



subroutine solang3( yday, hour, lat, lon, elevation, azimuth )
implicit none

!  calculate solar elevation angle and azimuthal angle
!   using day_of_year, and time (hour)
! latitude and longitude of the place is set below
!  also latitude of the place for time   (JST)
!    elevation, azimuth  in radian
!  lat lon in radian


!  from Campbell, Norman 2003

 double precision,parameter :: pi = 3.14159265358979323846
 integer :: yday
 real    :: hour, elevation, azimuth

 real    :: lat, lon
 !real    :: lat = 35+10.0/60.0, lon = 136+57.9/60.0    ! Nagoya (JMA)
 real    :: lat1, lon1
 real    :: timelon = 135.0

 real(kind=8)    :: declination ! Solar declination  (SEKII)
 real(kind=8)    :: sin_declination, cos_declination ! sin and cos of Solar declination  (SEKII)
 real(kind=8)    :: f, f_yday, ff
 real(kind=8)    :: A, B, et, et1

 real    :: localtime, t
 real(kind=8)    :: coszenith, sinkai, coskai, kai, kai1

  f = pi/180.0_8

  lon1 = lon/f 

   A = 2 * pi * yday / 365

!    declination = 0.006918 - 0.399912*cos(A)+0.070257*sin(A) &
!       -0.006758*cos(2*A) + 0.000907*sin(2*A) & 
!       - 0.002697*cos(3*A) + 0.001480 * sin(3*A) 

   sin_declination = 0.39785_8 * sin(f*(278.97_8 + 0.9856_8 * yday &
      + 1.9165_8*sin(f*356.6_8+0.9856_8*yday)))
   
   if( 1.0_8 - sin_declination*sin_declination > dble(0) ) then
     cos_declination = sqrt( 1.0_8 - sin_declination*sin_declination)
   else
     cos_declination = 0.0_8
   end if

!   write(*,*) 'DEC', declination/f, asin(sin_declination)/f
!   declination = asin(sin_declination)

   localtime = hour + (lon1-timelon)/15

    B = 2*pi*yday/365
    et1 = (0.000075 + 0.001868 * cos(B)  -0.032077*sin(B) &
           -0.014615 * cos(2*B) -0.040849 * sin(2*B) ) * 12 / pi

    ff = 279.575_8+0.9856_8*yday

     et = (-104.7_8*sin(f*ff)+596.2_8*sin(f*2*ff)+4.3_8*sin(f*3*ff)-12.7_8*sin(f*4*ff) &
             -429.3_8*cos(f*ff)-2.0_8*cos(f*2*ff)+19.3_8*cos(f*3*ff))/dble(3600)

! write(*,*) 'ET=', et1, et
    t = 15 * ( localtime + et - 12 )

    coszenith = sin_declination*sin(lat) &
                            +cos_declination*cos(lat)*cos(f*t)
    elevation = 0.5*pi - acos(coszenith)

    sinkai = cos_declination*sin(f*t)/sin(0.5_8*pi-elevation)
    coskai = (-cos(lat)*sin_declination+sin(lat) &
            *cos_declination*cos(f*t))/sin(0.5_8*pi-elevation)

if(abs(sinkai)<= 1.0_8) then
    kai = real(asin( sinkai )) 
else if(sinkai> 1.0_8) then
    kai = real( pi/2.0_8 )
else if(sinkai< -1.0_8) then
    kai = real (2.0_8*pi - pi/2.0_8)
else
    write(0,*) 'Never here'
end if
!    write(*,*) 'kai', sinkai, coskai, kai

    if( coskai < 0 ) then
    	kai1 = pi - kai
    else if( coskai > 0 .AND. sinkai < 0 ) then
    	kai1 = 2*pi+kai
    else
    	kai1 = kai
    end if

    kai1 = kai1+pi
    if(kai1>2*pi) then
    	kai1 = kai1 - 2*pi
    end if

    azimuth = kai1

!    write(*,*) 'delta', (180/pi)*declination, &
!    			'time angle', (180/pi)*t, 'et', et

end subroutine
    	

