*"%表題   日射量, 日照時間の一定時間間隔, 年平均値の
*"%       緯度分布のグラフ(一般惑星用)
*%
*"%履歴   91/02/28 佐藤正樹       C:\FORT\NISSHA\NISSHAG.FOR
*"%       91/03/07 佐藤正樹       C:\FORT\NISSHA\NISSHAG.FOR
*"%       91/03/14 佐藤正樹       C:\FORT\NISSHA\NISSHAG.FOR
*"%       91/03/15 佐藤正樹       C:\FORT\NISSHA\NISSHAG.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  NISSHAG

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

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

      REAL  ISOLAR  !"太陽の全日射量
      REAL  RLONG   !"軌道長半径
      REAL  ECCENT  !"離心率
      REAL  PHI0    !"近日点黄経(rad)
      REAL  THETAP  !"赤道傾斜角(rad)
      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) !"日照時間



      INTEGER  NUMBER                   !"グラフ出力数
      INTEGER  LTIME                    !"作画時間次元間隔
      REAL  FG(-NLAT:NLAT)              !"グラフ出力用日射量
      REAL  H0G(-NLAT:NLAT)             !"グラフ出力用日照時間
      REAL  FNEN(-NLAT:NLAT)            !"年平均日射量
      REAL  H0NEN(-NLAT:NLAT)           !"年平均日照時間



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

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

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



*"    < 描画パラメーター >

      LOGICAL  LFSC                 !"F 軸のスケーリングを外部指定する
      LOGICAL  LH0SC                !"H0軸のスケーリングを外部指定する
      REAL  FMIN,     FMAX,     DF  !"F 軸のスケーリングパラメーター
      REAL  H0MIN,    H0MAX,    DH0 !"H0軸のスケーリングパラメーター



*     < NAMELIST >

      NAMELIST / FLPARA /   HFLDAT
      NAMELIST / FHPARA /   LFSC,     FMIN,     FMAX,     DF, 
     &                      LH0SC,    H0MIN,    H0MAX,    DH0
      NAMELIST / NUMPARA /  NUMBER

*-----------------------------------------------------------------------
*"    < 入力データファイル名の入力 >

      HFLDAT = 'NISSHA.DAT'

      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)

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

      LFSC = .FALSE.
      FMIN = 0.
      FMAX = ISOLAR / ( 4.D0 * PI * RLONG ** 2 )
      DF   = FMAX / 10.
      LH0SC = .FALSE.
      H0MIN = 0.
      H0MAX = DTIME
      DH0   = DTIME / 10.

      NUMBER = 1



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


      OPEN(IPR,FILE=HIPFL)
        READ(IPR,FHPARA)
        READ(IPR,NUMPARA)
      CLOSE(IPR)


      IF( NUMBER .GT. NTIME ) THEN
        WRITE(6,*) '出力する数が大きすぎます : NUMBER=', NUMBER
        STOP
      END IF

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

      DO 1000 ILAT = -NLAT, NLAT
      
        FNEN(ILAT)  = 0.
        H0NEN(ILAT) = 0.

 1000 CONTINUE

      CNEN = 0



*"    < 単位の変更 >

      DO 100 ILAT = -NLAT, NLAT
        LAT(ILAT) = LAT(ILAT) / PI * 180.
        DO 200 ITIME = 0, NTIME
          H0(ITIME, ILAT) = H0(ITIME, ILAT) / PI * DTIME
 200    CONTINUE
 100  CONTINUE

*-----------------------------------------------------------------------
*"    < 出力間隔の設定 >

      LTIME = INT( REAL( NTIME ) / REAL( NUMBER ) )



*"    < 平均化 >

      DO 10 ILAT = -NLAT, NLAT
        DO 20 ITIME = 0, NTIME-1

          FNEN(ILAT)  = FNEN(ILAT)  + F(ITIME, ILAT)
          H0NEN(ILAT) = H0NEN(ILAT) + H0(ITIME, ILAT)

 20     CONTINUE

        FNEN(ILAT)  = FNEN(ILAT)  / REAL( NTIME )
        H0NEN(ILAT) = H0NEN(ILAT) / REAL( NTIME )

 10   CONTINUE

*-----------------------------------------------------------------------
*"    < 描画 >

      WRITE(6,*) 'workstation no.? '
      READ (5,*) IWS

      CALL GROPN(IWS)


      CALL GRFRM
      CALL USSVPT( 0.15, 0.8, 0.2, 0.85 )

      CALL SGSTXC( 0 )
      CALL SGSTXS( 0.03 )
      CALL SGSTXI( 5 )
      CALL SGTXV(   0.5,  0.95,  'Solar Radiation')
*                                 12345678901234567890


      CALL UZPSET('INNER', -1)           !"外側に向かって目盛を打つ

      CALL USPSET('XMIN', -90.)
      CALL USPSET('XMAX',  90.)
      CALL USPSET('DXL', 30.)
      CALL USCSET('CXTTL' , 'latitude')



      IF( LFSC ) THEN
        CALL USPSET('YMIN', FMIN )
        CALL USPSET('YMAX', FMAX )
        CALL USPSET('DYL', DF )
      ELSE
        CALL USSYEX((2*NLAT+1)*(NTIME+1), F)
      END IF

      CALL USCSET('CYTTL' , 'solar radiation')
      CALL USCSET('CYUNIT', 'W/m|2"')

      CALL USAXIS

      DO 2000 ITIME = 0, NTIME, LTIME
        DO 2010 ILAT = -NLAT, NLAT
          FG(ILAT)=F(ITIME, ILAT)
 2010   CONTINUE
        CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), FG(-NLAT), 1, 3)
 2000 CONTINUE
      CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), FNEN(-NLAT), 1, 5)



      CALL GRFRM
      CALL USSVPT( 0.15, 0.8, 0.2, 0.85 )

      CALL SGSTXC( 0 )
      CALL SGSTXS( 0.03 )
      CALL SGSTXI( 5 )
      CALL SGTXV(   0.5,  0.95,  'Duration of Sunshine')
*                                 12345678901234567890



      IF( LH0SC ) THEN
        CALL USPSET('YMIN', H0MIN )
        CALL USPSET('YMAX', H0MAX )
        CALL USPSET('DYL', DH0 )
      ELSE
        CALL USSYEX((2*NLAT+1)*(NTIME+1), H0)
      END IF

      CALL USCSET('CYTTL' , 'duration of sunshine')
      CALL USCSET('CYUNIT', 'day')

      CALL USAXIS

      DO 2100 ITIME = 0, NTIME, LTIME
        DO 2110 ILAT = -NLAT, NLAT
          H0G(ILAT)=H0(ITIME, ILAT)
 2110   CONTINUE
        CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), H0G(-NLAT), 1, 3)
 2100 CONTINUE
      CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), H0NEN(-NLAT), 1, 5)



      CALL GRCLS

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

C      WRITE(6,*) 'FNEN', FNEN

*-----------------------------------------------------------------------
      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
