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
 integer :: i, j

 fact = 3.1415/180

 year = 2008; month = 7; day = 1

  yday = day_of_year(year, month, day)

  do i=0, 24
  hour = i
  call  solang2( yday, hour, elevation, azimuth )

   write(*,*) year, month, day, hour, elevation/fact, azimuth/fact
  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 solang2( yday, hour, 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

!  from http://ffpsc.agr.kyushu-u.ac.jp/forman/muratac/solar/solpos1.html
!      it say 'from 6s'

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

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

 real    :: declination ! Solar declination  (SEKII)
 real    :: A, B, et

 real    :: localtime, t
 real    :: coszenith, sinkai, coskai, kai, kai1


   lat1 = lat * pi/180
   lon1 = lon * pi/180

   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) 

   localtime = hour + (lon-timelon)/15

    B = 2*pi*yday/365
    et = (0.000075 + 0.001868 * cos(B)  -0.032077*sin(B) &
           -0.014615 * cos(2*B) -0.040849 * sin(2*B) ) * 12 / pi
    t = 15 * (pi/180) * ( localtime + et - 12 )

    coszenith = sin(declination)*sin(lat1) &
                            +cos(declination)*cos(lat1)*cos(t)
    elevation = 0.5*pi - acos(coszenith)

    sinkai = cos(declination)*sin(t)/sin(0.5*pi-elevation)
    coskai = (-cos(lat1)*sin(declination)+sin(lat1) &
            *cos(declination)*cos(t))/sin(0.5*pi-elevation)
    kai = asin( sinkai ) ;

    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
    	

