*"%表題   日射量, 日照時間の時間 - 緯度面での描画(地球用)
*%
*"%履歴   91/02/26 佐藤正樹       C:\FORT\NISSHA\NISSHAZU.FOR
*"%       91/03/14 佐藤正樹       C:\FORT\NISSHA\TIKYUZU.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  TIKYUZU

      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)    !"緯度
      DOUBLE PRECISION  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 )
      REAL  F2(0:NTIME, -NLAT:NLAT)  !"日平均日射量(修正後)
      REAL  H02(0:NTIME, -NLAT:NLAT) !"日照時間(修正後)

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



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

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

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



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

      PARAMETER( N1 = 0, N2 = 13, NC = 2 )
      CHARACTER*(NC)  CH(N2)   !"X軸のラベル
C      REAL  UX1(N1)
      REAL  UX2(N2)            !"X軸の目盛の位置

      DATA  CH  / '1', '2', '3', '4', '5', '6',
     &            '7', '8', '9', '10', '11', '12', ' ' /
      DATA  UX2 / 0., 31., 59., 90., 120., 151.,
     &            181., 212., 243., 273., 304., 334., 365. /

      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 = 1000.
      DF   = 50.
      LH0SC = .FALSE.
      H0MIN = 0.
      H0MAX = 24.
      DH0   = 2.



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

*-----------------------------------------------------------------------
*"    < 時間の調整 :  近日点通過の日付を 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 )
C        WRITE(6,*) 'MTIME', MTIME
        
        DO 20 ILAT = -NLAT, NLAT
          
          F2(MTIME, ILAT)  = F(ITIME, ILAT)
          H02(MTIME, ILAT) = H0(ITIME, ILAT) / PI * DTIME
          
 20     CONTINUE
 10   CONTINUE


      DO 30 ILAT = -NLAT, NLAT
        
        F2(NTIME, ILAT)  = F2(0, ILAT)
        H02(NTIME, ILAT) = H02(0, ILAT)
        
 30   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 USCSET('CXSIDE', ' ')         !"X軸を描かない
      CALL USPSET('XMIN', 0.)
      CALL USPSET('XMAX', 365.)

      CALL USAXIS

      CALL UZPSET('LBTWN', .TRUE.)       !"目盛の間にラベルを書く
      CALL UXAXLB('B', UX2, N1, UX2, CH, NC, N2) !"X軸下側
      CALL UXAXNM('T', UX2, N1, UX2, N2)         !"X軸上側
      CALL UXPTTL('B', 1, 'manth', 0)            !"X軸タイトル

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

      CALL UDCNTR(F2, 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,  'hour')
*                                12345678901234567890


      CALL USAXIS
      CALL UZPSET('LBTWN', .TRUE.)       !"目盛の間にラベルを書く
      CALL UXAXLB('B', UX2, N1, UX2, CH, NC, N2) !"X軸下側
      CALL UXAXNM('T', UX2, N1, UX2, N2)         !"X軸上側
      CALL UXPTTL('B', 1, 'manth', 0)            !"X軸タイトル

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

      CALL UDCNTR(H02, 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
