*"%表題   季節平均, 年平均日射量の計算(地球用)
*%
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\KISETU.FOR
*"%       91/03/14 佐藤正樹       C:\FORT\NISSHA\KISETU.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  KISETU

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

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

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

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

      REAL  LAT(-NLAT:NLAT)    !"緯度
      REAL  TIME(0:NTIME)
     &                         !"時間
     &                         !"(近日点から測る, YTIMEで規格化)

      REAL  R(0:NTIME)         !"太陽からの距離(RLONGで規格化)
      REAL  SOLAR(0:NTIME)     !"太陽定数
      REAL  DELTAS(0:NTIME)    !"太陽傾斜角
      REAL  PHI(0:NTIME)       !"近日点から測った黄経
      REAL  XI(0:NTIME)        !"楕円をあらわすパラメーター
      REAL  F(0:NTIME, -NLAT:NLAT)  !"日平均日射量
      REAL  H0(0:NTIME, -NLAT:NLAT) !"日照時間



      REAL  TKIN              !"近日点通過の日付 ( 1月5日2h )

      INTEGER  IKIN            !"近日点の時間次元修正値
      INTEGER  MTIME           !"時間次元修正値
      REAL     TTIME           !"時間修正値



      PARAMETER( NTUKI= 12 )           !"月数
      PARAMETER( NKISE= 4 )            !"季節数

      REAL  FTUKI(-NLAT:NLAT, 1:NTUKI)  !"月平均日射量
      REAL  FKISE(-NLAT:NLAT, 1:NKISE)  !"季節平均日射量
      REAL  FNEN(-NLAT:NLAT)            !"年平均日射量
      REAL  H0TUKI(-NLAT:NLAT, 1:NTUKI) !"月平均日照時間
      REAL  H0KISE(-NLAT:NLAT, 1:NKISE) !"季節平均日照時間
      REAL  H0NEN(-NLAT:NLAT)           !"年平均日照時間

      REAL  TTUKI(0:NTUKI)              !"月
      REAL  TKISE(0:NKISE+1)            !"季節

      INTEGER  CTUKI(1:NTUKI)           !"月のデータ数
      INTEGER  CKISE(1:NKISE)           !"季節のデータ数
      INTEGER  CNEN                     !"年のデータ数

      DATA  TTUKI / 0., 31., 59., 90., 120., 151.,
     &            181., 212., 243., 273., 304., 334., 365. /
      DATA  TKISE / 0., 59., 151., 243., 334., 365. /



      INTEGER  ITIME           !"DO ループ
      INTEGER  ILAT            !"DO ループ
      INTEGER  ITUKI           !"DO ループ
      INTEGER  IKISE           !"DO ループ
      INTEGER  LKISE           !"DO ループ



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

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

      INTEGER   IR   / 10 /          !"入力ユニット番号
      CHARACTER HFLDAT*13            !"入力ファイル名
      LOGICAL   LIEX                 !"入力ファイルの存在

      INTEGER   OR   / 20 /          !"出力ユニット番号
      CHARACTER HOFL*13 / 'KISETU.DAT' /  !"出力ファイル名
      LOGICAL   LOEX                  !"出力ファイルの存在


*     < NAMELIST >

      NAMELIST / FLPARA /   HFLDAT

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

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

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

*-----------------------------------------------------------------------
*"    < 入力データファイルのオープン >

      INQUIRE(FILE=HFLDAT, EXIST=LIEX)
      IF( .NOT. LIEX ) THEN
        WRITE(6,*) '入力ファイルが存在しません'
        WRITE(6,*) ' ファイル名: ', HFLDAT
        STOP
      END IF
      
      
      OPEN ( IR, FILE=HFLDAT, STATUS='OLD' )

*"    < 入力 >

      CALL  INISSHA
     O           ( IR,
     O             ISOLAR,   RLONG,    ECCENT,   PHI0,     THETAP,
     O             YTIME,    DTIME,
     O             SOLARC,   XI0,      TIME0,
     O             LAT,      TIME,
     O             R,        SOLAR,    DELTAS,   PHI,      XI,
     O             F,        H0,
     D             NLAT,     NTIME )

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

      CLOSE(IR)

