*"%表題   日射の計算(春分点を基準に日射量だけ計算する)
*%
*"%履歴   91/02/22 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*"%       91/03/14 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*"%       91/03/21 佐藤正樹       C:\FORT\NISSHA\NISSHA2.FOR
*"%         F は倍精度だと領域オーバーになってしまうので
*"%         単精度に変更した
*-----------------------------------------------------------------------
      PROGRAM  NISSHA2

      PARAMETER ( NLAT = 90 )   !"緯度配列の次元
      PARAMETER ( NTIME = 60 )  !"時間配列の次元

      DOUBLE PRECISION  PI
      PARAMETER( PI = 3.141592653589793 D 0 ) !"円周率

      DOUBLE PRECISION  ISOLAR  !"太陽の全日射量
      DOUBLE PRECISION  RLONG   !"軌道長半径
      DOUBLE PRECISION  ECCENT  !"離心率
      DOUBLE PRECISION  PHI0    !"近日点黄経(rad)
      DOUBLE PRECISION  THETAP  !"赤道傾斜角(rad)
      DOUBLE PRECISION  YTIME   !"公転周期
      DOUBLE PRECISION  DTIME   !"自転周期

      DOUBLE PRECISION  PHI0A   !"近日点黄経(度)
      DOUBLE PRECISION  THETAPA !"赤道傾斜角(度)

      DOUBLE PRECISION  SOLARC  !"遠日点(R=RLONG)での太陽定数
      DOUBLE PRECISION  XI0     !"XI の春分点での値 * (-1)
      DOUBLE PRECISION  TIME0   !"春分点での時間 * (-1)

      DOUBLE PRECISION  LAT(-NLAT:NLAT)  !"緯度
      DOUBLE PRECISION  TIME(0:NTIME) 
     &                          !"時間(春分点から測る, YTIMEで規格化)

      DOUBLE PRECISION  R      !"太陽からの距離(RLONGで規格化)
      DOUBLE PRECISION  SOLAR  !"太陽定数
      DOUBLE PRECISION  DELTAS !"太陽傾斜角
      DOUBLE PRECISION  PHI    !"近日点から測った黄経
      DOUBLE PRECISION  XI     !"楕円をあらわすパラメーター
      DOUBLE PRECISION  H0     !"日照時間

      DOUBLE PRECISION  PHIHARU(0:NTIME)        !"春分点から測った黄経
C"     DOUBLE PRECISION  F(0:NTIME, -NLAT:NLAT)  ! 日平均日射量
      REAL  F(0:NTIME, -NLAT:NLAT)  !"日平均日射量

      DOUBLE PRECISION  FTIME   !"関数副プログラム
      DOUBLE PRECISION  FXI     !"関数副プログラム

      INTEGER           ILAT    !"DO ループで用いる緯度配列の次元
      INTEGER           ITIME   !"DO ループで用いる時間配列の次元
      DOUBLE PRECISION  TAIYOU  !"一時変数, 太陽の出入り



*"    < 入出力パラメーター >

      INTEGER   IPR  / 9 /           !"入力ユニット番号
      CHARACTER HIPFL*13 / 'NISSHA.PAR' / 
     &                               !"入力パラメーターファイル名
      LOGICAL   LIPEX                !"入力パラメーターファイルの存在


      INTEGER   OR   / 20 /          !"出力ユニット番号
      CHARACTER HFLDAT*13              !"出力ファイル名

      LOGICAL   LOEX                  !"出力ファイルの存在



*     < NAMELIST >

      NAMELIST / FLPARA /   HFLDAT
      NAMELIST / IPARA /    ISOLAR,   RLONG,    ECCENT,
     &                      PHI0A,    THETAPA,
     &                      YTIME,    DTIME

*-----------------------------------------------------------------------
*"    < 外部パラメーター >
*"    << 地球のデータ >>

      HFLDAT = 'NISSHA.DAT'

      ISOLAR  = 3.85  D 26                 ! [W]
      RLONG   = 1.496 D 11                 ! [m]
      ECCENT  = 0.0167 D 0
      PHI0A   = 102.908 D 0                !"[度]
      THETAPA = 23.44   D 0                !"[度]
      YTIME   = 365.0 D 0                  !"[日]
      DTIME   = 24.0 D 0                   !"[時間]


