太陽高度角と方位角、散乱、直達への分配

太陽高度角と方位角の計算

研究室の皆さん

玉川です。太陽高度角と方位角を計算するプログラムを書いてみました。

計算方法は、  太陽天頂角・太陽方位角の算出方法  (九州大学の村上拓彦氏のWWW頁)によります。



以下につけたようなもので

最初の program 〜 end までは、テスト用のメインプログラムです。

その下以降を自分のソフトに組み込んで使って下さい。

使い方は、見れば分かりますよね。サブルーチン自体の
計算結果はラディアンですので、注意して下さい。

2008.1.18



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
    	




 散乱、直達への分配 

 日射量を読んで、散乱、直達、太陽高度、太陽方位角 を
出力するプログラム例 
 (無テストなので、注意して使用の事)
放射の分配は、
Spitters, C.J.T., Toussaint, H.A.J.M. and Goudrian, J.
Separating the diffuse and direct component of global radiation and its implicat
ions for modeling canopy photosynthesis, part I. components of incoming radiatio
n
Agricultural and forest meteorology, 38,(1986) 217-229

による。

2008.1.21 

関数部分だけ C 版も書いてみました。 ここ  、
上のC版と合わせると使えると思います。(2009.4.1)