*-----------------------------------------------------------------------
*"    < 初期化 >

      DO 1000 ILAT = -NLAT, NLAT
      
        DO 1010 ITUKI = 1, NTUKI
          FTUKI(ILAT, ITUKI)  = 0.
          H0TUKI(ILAT, ITUKI) = 0.
 1010   CONTINUE

        DO 1020 IKISE = 1, NKISE
          FKISE(ILAT, IKISE)  = 0.
          H0KISE(ILAT, IKISE) = 0.
 1020   CONTINUE

        FNEN(ILAT)  = 0.
        H0NEN(ILAT) = 0.

 1000 CONTINUE



      DO 1030 ITUKI = 1, NTUKI
        CTUKI(ITUKI)  = 0
 1030 CONTINUE

      DO 1040 IKISE = 1, NKISE
        CKISE(IKISE)  = 0
 1040 CONTINUE

      CNEN = 0



*"    < 単位の変更 >

      DO 200 ILAT = -NLAT, NLAT
        LAT(ILAT) = LAT(ILAT) / PI * 180.
 200  CONTINUE



*-----------------------------------------------------------------------
*"    < 時間の調整 :  近日点通過の日付を 1月5日2hとする(理科年表'90) >

      TKIN = ( 5. + 2. / 24 ) / 365.
      IKIN  = NTIME * MOD( TKIN, 1.D0 )



*"    < 平均化 >

      DO 10 ITIME = 0, NTIME
        
        MTIME = MOD( ITIME + IKIN, NTIME )
        TTIME = REAL(MTIME) / REAL(NTIME) * 365.
        
        DO 20 ITUKI = 1, NTUKI
          
          IF(       ( TTIME .GE. TTUKI(ITUKI-1) ) 
     &        .AND. ( TTIME .LT. TTUKI(ITUKI)   ) ) THEN
            
C            WRITE(6,*) 'ITIME', ITIME
C            WRITE(6,*) 'ITUKI', ITUKI
            
            CTUKI(ITUKI) = CTUKI(ITUKI) + 1
            DO 30 ILAT = -NLAT, NLAT
              
              FTUKI(ILAT, ITUKI)  = FTUKI(ILAT, ITUKI)  + F(ITIME, ILAT)
              H0TUKI(ILAT, ITUKI) = H0TUKI(ILAT, ITUKI) 
     &                            + H0(ITIME, ILAT) / PI * DTIME
              
 30         CONTINUE
            
          END IF
        
 20     CONTINUE
 
 
 
        DO 40 IKISE = 1, NKISE+1
          
          IF(       ( TTIME .GE. TKISE(IKISE-1) ) 
     &        .AND. ( TTIME .LT. TKISE(IKISE)   ) ) THEN
            
            LKISE = IKISE
            IF( IKISE .EQ. NKISE+1 ) THEN
              LKISE = 1
            END IF
            
C            WRITE(6,*) 'ITIME', ITIME
C            WRITE(6,*) 'LKISE', LKISE
            
            CKISE(LKISE) = CKISE(LKISE) + 1
            DO 50 ILAT = -NLAT, NLAT
              
              FKISE(ILAT, LKISE)  = FKISE(ILAT, LKISE)  + F(ITIME, ILAT)
              H0KISE(ILAT, LKISE) = H0KISE(ILAT, LKISE) 
     &                            + H0(ITIME, ILAT) / PI * DTIME
              
 50         CONTINUE
            
          END IF
        
 40     CONTINUE
 
 
 
        CNEN = CNEN + 1
        DO 60 ILAT = -NLAT, NLAT
 
          FNEN(ILAT)  = FNEN(ILAT)  + F(ITIME, ILAT)
          H0NEN(ILAT) = H0NEN(ILAT) 
     &                      + H0(ITIME, ILAT) / PI * DTIME
          
 60     CONTINUE
 
 10   CONTINUE



      DO 100 ILAT = -NLAT, NLAT
        
        DO 110 ITUKI = 1, NTUKI
        
          FTUKI(ILAT, ITUKI)  = FTUKI(ILAT, ITUKI)  / CTUKI(ITUKI)
          H0TUKI(ILAT, ITUKI) = H0TUKI(ILAT, ITUKI) / CTUKI(ITUKI)
        
 110    CONTINUE
 
        DO 120 IKISE = 1, NKISE
        
          FKISE(ILAT, IKISE)  = FKISE(ILAT, IKISE)  / CKISE(IKISE)
          H0KISE(ILAT, IKISE) = H0KISE(ILAT, IKISE) / CKISE(IKISE)
        
 120    CONTINUE
        
        FNEN(ILAT)  = FNEN(ILAT)  / CNEN
        H0NEN(ILAT) = H0NEN(ILAT) / CNEN
        
 
 100  CONTINUE

