*"%表題 日射量, 日照時間の一定期間平均(一般惑星用) *"% 春分の日を基準に平均をとる *% *"%履歴 91/03/08 佐藤正樹 C:\FORT\NISSHA\HEKIN2.FOR *"% 91/03/14 佐藤正樹 C:\FORT\NISSHA\HEKIN.FOR *% *% *----------------------------------------------------------------------- PROGRAM HEKIN 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 !"平均数 REAL PHIHARU !"春分の日の黄経 REAL DTH !"平均時間間隔 REAL TH(1:NTIME) !"平均をとる時間 REAL FH(-NLAT:NLAT, 1:NTIME) !"平均日射量 REAL H0H(-NLAT:NLAT, 1:NTIME) !"平均日照時間 REAL FNEN(-NLAT:NLAT) !"年平均日射量 REAL H0NEN(-NLAT:NLAT) !"年平均日照時間 INTEGER CH(0:NTIME) !"平均をとるデータ数 INTEGER CNEN !"年のデータ数 INTEGER ILAT !"DO ループで用いる緯度配列の次元 INTEGER ITIME !"DO ループで用いる時間配列の次元 INTEGER INUM !"DO ループで用いる平均数 INTEGER LNUM !"平均数として一時的に用いる *" < 入出力パラメーター > 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 1010 ITIME = 1, NTIME DO 1020 ILAT = -NLAT, NLAT FH(ILAT, ITIME) = 0. H0H(ILAT, ITIME) = 0. 1020 CONTINUE CH(ITIME) = 0 1010 CONTINUE 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 *----------------------------------------------------------------------- *" < 平均をとる時間の設定 > PHIHARU = MOD( 3. * PI - PHI0, 2. * PI ) DO 300 ITIME = 0, NTIME-1 IF( ( PHI(ITIME) .GE. PHIHARU ) & .AND. ( PHI(ITIME+1) .GE. PHIHARU ) ) THEN TH(1) = TIME(ITIME) & + ( PHIHARU - PHI(ITIME) ) & / ( PHI(ITIME+1) - PHI(ITIME) ) & * ( TIME(ITIME+1) - TIME(ITIME) ) GOTO 305 END IF 300 CONTINUE 305 CONTINUE DTH = TIME(NTIME) / REAL( NUMBER ) DO 310 INUM = 2, NUMBER TH(INUM) = MOD( TH(1) + DTH * REAL( INUM - 1 ), TIME(NTIME) ) 310 CONTINUE *" < テスト出力1 > WRITE(6,*) '平均数 : NUMBER=', NUMBER WRITE(6,*) '春分の日の黄経 : PHIHARU', PHIHARU WRITE(6,*) '平均時間間隔 : DTH', DTH WRITE(6,*) 'TH' WRITE(6, 600) (TH(INUM), INUM=1,NUMBER) 600 FORMAT(2X, 5F13.4) *----------------------------------------------------------------------- *" < 平均化 > DO 10 ITIME = 0, NTIME DO 20 INUM = 1, NUMBER IF( ( ABS( TIME(ITIME) - TH(INUM) ) .LT. DTH / 2. ) & .OR. ( ABS( TIME(ITIME) - TH(INUM) + TIME(NTIME) ) & .LT. DTH / 2. ) & .OR. ( ABS( TIME(ITIME) - TH(INUM) - TIME(NTIME) ) & .LT. DTH / 2. ) & ) THEN LNUM = INUM GOTO 25 END IF 20 CONTINUE 25 CONTINUE C WRITE(6,*) 'ITIME', ITIME C WRITE(6,*) 'LNUM', LNUM CH(LNUM) = CH(LNUM) + 1 DO 30 ILAT = -NLAT, NLAT FH(ILAT, LNUM) = FH(ILAT, LNUM) + F(ITIME, ILAT) H0H(ILAT, LNUM) = H0H(ILAT, LNUM) + H0(ITIME, ILAT) 30 CONTINUE CNEN = CNEN + 1 DO 40 ILAT = -NLAT, NLAT FNEN(ILAT) = FNEN(ILAT) + F(ITIME, ILAT) H0NEN(ILAT) = H0NEN(ILAT) + H0(ITIME, ILAT) 40 CONTINUE 10 CONTINUE DO 50 ILAT = -NLAT, NLAT DO 60 INUM = 1, NUMBER FH(ILAT, INUM) = FH(ILAT, INUM) / REAL( CH(INUM) ) H0H(ILAT, INUM) = H0H(ILAT, INUM) / REAL( CH(INUM) ) 60 CONTINUE FNEN(ILAT) = FNEN(ILAT) / REAL( CNEN ) H0NEN(ILAT) = H0NEN(ILAT) / REAL( CNEN ) 50 CONTINUE *----------------------------------------------------------------------- *" < テスト出力2 > WRITE(6,*) 'ILAT=',-NLAT,NLAT,20 WRITE(6,*) '平均' WRITE(6,*) 'F' DO 3000 INUM = 1, NUMBER WRITE(6,*) & INUM, & FH(-NLAT,INUM), FH( NLAT,INUM), FH( 20 ,INUM) 3000 CONTINUE WRITE(6,*) 'H0' DO 3005 INUM = 1, NUMBER WRITE(6,*) & INUM, & H0H(-NLAT,INUM), H0H( NLAT,INUM), H0H( 20 ,INUM) 3005 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 ) *----------------------------------------------------------------------- *" < 描画 > 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 INUM = 1, NUMBER CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), FH(-NLAT, INUM), 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', 'hour') CALL USAXIS DO 2100 INUM = 1, NUMBER CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), H0H(-NLAT, INUM), 1, 3) 2100 CONTINUE CALL SGPLZU( 2*NLAT+1, LAT(-NLAT), H0NEN(-NLAT), 1, 5) 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