Attribute VB_Name = "Module1"
Public n As Integer, m As Integer, par(20) As Single, area As Single, uh() As Single
Public dt As Single, p() As Single, ep() As Single, qr() As Single, w(3) As Single
Public fr As Single, s As Single, qrss0 As Single, qrg0 As Single
'參數說明：apr1 上層張力水容量wum ，par2下層張力水容量 wlm，par3深層張力水容量wdm，par4 蒸發能力折算系數
'par5 深層蒸發系數c，par6 張力水蓄水容量系數b，par7不透水面積比率imp1，par8 自由水蓄水容量sm，par9 自由水蓄水容量指數ex
'par10地下水出流系數kg，par11壤中流出流系數kss，par12地下水出流系數kkg，par13 壤中流出流系數kkss，area 單元面積，uh 無因次單位線
'dt 時段步長 ，p降雨系列，ep 蒸發皿蒸發能力，qr 單元出流，w土壤含水量1上層，2下層，3深層
'fr 初始產流面積，s初始自由水深，qrss0初始壤中流流量，qrg0 初始地下水徑流量
Public Sub xajmx(n, m, par, area, uh, dt, p, ep, qr, w, fr, s, qrss0, qrg0)
      Dim d As Integer
      Dim kssd As Single, kgd As Single, kg As Single, kss As Single, kkss As Single, kkg As Single
      Dim kc As Single, imp1 As Single, e(3) As Single, wm(3) As Single
      For i = 1 To 3
         wm(i) = par(i)
      Next i
      kc = par(4)
      c = par(5)
      b = par(6)
      imp1 = par(7)
      sm = par(8)
      ex = par(9)
      kg = par(10)
      kss = par(11)
      kkg = par(12)
      kkss = par(13)
      For i = 1 To n
      qr(i) = 0
      Next i
      u = area / 3.6 / dt
      If dt <= 24 Then
        d = 24 / dt
        ci = kkss ^ (1 / d)
        cg = kkg ^ (1 / d)
        kssd = (1 - (1 - (kg + kss)) ^ (1 / d)) / (1 + kg / kss)
        kgd = kssd * kg / kss
      Else
      MsgBox "計算時段長度不合適"
      End If
      For i = 1 To n
          If ep(i) < 0 Then ep(i) = 0
          If p(i) < 0 Then p(i) = 0
          ep(i) = ep(i) * kc
          wm0 = wm(1) + wm(2) + wm(3)
          w0 = w(1) + w(2) + w(3)
          pe = p(i) - ep(i)
           r = 0
          rimp = 0
          If pe > 0 Then
             wmm = (1 + b) * wm0 / (1 - imp1)
             If wm0 - w0 <= 0.0001 Then
                a = wmm
             Else
             a = wmm * (1 - (1 - w0 / wm0) ^ (1 / (1 + b)))
             End If
             If pe + a < wmm Then
                r = pe - wm0 + w0 + wm0 * ((1 - (pe + a) / wmm) ^ (1 + b))
             Else
                r = pe - (wm0 - w0)
             End If
             rimp = pe * imp1
          End If
            ' c____________________________________________________________
            If w(1) + p(i) > ep(i) Then
                e(1) = ep(i)
                e(2) = 0
                e(3) = 0
            Else
                e(1) = w(1) + p(i)
                e(2) = (ep(i) - e(1)) * w(2) / wm(2)
                If w(2) <= c * wm(2) Then
                    e(2) = c * (ep(i) - e(1))
                    e(3) = 0
                If w(2) >= c * (ep(i) - e(1)) Then
                    e(2) = c * (ep(i) - e(1))
                    e(3) = 0
                Else
                    e(2) = w(2)
                    e(3) = c * (ep(i) - e(1) - e(2))
                End If
                End If
            End If
            
           w(1) = w(1) + p(i) - r - e(1)
           w(2) = w(2) - e(2)
           w(3) = w(3) - e(3)
           If w(1) > wm(1) Then
              w(2) = w(1) - wm(1) + w(2)
              w(1) = wm(1)
              If w(2) > wm(2) Then
                 w(3) = w(3) + w(2) - wm(2)
                 w(2) = wm(2)
              End If
            End If
            X = fr
            If pe <= 0 Then
               rs = 0
               rss = s * kssd * fr
               rgd = s * fr * kgd
               s = s - (rss + rg) / fr
            Else
            fr = r / pe
            s = X * s / fr
            ss = s
            q = r / fr
            nn = Int(q / 5#) + 1
            q = q / nn
            kssdd = (1 - (1 - (kgd + kssd)) ^ (1 / nn)) / (1 + kgd / kssd)
            kgdd = kssdd * kgd / kssd
            rs = 0
            rss = 0
            rg = 0
            smm = (1 + ex) * sm
            If ex < 0.001 Then
               smmf = smm
            Else
               smmf = smm * (1 - (1 - fr) ^ (1 / ex))
            End If
            smf = smmf / (1 + ex)
            For j = 1 To nn
                If s > smf Then s = smf
                au = smmf * (1 - (1 - s / smf) ^ (1 / (1 + ex)))
                If q + au <= 0 Then
                   rsd = 0
                   rssd = 0
                   rgd = 0
                   s = 0
                ElseIf q + au >= smmf Then
                    rsd = (q + s - smf) * fr
                    rssd = smf * kssdd * fr
                    rgd = smf * fr * kgdd
                    s = smf - (kssd + kgd) / fr
                    ElseIf q + au < smmf Then
                         rsd = (q - smf + s + smf * (1 - (q + au) / smmf) ^ (1 + ex)) * fr
                         rssd = (s + q - rsd / fr) * kssdd * fr
                         rgd = (s + q - rsd / fr) * kgdd * fr
                         s = s + q - (rsd + rssd + rgd) / fr
                    End If
                    rs = rs + rsd
                    rss = rss + rssd
                    rg = rg + rgd
            Next j
        End If
        rs = rs * (1 - imp1)
        rss = rss * (1 - imp1)
        rg = rg * (1 - imp1)
        qrs = (rs + rimp) * u
        qrss = qrss0 * ci + rss * (1 - ci) * u
        qrg = qrg0 * cg + rg * (1 - cg) * u
        qtr = qrs + qrss + qrg
        For j = 1 To m
            If i + j - 1 <= n Then
               qr(i + j - 1) = qr(i + j - 1) + qtr * uh(j)
            End If
        Next j
        qrss0 = qrss
        qrg0 = qrg
        Next i
End Sub