! ! ------------------------------------------------------------- ! Sample program of parameter calibration (2-series tank model) ! ------------------------------------------------------------- ! !---------------------------------------------------------------------- ! Module for Common Data !---------------------------------------------------------------------- module common_data implicit none integer, parameter, public :: ndmax = 100000 integer, public :: nd ! number of data real, public :: p(ndmax) ! precipitation (mm/day) real, public :: e(ndmax) ! evapotranspiration (mm/day) real, public :: q(ndmax) ! discharge (mm/day) end module common_data !---------------------------------------------------------------------- ! Main Program !---------------------------------------------------------------------- program test4 use common_data implicit none integer, parameter :: nmax = 100 integer :: n, i real :: xinit(1:nmax), xmin(1:nmax), xmax(1:nmax) real :: x(1:nmax), rmse character(len=3) :: xname(1:nmax) real :: qc(1:ndmax) data n / 8 / data xname(1),xmin(1),xmax(1),xinit(1) / 'a11', 0., 1., 0.1 / data xname(2),xmin(2),xmax(2),xinit(2) / 'a12', 0., 1., 0.01 / data xname(3),xmin(3),xmax(3),xinit(3) / 'a2 ', 0., 1., 0.01 / data xname(4),xmin(4),xmax(4),xinit(4) / 'b ', 0., 1., 0.01 / data xname(5),xmin(5),xmax(5),xinit(5) / 'z11', 0., 400., 50. / data xname(6),xmin(6),xmax(6),xinit(6) / 'z12', 0., 200., 10. / data xname(7),xmin(7),xmax(7),xinit(7) / 'h1 ', 0., 500., 0. / data xname(8),xmin(8),xmax(8),xinit(8) / 'h2 ', 0., 5000., 300. / !====== Read input date ====== open(1, file='indata.txt', status='old') do i=1, ndmax read(1, *, end=99) p(i), e(i), q(i) end do 99 continue close(1) nd = i - 1 !====== Optimize ====== call optimize(n, xmin, xmax, xinit, x, rmse) !====== Calculate discharge with optimum parameters ====== call run_tank2(x, nd, p, e, qc) !====== Print result ====== do i=1, n print *, xname(i), ' =', x(i) end do print * print *, 'RMSE =', rmse print * !====== Write output date ====== open(2, file='outdata.txt') do i=1, nd write(2, *) p(i), e(i), q(i), qc(i) end do close(2) end program test4 !---------------------------------------------------------------------- ! Run 2-seriese Tank Model !---------------------------------------------------------------------- subroutine run_tank2(x, nd, p, e, q) implicit none real, intent(in) :: x(1:8) ! parameters integer, intent(in) :: nd ! number of data real, intent(in) :: p(1:nd) ! rainfall (mm/day) real, intent(in) :: e(1:nd) ! evapotranspiration (mm/day) real, intent(out) :: q(1:nd) ! simulated discharge (mm/day) real :: a11, a12, a2, b, z11, z12, h1, h2 real :: dt, q11, q12, q2, r integer :: nt, i, j a11 = x(1); a12 = x(2); a2 = x(3); b = x(4) z11 = x(5); z12 = x(6); h1 = x(7); h2 = x(8) nt = 6 dt = 1. / real(nt) do i=1, nd q(i) = 0.0 do j=1, nt q11 = a11 * max(h1 - z11, 0.0) ! surface flow q12 = a12 * max(h1 - z12, 0.0) ! sub-surface flow r = b * h1 ! infiltration q2 = a2 * h2 ! base flow h1 = h1 + (p(i) - e(i) - q11 - q12 - r) * dt h2 = h2 + (r - q2) * dt h2 = h2 + min(h1, 0.0) ! evaporate from h2 if h1 is empty h1 = max(h1, 0.0) h2 = max(h2, 0.0) q(i) = q(i) + q11 + q12 + q2 end do end do end subroutine run_tank2 !---------------------------------------------------------------------- ! Evaluation Function !---------------------------------------------------------------------- function ev_tank2(x, ev) result(r) use common_data implicit none real, intent(in) :: x(:) ! parameter list real, intent(out) :: ev ! evaluation value integer :: r ! number of constrained condition real :: qc(1:nd) real :: rmse !====== Check constraint conditions ====== r = 0 if(x(1) < x(2)) r = r + 1 ! a11 < a12 if(x(5) < x(6)) r = r + 1 ! z11 < z12 if(x(1) + x(2) + x(4) > 1.0) r = r + 1 ! a11+a12+b > 1 !====== Run model and evaluate ====== if(r == 0) then call run_tank2(x, nd, p, e, qc) ev = rmse(nd, q, qc) else ev = 0.0 end if end function ev_tank2 !---------------------------------------------------------------------- ! Optimization !---------------------------------------------------------------------- subroutine optimize(n, xmin, xmax, xinit, x, rmse) use common_data use opti_sceua use opti_pso implicit none integer, intent(in) :: n real, intent(in) :: xmin(1:n), xmax(1:n), xinit(1:n) real, intent(out) :: x(1:n), rmse character(len=80) :: pfname integer :: algorithm external :: ev_tank2 namelist /ctrl_opti/ algorithm !====== Read control parameter file ====== pfname = 'param_tank2.txt' open(3, file=pfname, status='old') read(3, ctrl_opti) close(3) !====== Optimize ====== select case (algorithm) case(1) call do_sceuaf(n, xmin, xmax, xinit, sce_minimize, ev_tank2, & & pfname, x, rmse) case(2) call do_psof(n, xmin, xmax, xinit, pso_minimize, ev_tank2, & & pfname, x, rmse) case default print *, 'invalid algorithm' stop end select end subroutine optimize !---------------------------------------------------------------------- ! Root Mean Square Error !---------------------------------------------------------------------- function rmse(n, x, y) result(r) implicit none integer, intent(in) :: n real, intent(in) :: x(1:n), y(1:n) real :: r real :: se(1:n) integer :: i forall(i=1: n) se(i) = (x(i) - y(i))**2 r = sqrt(sum(se(1:n) / real(n))) end function rmse