*"%表題   日平均日射量の作画
*%
*"%履歴   91/02/27 佐藤正樹       C:\FORT\NISSHA\NISSHAZ2.FOR
*"%       91/03/07 佐藤正樹       C:\FORT\NISSHA\NISSHAZ2.FOR
*"%       91/03/14 佐藤正樹       C:\FORT\NISSHA\NISSHAZU.FOR
*"%       91/03/15 佐藤正樹       C:\FORT\NISSHA\NISSHAZU.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  NISSHAZU

      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) !"日照時間



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

      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

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

      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.



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

*-----------------------------------------------------------------------
*"    < 単位の調整 >

      DO 10 ITIME = 0, NTIME
        DO 20 ILAT = -NLAT, NLAT
          H0(ITIME, ILAT) = H0(ITIME, ILAT) / PI * DTIME
 20     CONTINUE
 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 SGSTXC( 1 )
      CALL SGSTXS( 0.02 )
      CALL SGSTXI( 3 )
      CALL SGTXV(   0.8,  0.9,  'W/m|2"')
*                                12345678901234567890


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

      CALL USPSET('YMIN', -90.)
      CALL USPSET('YMAX',  90.)
      CALL USPSET('DYL', 30.)
      CALL USCSET('CYTTL' , 'latitude')

      CALL USPSET('XMIN', 0.)
      CALL USPSET('XMAX', 1.)
      CALL USPSET('DXL', 0.2)
      CALL USCSET('CXTTL' , 'time / year')

      CALL USAXIS

      IF( LFSC ) THEN
        CALL UDGCLA(FMIN, FMAX, DF)        !"コンター間隔
      END IF

      CALL UDCNTR(F, NTIME+1, NTIME+1, 2*NLAT+1)



      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


      CALL SGSTXC( 1 )
      CALL SGSTXS( 0.02 )
      CALL SGSTXI( 3 )
      CALL SGTXV(   0.8,  0.9,  'day')
*                                12345678901234567890


      CALL USAXIS

      IF( LH0SC ) THEN
        CALL UDGCLA(H0MIN, H0MAX, DH0)        !"コンター間隔
      END IF

      CALL UDCNTR(H0, NTIME+1, NTIME+1, 2*NLAT+1)


      CALL GRCLS

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