*-----------------------------------------------------------------------
*"    < テスト出力 >

      WRITE(6,*) 'ILAT=',-NLAT,NLAT,20
      WRITE(6,*) '月'
      WRITE(6,*) 'F'
      DO 3000 ITUKI = 1, NTUKI
        WRITE(6,*) 
     &  ITUKI, 
     &  FTUKI(-NLAT,ITUKI), FTUKI( NLAT,ITUKI), FTUKI( 20 ,ITUKI)
 3000 CONTINUE

      WRITE(6,*) 'H0'
      DO 3005 ITUKI = 1, NTUKI
        WRITE(6,*) 
     &  ITUKI, 
     &  H0TUKI(-NLAT,ITUKI), H0TUKI( NLAT,ITUKI), H0TUKI( 20 ,ITUKI)
 3005 CONTINUE

      WRITE(6,*) '季節'
      WRITE(6,*) 'F'
      DO 3010 IKISE = 1, NKISE
        WRITE(6,*) 
     &  IKISE, 
     &  FKISE(-NLAT,IKISE), FKISE( NLAT,IKISE), FKISE( 20 ,IKISE)
 3010 CONTINUE

      WRITE(6,*) 'H0'
      DO 3015 IKISE = 1, NKISE
        WRITE(6,*) 
     &  IKISE, 
     &  H0KISE(-NLAT,IKISE), H0KISE( NLAT,IKISE), H0KISE( 20 ,IKISE)
 3015 CONTINUE

      WRITE(6,*) '年'
      WRITE(6,*) 'F'
      WRITE(6,*) 
     &  FNEN(-NLAT), FNEN( NLAT), FNEN( 20 )
      WRITE(6,*) 'H0'
      WRITE(6,*) 
     &  H0NEN(-NLAT), H0NEN( NLAT), H0NEN( 20 )

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

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

      WRITE(6,*) '出力します : ファイル', HOFL
      OPEN ( OR, FILE=HOFL )

*"    < 出力 >

      CALL  OHEKIN
     I           ( OR,
     I             ISOLAR,   RLONG,    ECCENT,   PHI0,     THETAP,
     I             YTIME,    DTIME,
     I             LAT,      TIME,
     I             TTUKI,    TKISE, 
     I             FTUKI,    FKISE,    FNEN, 
     I             H0TUKI,   H0KISE,   H0NEN, 
     D             NLAT,     NTIME,    NTUKI,   NKISE )

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

      CLOSE(OR)

*-----------------------------------------------------------------------
      STOP
      END