*"    < パラメーターの入力 >

      INQUIRE(FILE=HIPFL, EXIST=LIPEX)
      IF( .NOT. LIPEX ) THEN
        WRITE(6,*) '入力パラメーターファイルが存在しません'
        WRITE(6,*) ' ファイル名: ', HIPFL
        STOP
      END IF

      OPEN(IPR,FILE=HIPFL)
        READ(IPR,FLPARA)
        READ(IPR,IPARA)
      CLOSE(IPR)

      WRITE(6,*) '外部パラメーター'
      WRITE(6,*) 'ISOLAR', ISOLAR
      WRITE(6,*) 'RLONG', RLONG
      WRITE(6,*) 'ECCENT', ECCENT
      WRITE(6,*) 'PHI0A', PHI0A
      WRITE(6,*) 'THETAPA', THETAPA
      WRITE(6,*) 'YTIME', YTIME
      WRITE(6,*) 'DTIME', DTIME

      PHI0    = PHI0A / 180.D0 * PI
      THETAP  = THETAPA / 180.D0 * PI

*-----------------------------------------------------------------------
*"    < 定数値の計算 >

      SOLARC  = ISOLAR / ( 4.D0 * PI * RLONG ** 2 )

      XI0     = 2.D0 * MOD(
     &          ATAN(   SQRT( ( 1.D0 - ECCENT ) / ( 1.D0 + ECCENT ) )
     &                * TAN( PHI0 / 2.D0 )
     &              ) + PI, 
     &          PI
     &          )

      TIME0   = FTIME( XI0, ECCENT )

C      WRITE(6,*) 'SOLARC(W/m2)', SOLARC
C      WRITE(6,*) 'XI0(rad)', XI0
C"     WRITE(6,*) 'XI0(度)', XI0 * 180.D0 / PI
C      WRITE(6,*) 'TIME0', TIME0
C"     WRITE(6,*) 'TIME0(日)', TIME0 * YTIME

*-----------------------------------------------------------------------
*"    < 緯度の設定 >

      LAT(0) = 0.D0
      DO 10 ILAT = 1, NLAT

        LAT(ILAT)  = PI/2.D0 * REAL( ILAT ) / REAL( NLAT )
        LAT(-ILAT) = - LAT(ILAT)

 10   CONTINUE



*"    < 時間の設定 : 春分点を基準にする >

      DO 20 ITIME = 0, NTIME

        TIME(ITIME) = - TIME0 + 1.D0 * REAL( ITIME ) / REAL( NTIME )
        TIME(ITIME) = MOD( TIME(ITIME) + 1.D0, 1.D0 )

 20   CONTINUE

*-----------------------------------------------------------------------
*"    < 計算> 

      DO 100 ITIME = 0, NTIME

        WRITE(6,*) 'ITIME', ITIME

        XI = FXI( TIME(ITIME), ECCENT )
        
        PHI= 2.D0 * ATAN( 
     &                SQRT( ( 1.D0 + ECCENT ) / ( 1.D0 - ECCENT ) )
     &              * TAN( XI / 2.D0 )
     &              )
        
        R = 1.D0 - ECCENT * COS( XI )
        
        SOLAR = SOLARC / R ** 2
        
        DELTAS = ASIN( - SIN( THETAP )
     &                        * SIN( PHI0 + PHI )
     &               )
        
        
        
        DO 200 ILAT = -NLAT, NLAT
          
C          WRITE(6,*) 'ILAT', ILAT
          
          TAIYOU = TAN( LAT(ILAT) ) * TAN( DELTAS )
          
          IF( ABS( TAIYOU ) .LT. 1.D0 ) THEN
            H0 = ACOS( - TAIYOU )
          ELSE IF( TAIYOU .LE. -1.D0 ) THEN
            H0 = 0.D0
          ELSE IF( TAIYOU .GE. 1.D0 ) THEN
            H0 = PI
          ENDIF
          
          F(ITIME, ILAT)
     &    = SOLAR / PI
     &      * (   SIN( H0 ) 
     &          * COS( LAT(ILAT) ) * COS( DELTAS )
     &          + H0
     &          * SIN( LAT(ILAT) ) * SIN( DELTAS )
     &        )

C          WRITE(6,'('' '', A, 5F13.4)')
C"    &          'F(0:NTIME, -NLAT:NLAT)  ! 日平均日射量', 
C     &           F(ITIME, ILAT)

        PHIHARU(ITIME) = PHI + PHI0
        PHIHARU(ITIME) = MOD( PHIHARU(ITIME) + 2.D0 * PI, 2.D0 * PI )

 200    CONTINUE
 100  CONTINUE

*-----------------------------------------------------------------------
*"    < 出力ファイルのオープン >

      INQUIRE(FILE=HFLDAT, EXIST=LOEX)
      IF( LOEX ) THEN
        OPEN ( OR, FILE=HFLDAT, STATUS='OLD' )
        CLOSE( OR, STATUS='DELETE' )
      END IF

      OPEN ( OR, FILE=HFLDAT )

