!   read   day hour Qs(MJ/hour)
!   output  day hour Qs  Qs_direct  Qs_diffuse Solar_Angle Solar_Adimuth_Angle
!       Qses in W/m2
!
!

program Q50_Sr_Diff_Dir
implicit none

 integer :: year, month, day
 real    :: hour
 integer :: yday
 real    :: elevation, azimuth
 real    :: qs, qsdir, qsdif
 integer :: ios
 integer :: day_of_year

 real, parameter :: factor=(1.0e6/3600.)    !   MJ/hour --> W/m2

   year = 2005
   month = 8

   do while( 1==1  )
     read(*,*, END=900) day, hour,  qs
     qs = qs*factor
      yday = day_of_year(year, month, day)
      call  solang2( yday, hour, elevation, azimuth )
      call  dirdif( qs, elevation, qsdir, qsdif )
     write(*,*) day, hour, qs,  qsdir, qsdif, elevation, azimuth
  end do
  900 stop

end

subroutine dirdif( qs, elevation, qsdir, qsdif )
!  from de Jong 1980
!     via Spitters et. al. 1986


 real :: qs, elevation, qsdir, qsdif
 real, parameter ::  qs00=1370.0    ! Solar constant (W/m2)
 real, parameter :: pi=3.14159265358979323846
 real :: qs0, ratio, k, r



	if( elevation > 0 .AND. elevation < pi .AND. qs > 0) then
	  qs0 = qs00*sin( elevation )
	  ratio = qs/qs0
	    r = 0.847-1.61*sin(elevation)  &
	    		+ 1.04*sin(elevation)*sin(elevation)
	    k = (1.47-r)/1.66

	  if(ratio<=0.22) then
	    qsdif=qs
	    qsdir=0
	  else if(ratio>0.22 .AND. ratio <=0.35) then
	    qsdif=qs*(1.0-6.4*(ratio-0.22)*(ratio-0.22))
	    qsdir=qs-qsdif
	  else if(ratio>0.35 .AND. ratio<=k) then
	    qsdif=qs*(1.47-1.66*ratio)
	    qsdir=qs-qsdif
	  else
	    qsdif=r
	    qsdir=qs-qsdif
	  end if


	else
	  qsdir = 0
	  qsdif = 0
	end if

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
    	

