表層部分の土壌への浸入と蒸発散量の計算のためのプログラム


C***** ﾊﾟﾗﾒ-ﾀの値
        SMMAX		：最大土壌水分量
        PET		：蒸発散比
	COFHR		：Horton式の逓減係数 COef Horton Reduction
	SM18		：地下への浸透を始める土壌水分量 SoilMoisture pf1.8
	INFFS		：最終浸入能＝地下浸透量 INFiltration Final
	COFEV		：蒸発散係数 COef EVapotranspiration 

      	INFF18		：地下への浸透を始める浸入能　INFiltration pf1.8
    	INFMAX		：土壌が最も乾燥したときの浸入能
                          SMMAXから逆算する

        SW18=SMMAX*SM18

      CALL CALFINFF18(INFFS,COFHR,SW18,INFF18)	    ：土壌水分量SW18からINFINFF18を求める

      CALL CALFMAX(INFF18,INFFS,COFHR,SMMAX,INFMAX) ：土壌水分量SMMAXから浸入能を計算する

      OLA=PET/0.01	：蒸発散計算のための定数
      ES=EXP(-COFHR)
      ES1=1.0-ES
C***** 時間（DO 1）
      DO 1 J=1,24
C*** 現在の土壌水分(SW)の計算
        SW=0.0
        CALL SWCAL(INF(J),SW,INFF18,FS,COFHR,INFMAX)
C*** 蒸発散比(ETR)の算出と蒸発散量(ActE)の決定
        EB=EXP(-COFEV*SW)
        ETR=PET/(1.0+OLA*EB)
        ActE=PE*ETR
C*** 浸入能と地表流（SurfQ）の計算
        RAIN=RAIN-ActE
        IF(RAIN.GT.0.0) THEN
C                                      ** 浸入能の減衰と地表流の計算
             CALL INFILS(RAIN,J,ActE)
C                                      ** 地表中(有効雨量) = 降雨量 − 浸み込んだ量
C            ただし，単位時間降雨量が単位時間に浸透しうる最大量よりも大きくない
C            と地表流は発生しない．・・・すべて土壌中にしみ込んでしまう．
             SurfQ=RAIN-(INFFS+(INF(J)-INFFS)*ES1/COFHR)
             IF(SurfQ.LT.0.0) SurfQ=0.0

        ELSE
C                                      ** 蒸発による浸入能の回復
            IF(INF(J).LE.INFMAX) THEN
C                                      ** 浸入能回復可能
                 CALL INFILS(RAIN,J,ActE)
                 SurfQ=0.0
             ELSE
C                                      ** 乾燥しきっているため浸入能回復不可能
                 INF(J+1)=INFMAX
                 ActE=0.0
                 SurfQ=0.0
            ENDIF
        ENDIF

C*** 地下への浸透計算(InfQ：地下への浸透量)

        IF(INF(J+1).GE.INFF18) THEN
            InfQ=0.0
        ELSE
            InfQ=INFFS
        ENDIF

    1 CONTINUE
      END

C ***** 地表面における浸入能計算ｻﾌﾞﾙ-ﾁﾝ
      SUCOFHRROUTINE INFILS(RR,K,ActE)
C ******************************************
      RI=RR
           IF (INF(K).LT.INFF18) RI=RI-INFFS
      IF(RI.LT.0.0) THEN
C                       ***** RI(<0)に相当する土壌水分だけ浸入能回復
         CALL INFIL2(RI,K)
         IF(INF(K+1).GT.INFMAX) THEN
C                       ***** ただし土壌が乾燥しきっていたら回復でき
C                             ないので，修正する
            ActE=(INFMAX-INF(K))/COFHR-LOG((INF(K)-INFFS)/(INFMAX-
     +            INFFS))*INFFS/COFHR
            INF(K+1)=INFMAX
         ENDIF
      ELSE
C                       *** RI(>0)に相当する土壌水分だけ浸入能減衰

