*"%表題 日射の計算(春分点を基準に日射量だけ計算する) *% *"%履歴 91/02/22 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *"% 91/03/14 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *"% 91/03/21 佐藤正樹 C:\FORT\NISSHA\NISSHA2.FOR *"% F は倍精度だと領域オーバーになってしまうので *"% 単精度に変更した *----------------------------------------------------------------------- PROGRAM NISSHA2 PARAMETER ( NLAT = 90 ) !"緯度配列の次元 PARAMETER ( NTIME = 60 ) !"時間配列の次元 DOUBLE PRECISION PI PARAMETER( PI = 3.141592653589793 D 0 ) !"円周率 DOUBLE PRECISION ISOLAR !"太陽の全日射量 DOUBLE PRECISION RLONG !"軌道長半径 DOUBLE PRECISION ECCENT !"離心率 DOUBLE PRECISION PHI0 !"近日点黄経(rad) DOUBLE PRECISION THETAP !"赤道傾斜角(rad) DOUBLE PRECISION YTIME !"公転周期 DOUBLE PRECISION DTIME !"自転周期 DOUBLE PRECISION PHI0A !"近日点黄経(度) DOUBLE PRECISION THETAPA !"赤道傾斜角(度) DOUBLE PRECISION SOLARC !"遠日点(R=RLONG)での太陽定数 DOUBLE PRECISION XI0 !"XI の春分点での値 * (-1) DOUBLE PRECISION TIME0 !"春分点での時間 * (-1) DOUBLE PRECISION LAT(-NLAT:NLAT) !"緯度 DOUBLE PRECISION TIME(0:NTIME) & !"時間(春分点から測る, YTIMEで規格化) DOUBLE PRECISION R !"太陽からの距離(RLONGで規格化) DOUBLE PRECISION SOLAR !"太陽定数 DOUBLE PRECISION DELTAS !"太陽傾斜角 DOUBLE PRECISION PHI !"近日点から測った黄経 DOUBLE PRECISION XI !"楕円をあらわすパラメーター DOUBLE PRECISION H0 !"日照時間 DOUBLE PRECISION PHIHARU(0:NTIME) !"春分点から測った黄経 C" DOUBLE PRECISION F(0:NTIME, -NLAT:NLAT) ! 日平均日射量 REAL F(0:NTIME, -NLAT:NLAT) !"日平均日射量 DOUBLE PRECISION FTIME !"関数副プログラム DOUBLE PRECISION FXI !"関数副プログラム INTEGER ILAT !"DO ループで用いる緯度配列の次元 INTEGER ITIME !"DO ループで用いる時間配列の次元 DOUBLE PRECISION TAIYOU !"一時変数, 太陽の出入り *" < 入出力パラメーター > INTEGER IPR / 9 / !"入力ユニット番号 CHARACTER HIPFL*13 / 'NISSHA.PAR' / & !"入力パラメーターファイル名 LOGICAL LIPEX !"入力パラメーターファイルの存在 INTEGER OR / 20 / !"出力ユニット番号 CHARACTER HFLDAT*13 !"出力ファイル名 LOGICAL LOEX !"出力ファイルの存在 * < NAMELIST > NAMELIST / FLPARA / HFLDAT NAMELIST / IPARA / ISOLAR, RLONG, ECCENT, & PHI0A, THETAPA, & YTIME, DTIME *----------------------------------------------------------------------- *" < 外部パラメーター > *" << 地球のデータ >> HFLDAT = 'NISSHA.DAT' ISOLAR = 3.85 D 26 ! [W] RLONG = 1.496 D 11 ! [m] ECCENT = 0.0167 D 0 PHI0A = 102.908 D 0 !"[度] THETAPA = 23.44 D 0 !"[度] YTIME = 365.0 D 0 !"[日] DTIME = 24.0 D 0 !"[時間] *" < パラメーターの入力 > INQUIRE(FILE=HIPFL, EXIST=LIPEX) IF( .NOT. LIPEX ) THEN WRITE(6,*) '入力パラメーターファイルが存在しません' WRITE(6,*) ' ファイル名: ', HIPFL STOP END IF OPEN(IPR,FILE=HIPFL) READ(IPR,FLPARA) READ(IPR,IPARA) CLOSE(IPR) WRITE(6,*) '外部パラメーター' WRITE(6,*) 'ISOLAR', ISOLAR WRITE(6,*) 'RLONG', RLONG WRITE(6,*) 'ECCENT', ECCENT WRITE(6,*) 'PHI0A', PHI0A WRITE(6,*) 'THETAPA', THETAPA WRITE(6,*) 'YTIME', YTIME WRITE(6,*) 'DTIME', DTIME PHI0 = PHI0A / 180.D0 * PI THETAP = THETAPA / 180.D0 * PI *----------------------------------------------------------------------- *" < 定数値の計算 > SOLARC = ISOLAR / ( 4.D0 * PI * RLONG ** 2 ) XI0 = 2.D0 * MOD( & ATAN( SQRT( ( 1.D0 - ECCENT ) / ( 1.D0 + ECCENT ) ) & * TAN( PHI0 / 2.D0 ) & ) + PI, & PI & ) TIME0 = FTIME( XI0, ECCENT ) C WRITE(6,*) 'SOLARC(W/m2)', SOLARC C WRITE(6,*) 'XI0(rad)', XI0 C" WRITE(6,*) 'XI0(度)', XI0 * 180.D0 / PI C WRITE(6,*) 'TIME0', TIME0 C" WRITE(6,*) 'TIME0(日)', TIME0 * YTIME *----------------------------------------------------------------------- *" < 緯度の設定 > LAT(0) = 0.D0 DO 10 ILAT = 1, NLAT LAT(ILAT) = PI/2.D0 * REAL( ILAT ) / REAL( NLAT ) LAT(-ILAT) = - LAT(ILAT) 10 CONTINUE *" < 時間の設定 : 春分点を基準にする > DO 20 ITIME = 0, NTIME TIME(ITIME) = - TIME0 + 1.D0 * REAL( ITIME ) / REAL( NTIME ) TIME(ITIME) = MOD( TIME(ITIME) + 1.D0, 1.D0 ) 20 CONTINUE *----------------------------------------------------------------------- *" < 計算> DO 100 ITIME = 0, NTIME WRITE(6,*) 'ITIME', ITIME XI = FXI( TIME(ITIME), ECCENT ) PHI= 2.D0 * ATAN( & SQRT( ( 1.D0 + ECCENT ) / ( 1.D0 - ECCENT ) ) & * TAN( XI / 2.D0 ) & ) R = 1.D0 - ECCENT * COS( XI ) SOLAR = SOLARC / R ** 2 DELTAS = ASIN( - SIN( THETAP ) & * SIN( PHI0 + PHI ) & ) DO 200 ILAT = -NLAT, NLAT C WRITE(6,*) 'ILAT', ILAT TAIYOU = TAN( LAT(ILAT) ) * TAN( DELTAS ) IF( ABS( TAIYOU ) .LT. 1.D0 ) THEN H0 = ACOS( - TAIYOU ) ELSE IF( TAIYOU .LE. -1.D0 ) THEN H0 = 0.D0 ELSE IF( TAIYOU .GE. 1.D0 ) THEN H0 = PI ENDIF F(ITIME, ILAT) & = SOLAR / PI & * ( SIN( H0 ) & * COS( LAT(ILAT) ) * COS( DELTAS ) & + H0 & * SIN( LAT(ILAT) ) * SIN( DELTAS ) & ) C WRITE(6,'('' '', A, 5F13.4)') C" & 'F(0:NTIME, -NLAT:NLAT) ! 日平均日射量', C & F(ITIME, ILAT) PHIHARU(ITIME) = PHI + PHI0 PHIHARU(ITIME) = MOD( PHIHARU(ITIME) + 2.D0 * PI, 2.D0 * PI ) 200 CONTINUE 100 CONTINUE *----------------------------------------------------------------------- *" < 出力ファイルのオープン > INQUIRE(FILE=HFLDAT, EXIST=LOEX) IF( LOEX ) THEN OPEN ( OR, FILE=HFLDAT, STATUS='OLD' ) CLOSE( OR, STATUS='DELETE' ) END IF OPEN ( OR, FILE=HFLDAT ) *" < 出力 > CALL ONISSHA2 I ( OR, I ISOLAR, RLONG, ECCENT, PHI0, THETAP, I YTIME, DTIME, I SOLARC, XI0, TIME0, I LAT, TIME, PHIHARU, F, D NLAT, NTIME ) *" < ファイルのクローズ > CLOSE(OR) *----------------------------------------------------------------------- STOP END *======================================================================= *"%表題 軌道関数 時間 *% *"%履歴 91/02/23 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *"%履歴 91/02/26 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *% *% *----------------------------------------------------------------------- FUNCTION FTIME( XI, ECCENT ) DOUBLE PRECISION FTIME !"時間 DOUBLE PRECISION XI !"パラメーター DOUBLE PRECISION ECCENT !"離心率 DOUBLE PRECISION PI PARAMETER( PI = 3.141592653589793 D 0 ) !"円周率 *----------------------------------------------------------------------- FTIME = 1.D0 / 2.D0 / PI * ( XI - ECCENT * SIN( XI ) ) C WRITE(6,*) 'SUB : FTIME' C WRITE(6,*) 'XI, XI0, ECCENT, FTIME' C WRITE(6,*) XI, XI0, ECCENT, FTIME *----------------------------------------------------------------------- RETURN END *======================================================================= *"%表題 XI の計算 *% *"%履歴 91/02/22 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *"%履歴 91/02/26 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *% *% *----------------------------------------------------------------------- FUNCTION FXI( TIME, ECCENT ) DOUBLE PRECISION FXI !"楕円をあらわすパラメーター DOUBLE PRECISION TIME !"時間 DOUBLE PRECISION ECCENT !"離心率 DOUBLE PRECISION XI1 !"XI の値の初期値1 DOUBLE PRECISION XI2 !"XI の値の初期値2 DOUBLE PRECISION FTIME !"関数副プログラム : 時間 EXTERNAL FTIME DOUBLE PRECISION FGYAKU !"関数副プログラム : 逆関数 DOUBLE PRECISION PI PARAMETER( PI = 3.141592653589793 D 0 ) !"円周率 DOUBLE PRECISION ATIME !"標準化した時間 *----------------------------------------------------------------------- *" < TIME の標準化 > ATIME = MOD( MOD( TIME, 1.D0 ) + 1.D0, 1.D0 ) *" < XI の初期値の設定 > XI1 = 0.D0 XI2 = 2.D0 * PI *" < XI を求める > FXI = FGYAKU( FTIME, ATIME, XI1, XI2, ECCENT ) C WRITE(6,*) 'FXI, ATIME, XI1, XI2, ECCENT' C WRITE(6,*) FXI, ATIME, XI1, XI2, ECCENT *----------------------------------------------------------------------- RETURN END *======================================================================= *"%表題 逆関数 *% *"%履歴 91/02/23 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *"%履歴 91/02/26 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *% *% *" < X1, X2 との間の2分法によって X を求める > *----------------------------------------------------------------------- FUNCTION FGYAKU( FUNC, Y, X1, X2, PARA ) DOUBLE PRECISION FGYAKU !"逆関数 DOUBLE PRECISION FUNC !"関数 EXTERNAL FUNC DOUBLE PRECISION Y !"関数値 DOUBLE PRECISION X1 !"逆関数の値の初期値1 DOUBLE PRECISION X2 !"逆関数の値の初期値2 DOUBLE PRECISION PARA !"パラメーター *" < SUBROUTINE 内で定義される変数 > DOUBLE PRECISION Y0 !"Y の計算値 DOUBLE PRECISION XX1 ! X1 DOUBLE PRECISION XX2 ! X2 DOUBLE PRECISION XX3 ! X3 DOUBLE PRECISION SGN !"X1, X2 での関数値の差の符号 INTEGER I !"DO ループ INTEGER IMAX / 30 / !"I の最大値 *----------------------------------------------------------------------- *" < 範囲 X1, X2 の適合性のチェック > *" Y - FUNC(X1) と Y - FUNC(X2) の符号が逆である *" ことが仮定されている IF( ( Y - FUNC( X1, PARA ) ) & * ( Y - FUNC( X2, PARA ) ) & .GT. 0.0D0 ) THEN WRITE(6,*) 'SUBROUTINE FGYAKU エラー' WRITE(6,*) '範囲 X1, X2が適切でありません. ' STOP END IF *" < X1, X2 での値の大小関係を調べる > IF( FUNC( X1, PARA ) - FUNC( X2, PARA ) & .LE. 0.0D0 ) THEN SGN = 1.D0 ELSE SGN = - 1.D0 ENDIF *" < くりかえし > XX1 = X1 XX2 = X2 DO 10 I = 1, IMAX XX3 = ( XX1 + XX2 ) / 2.0 Y0 = FUNC( XX3, PARA ) C WRITE(6,*) 'I=', I, 'XX1=', XX1, 'XX2=', XX2 C WRITE(6,*) 'Y0=', Y0, 'Y-Y0=', Y - Y0 IF( ( Y - Y0 ) * SGN .LE. 0.0D0 ) THEN XX2 = XX3 ELSE XX1 = XX3 ENDIF 10 CONTINUE FGYAKU = XX3 *----------------------------------------------------------------------- RETURN END *======================================================================= *"%表題 日射の計算値の出力(簡易) *% *"%履歴 91/02/26 佐藤正樹 C:\FORT\NISSHA\NISSHA.FOR *"%履歴 91/03/21 佐藤正樹 C:\FORT\NISSHA\NISSHA2.FOR *% *% *----------------------------------------------------------------------- SUBROUTINE ONISSHA2 I ( OR, I ISOLAR, RLONG, ECCENT, PHI0, THETAP, I YTIME, DTIME, I SOLARC, XI0, TIME0, I LAT, TIME, PHIHARU, F, D NLAT, NTIME ) INTEGER NLAT !"緯度配列の次元 INTEGER NTIME !"時間配列の次元 INTEGER OR !"出力ユニット番号 DOUBLE PRECISION ISOLAR !"太陽の全日射量 DOUBLE PRECISION RLONG !"軌道長半径 DOUBLE PRECISION ECCENT !"離心率 DOUBLE PRECISION PHI0 !"近日点黄経 DOUBLE PRECISION THETAP !"赤道傾斜角 DOUBLE PRECISION YTIME !"公転周期 DOUBLE PRECISION DTIME !"自転周期 DOUBLE PRECISION SOLARC !"遠日点(R=RLONG)での太陽定数 DOUBLE PRECISION XI0 !"XI の春分点での値 * (-1) DOUBLE PRECISION TIME0 !"春分点での時間 * (-1) DOUBLE PRECISION LAT(-NLAT:NLAT) !"緯度 DOUBLE PRECISION TIME(0:NTIME) & !"時間(春分点から測る, YTIMEで規格化) DOUBLE PRECISION PHIHARU(0:NTIME) !"春分点から測った黄経 C" DOUBLE PRECISION F(0:NTIME, -NLAT:NLAT) ! 日平均日射量 REAL F(0:NTIME, -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) 'SOLARC ' WRITE(OR,630) SOLARC WRITE(OR,650) 'XI0 ' WRITE(OR,630) XI0 WRITE(OR,650) 'TIME0 ' WRITE(OR,630) TIME0 WRITE(OR,650) 'LAT ' WRITE(OR,620) LAT WRITE(OR,650) 'TIME ' WRITE(OR,620) TIME WRITE(OR,650) 'PHIHARU ' WRITE(OR,620) PHIHARU WRITE(OR,650) 'F ' WRITE(OR,620) F *----------------------------------------------------------------------- RETURN END