*"%表題 季節平均, 年平均日射量の計算(地球用) *% *"%履歴 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