C                       *** FSS 単位時間に土壌中にしみ込みうる土壌水分量の最大値
         FSS=INFFS+(INF(K)-INFFS)*ES1/COFHR
         IF(RI.GE.FSS) THEN
C                       ***** 単位降雨量が FSS を超えるとき
               IF(INF(K).GT.INFF18) THEN
C                       ***** Horton式そのままに浸入能は減衰する
                      INF(K+1)=INF(K)*ES+INFFS*ES1
                ELSE
C                       ***** 地下浸透が始まっていたら減衰量が変化するので
C                             地下浸透量を考慮した場合の減衰量を再度求める
                      R0=(INF(K)-INFFS)*ES1/COFHR
                      CALL INFIL2(R0,K)
                ENDIF
         ELSE
                CALL INFIL2(RI,K)
         ENDIF

       ENDIF
      RETURN
      END
C ***** 浸入能（積分の数値計算）ｻﾌﾞﾙ-ﾁﾝ
      SUCOFHRROUTINE INFIL2(RI,K)
C *************************************
C 解析的に次の時刻の浸入能を求められない場合は，数値計算によって求める
      IF(RI.EQ.0) GO TO 2
      DT1=1
      DO 1 I=1,10
      FSUM=FS*DT1-1.0/COFHR*(1-EXP(-COFHR*DT1))*(FS-INF(K))
    1 DT1=RI/FSUM*DT1
      INF(K+1)=FS+(INF(K)-FS)*EXP(-COFHR*DT1)
      GO TO 3
    2 INF(K+1)=INF(K)
    3 RETURN
      END
C ***** 浸入能最大値の計算ｻﾌﾞﾙ-ﾁﾝ
      SUCOFHRROUTINE CALFMAX(INFF18,INFFS,COFHR,SMMAX,INFMAX)
C ********************************************
        A=INFF18-INFFS
        DD=1.0
        X=INFF18
    1   FINF=-INFFS*LOG(A/X)+X
        XR=SMMAX*COFHR-FINF
        IF(ABS(XR).LE.0.01) GO TO 100
        IF(XR.GT.0.0) THEN
            X=X+DD*INFF18
          ELSE
            X=X-DD*INFF18
            DD=DD/10.0
            X=X+DD*INFF18
        ENDIF
        GO TO 1
  100 INFMAX=X+INFFS
      RETURN
      END
C ****** 土壌水分量の計算サブルーチン
      SUCOFHRROUTINE SWCAL(SJ,SW,INFF18,INFFS,COFHR,INFMAX)
C **********************************************
      INF=SJ
      SL=(INF-INFFS)/(INFMAX-INFFS)
      SF=(INFF18-INFFS)/(INFMAX-INFFS)
      IF(INF.GE.INFF18)THEN
        SW=(INFMAX-INF)/COFHR-INFFS*LOG(SL)/COFHR
      ELSE
        SW=(INFMAX-INF)/COFHR-INFFS*LOG(SF)/COFHR
      ENDIF
      SWMAX=(INFMAX-INFFS)/COFHR-LOG((INFF18-INFFS)/(INFMAX-INFFS))*INFFS/COFHR
      SW=SW/SWMAX*100.0
      IF(SW.LT.0.0)SW=0.0
      RETURN
      END

C ***** 地下浸透を開始する浸入能の計算ｻﾌﾞﾙ-ﾁﾝ
      SUCOFHRROUTINE CALFINFF18(INFFS,COFHR,SM18,INFF18)
C ********************************************
      DT1=1
      INF=INFFS*1.1
      DO 1 I=1,10
      FSUM=INFFS*DT1-1.0/COFHR*(1-EXP(-COFHR*DT1))*(INFFS-INF)
    1 DT1=-SW18/FSUM*DT1
      INFF18=INFFS+(INF-INFFS)*EXP(-COFHR*DT1)
      RETURN
      END