*=======================================================================
*"%表題   日射の計算値の入力
*%
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHAZU.FOR
*%
*%
*-----------------------------------------------------------------------
      SUBROUTINE  INISSHA
     O           ( IR,
     O             ISOLAR,   RLONG,    ECCENT,   PHI0,     THETAP,
     O             YTIME,    DTIME,
     O             SOLARC,   XI0,      TIME0,
     O             LAT,      TIME,
     O             R,        SOLAR,    DELTAS,   PHI,      XI,
     O             F,        H0,
     D             NLAT,     NTIME )

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

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

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

      REAL  LAT(-NLAT:NLAT)    !"緯度
      REAL  TIME(0:NTIME)
     &                         !"時間
     &                         !"(近日点から測る, YTIMEで規格化)

      REAL  R(0:NTIME)         !"太陽からの距離(RLONGで規格化)
      REAL  SOLAR(0:NTIME)     !"太陽定数
      REAL  DELTAS(0:NTIME)    !"太陽傾斜角
      REAL  PHI(0:NTIME)       !"近日点から測った黄経
      REAL  XI(0:NTIME)        !"楕円をあらわすパラメーター
      REAL  F(0:NTIME, -NLAT:NLAT)  !"日平均日射量
      REAL  H0(0:NTIME, -NLAT:NLAT) !"日照時間


      CHARACTER*13     HDUMMY  !"文字入力のダミー

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

 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  )

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

      READ(IR,650)     HDUMMY
      READ(IR,630)     ISOLAR

      READ(IR,650)     HDUMMY
      READ(IR,630)     RLONG

      READ(IR,650)     HDUMMY
      READ(IR,630)     ECCENT

      READ(IR,650)     HDUMMY
      READ(IR,630)     PHI0

      READ(IR,650)     HDUMMY
      READ(IR,630)     THETAP

      READ(IR,650)     HDUMMY
      READ(IR,630)     YTIME

      READ(IR,650)     HDUMMY
      READ(IR,630)     DTIME

      READ(IR,650)     HDUMMY
      READ(IR,630)     SOLARC

      READ(IR,650)     HDUMMY
      READ(IR,630)     XI0

      READ(IR,650)     HDUMMY
      READ(IR,630)     TIME0

      READ(IR,650)     HDUMMY
      READ(IR,620)     LAT

      READ(IR,650)     HDUMMY
      READ(IR,620)     TIME

      READ(IR,650)     HDUMMY
      READ(IR,620)     R

      READ(IR,650)     HDUMMY
      READ(IR,620)     SOLAR

      READ(IR,650)     HDUMMY
      READ(IR,620)     DELTAS

      READ(IR,650)     HDUMMY
      READ(IR,620)     PHI

      READ(IR,650)     HDUMMY
      READ(IR,620)     XI

      READ(IR,650)     HDUMMY
      READ(IR,620)     F

      READ(IR,650)     HDUMMY
      READ(IR,620)     H0

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

       RETURN
       END
*=======================================================================
*"%表題   日射の平均値の計算値の出力
*%
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHAZU.FOR
*%
*%
*-----------------------------------------------------------------------
      SUBROUTINE  OHEKIN
     I           ( OR,
     I             ISOLAR,   RLONG,    ECCENT,   PHI0,     THETAP,
     I             YTIME,    DTIME,
     I             LAT,      TIME,
     I             TTUKI,    TKISE, 
     I             FTUKI,    FKISE,    FNEN, 
     I             H0TUKI,   H0KISE,   H0NEN, 
     D             NLAT,     NTIME,    NTUKI,   NKISE )

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

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

      REAL  LAT(-NLAT:NLAT)    !"緯度
      REAL  TIME(0:NTIME)
     &                         !"時間
     &                         !"(近日点から測る, YTIMEで規格化)

      REAL  TTUKI(0:NTUKI)              !"月
      REAL  TKISE(0:NKISE+1)            !"季節

      REAL  FTUKI(-NLAT:NLAT, 1:NTUKI)  !"月平均日射量
      REAL  FKISE(-NLAT:NLAT, 1:NKISE)  !"季節平均日射量
      REAL  FNEN(-NLAT:NLAT)            !"年平均日射量
      REAL  H0TUKI(-NLAT:NLAT, 1:NTUKI) !"月平均日照時間
      REAL  H0KISE(-NLAT:NLAT, 1:NKISE) !"季節平均日照時間
      REAL  H0NEN(-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)   'LAT          '
      WRITE(OR,620)    LAT

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

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

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

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

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

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

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

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

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

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

       RETURN
       END