*"    < 出力 >

      CALL  ONISSHA2
     I           ( OR,
     I             ISOLAR,   RLONG,    ECCENT,   PHI0,     THETAP,
     I             YTIME,    DTIME,
     I             SOLARC,   XI0,      TIME0,
     I             LAT,      TIME,     PHIHARU,  F, 
     D             NLAT,     NTIME )

*"    < ファイルのクローズ >

      CLOSE(OR)

*-----------------------------------------------------------------------
      STOP
      END
*=======================================================================
*"%表題   軌道関数 時間
*%
*"%履歴   91/02/23 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*%
*%
*-----------------------------------------------------------------------
      FUNCTION FTIME( XI, ECCENT )

      DOUBLE PRECISION  FTIME   !"時間
      DOUBLE PRECISION  XI      !"パラメーター
      DOUBLE PRECISION  ECCENT  !"離心率

      DOUBLE PRECISION  PI
      PARAMETER( PI = 3.141592653589793 D 0 ) !"円周率

*-----------------------------------------------------------------------

      FTIME = 1.D0 / 2.D0 / PI * ( XI  - ECCENT * SIN( XI ) )

C      WRITE(6,*) 'SUB : FTIME'
C      WRITE(6,*) 'XI, XI0, ECCENT, FTIME'
C      WRITE(6,*)  XI, XI0, ECCENT, FTIME

*-----------------------------------------------------------------------
      RETURN
      END
*=======================================================================
*"%表題   XI の計算
*%
*"%履歴   91/02/22 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*%
*%
*-----------------------------------------------------------------------
      FUNCTION FXI( TIME, ECCENT )

      DOUBLE PRECISION  FXI     !"楕円をあらわすパラメーター

      DOUBLE PRECISION  TIME    !"時間
      DOUBLE PRECISION  ECCENT  !"離心率

      DOUBLE PRECISION  XI1     !"XI の値の初期値1
      DOUBLE PRECISION  XI2     !"XI の値の初期値2

      DOUBLE PRECISION  FTIME   !"関数副プログラム : 時間
      EXTERNAL          FTIME
      DOUBLE PRECISION  FGYAKU  !"関数副プログラム : 逆関数

      DOUBLE PRECISION  PI
      PARAMETER( PI = 3.141592653589793 D 0 ) !"円周率

      DOUBLE PRECISION  ATIME   !"標準化した時間

*-----------------------------------------------------------------------
*"    < TIME の標準化 >
      ATIME = MOD( MOD( TIME, 1.D0 ) + 1.D0, 1.D0 )



*"    < XI の初期値の設定 >

      XI1 = 0.D0
      XI2 = 2.D0 * PI



*"    < XI を求める >

      FXI = FGYAKU( FTIME, ATIME, XI1, XI2, ECCENT )

C      WRITE(6,*) 'FXI, ATIME, XI1, XI2, ECCENT'
C      WRITE(6,*)  FXI, ATIME, XI1, XI2, ECCENT

*-----------------------------------------------------------------------
      RETURN
      END
*=======================================================================
*"%表題   逆関数
*%
*"%履歴   91/02/23 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*%
*%
*"    < X1, X2 との間の2分法によって X を求める >
*-----------------------------------------------------------------------
      FUNCTION FGYAKU( FUNC, Y, X1, X2, PARA )

      DOUBLE PRECISION  FGYAKU  !"逆関数

      DOUBLE PRECISION  FUNC    !"関数
      EXTERNAL          FUNC

      DOUBLE PRECISION  Y       !"関数値
      DOUBLE PRECISION  X1      !"逆関数の値の初期値1
      DOUBLE PRECISION  X2      !"逆関数の値の初期値2

      DOUBLE PRECISION  PARA   !"パラメーター

*"    < SUBROUTINE 内で定義される変数 >

      DOUBLE PRECISION  Y0      !"Y の計算値
      DOUBLE PRECISION  XX1     ! X1
      DOUBLE PRECISION  XX2     ! X2
      DOUBLE PRECISION  XX3     ! X3
      
      DOUBLE PRECISION  SGN     !"X1, X2 での関数値の差の符号

      INTEGER   I            !"DO ループ
      INTEGER   IMAX / 30 /  !"I の最大値

*-----------------------------------------------------------------------
*"    < 範囲 X1, X2 の適合性のチェック >
*"    Y - FUNC(X1)  と Y - FUNC(X2) の符号が逆である
*"    ことが仮定されている

      IF( ( Y - FUNC( X1, PARA ) )
     &  * ( Y - FUNC( X2, PARA ) )
     &    .GT. 0.0D0 ) THEN
        WRITE(6,*) 'SUBROUTINE  FGYAKU エラー'
        WRITE(6,*) '範囲 X1, X2が適切でありません. '
        STOP
      END IF


