新安江模型.bas

可直接用于新安江模型计算

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