*"    < X1, X2 での値の大小関係を調べる >

      IF( FUNC( X1, PARA ) - FUNC( X2, PARA )
     &    .LE. 0.0D0 ) THEN
        SGN = 1.D0
      ELSE
        SGN = - 1.D0
      ENDIF


*"    < くりかえし >

      XX1 = X1
      XX2 = X2

      DO 10 I = 1, IMAX

        XX3 = ( XX1 + XX2 ) / 2.0

        Y0 = FUNC( XX3,  PARA )

C        WRITE(6,*) 'I=', I, 'XX1=', XX1, 'XX2=', XX2
C        WRITE(6,*) 'Y0=', Y0, 'Y-Y0=', Y - Y0

        IF( ( Y - Y0 ) * SGN .LE. 0.0D0 ) THEN
          XX2 = XX3
        ELSE
          XX1 = XX3
        ENDIF

 10   CONTINUE

      FGYAKU = XX3

*-----------------------------------------------------------------------
      RETURN
      END
*=======================================================================
*"%表題   日射の計算値の出力(簡易)
*%
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHA.FOR
*"%履歴   91/03/21 佐藤正樹       C:\FORT\NISSHA\NISSHA2.FOR
*%
*%
*-----------------------------------------------------------------------

      SUBROUTINE  ONISSHA2
     I           ( OR,
     I             ISOLAR,   RLONG,    ECCENT,   PHI0,     THETAP,
     I             YTIME,    DTIME,
     I             SOLARC,   XI0,      TIME0,
     I             LAT,      TIME,     PHIHARU,  F, 
     D             NLAT,     NTIME )

      INTEGER  NLAT            !"緯度配列の次元
      INTEGER  NTIME           !"時間配列の次元
      INTEGER  OR              !"出力ユニット番号

      DOUBLE PRECISION  ISOLAR  !"太陽の全日射量
      DOUBLE PRECISION  RLONG   !"軌道長半径
      DOUBLE PRECISION  ECCENT  !"離心率
      DOUBLE PRECISION  PHI0    !"近日点黄経
      DOUBLE PRECISION  THETAP  !"赤道傾斜角
      DOUBLE PRECISION  YTIME   !"公転周期
      DOUBLE PRECISION  DTIME   !"自転周期

      DOUBLE PRECISION  SOLARC  !"遠日点(R=RLONG)での太陽定数
      DOUBLE PRECISION  XI0     !"XI の春分点での値 * (-1)
      DOUBLE PRECISION  TIME0   !"春分点での時間 * (-1)

      DOUBLE PRECISION  LAT(-NLAT:NLAT)  !"緯度
      DOUBLE PRECISION  TIME(0:NTIME) 
     &                          !"時間(春分点から測る, YTIMEで規格化)

      DOUBLE PRECISION  PHIHARU(0:NTIME)        !"春分点から測った黄経
C"     DOUBLE PRECISION  F(0:NTIME, -NLAT:NLAT)  ! 日平均日射量
      REAL  F(0:NTIME, -NLAT:NLAT)  !"日平均日射量

*-----------------------------------------------------------------------

 600  FORMAT(2X, 5F13.4)
 610  FORMAT(2X,  F13.4)
 620  FORMAT(2X, 5E13.4)
 630  FORMAT(2X,  E13.4)
 640  FORMAT(2X,  I13  )
 650  FORMAT(2X,  A13  )

*-----------------------------------------------------------------------

      WRITE(OR,650)   'ISOLAR       '
      WRITE(OR,630)    ISOLAR

      WRITE(OR,650)   'RLONG        '
      WRITE(OR,630)    RLONG

      WRITE(OR,650)   'ECCENT       '
      WRITE(OR,630)    ECCENT

      WRITE(OR,650)   'PHI0         '
      WRITE(OR,630)    PHI0

      WRITE(OR,650)   'THETAP       '
      WRITE(OR,630)    THETAP

      WRITE(OR,650)   'YTIME        '
      WRITE(OR,630)    YTIME

      WRITE(OR,650)   'DTIME        '
      WRITE(OR,630)    DTIME

      WRITE(OR,650)   'SOLARC       '
      WRITE(OR,630)    SOLARC

      WRITE(OR,650)   'XI0          '
      WRITE(OR,630)    XI0

      WRITE(OR,650)   'TIME0        '
      WRITE(OR,630)    TIME0

      WRITE(OR,650)   'LAT          '
      WRITE(OR,620)    LAT

      WRITE(OR,650)   'TIME         '
      WRITE(OR,620)    TIME

      WRITE(OR,650)   'PHIHARU      '
      WRITE(OR,620)    PHIHARU

      WRITE(OR,650)   'F            '
      WRITE(OR,620)    F

*-----------------------------------------------------------------------

       RETURN
       END
