module bmproblem

  implicit none

  ! PARAMETERS
  integer, parameter :: nobs = 122, nqbound = 94, nzb = 39, nwidth = 39, nx = 100, nsv = nx !, nsv = 2 * nx
  real(kind=8), parameter :: qboundx = - 39.0d0

  ! LOCAL SCALARS
  integer :: i
  
  ! LOCAL ARRAYS
  real(kind=8) :: qboundb(nqbound),qboundc(nqbound),qboundd(nqbound),zbb(nzb),zbc(nzb),zbd(nzb), &
       widthb(nwidth),widthc(nwidth),widthd(nwidth),svx(nsv),svt(nsv)

  real(kind=8) :: qboundt(nqbound) = 3600.0d0 * (/ &
       0.00000E+00, 0.12000E+02, 0.24000E+02, 0.36000E+02, &
       0.48000E+02, 0.60000E+02, 0.72000E+02, 0.84000E+02, &
       0.96000E+02, 0.10800E+03, 0.12000E+03, 0.13200E+03, &
       0.14400E+03, 0.15600E+03, 0.16800E+03, 0.18000E+03, &
       0.19200E+03, 0.20400E+03, 0.21600E+03, 0.22800E+03, &
       0.24000E+03, 0.25200E+03, 0.26400E+03, 0.27600E+03, &
       0.28800E+03, 0.30000E+03, 0.31200E+03, 0.32400E+03, &
       0.33600E+03, 0.34800E+03, 0.36000E+03, 0.37200E+03, &
       0.38400E+03, 0.39600E+03, 0.40800E+03, 0.42000E+03, &
       0.43200E+03, 0.44400E+03, 0.45600E+03, 0.46800E+03, &
       0.48000E+03, 0.49200E+03, 0.50400E+03, 0.51600E+03, &
       0.52800E+03, 0.54000E+03, 0.55200E+03, 0.56400E+03, &
       0.57600E+03, 0.58800E+03, 0.60000E+03, 0.61200E+03, &
       0.62400E+03, 0.63600E+03, 0.64800E+03, 0.66000E+03, &
       0.67200E+03, 0.68400E+03, 0.69600E+03, 0.70800E+03, &
       0.72000E+03, 0.73200E+03, 0.74400E+03, 0.75600E+03, &
       0.76800E+03, 0.78000E+03, 0.79200E+03, 0.80400E+03, &
       0.81600E+03, 0.82800E+03, 0.84000E+03, 0.85200E+03, &
       0.86400E+03, 0.87600E+03, 0.88800E+03, 0.90000E+03, &
       0.91200E+03, 0.92400E+03, 0.93600E+03, 0.94800E+03, &
       0.96000E+03, 0.97200E+03, 0.98400E+03, 0.99600E+03, &
       0.10080E+04, 0.10200E+04, 0.10320E+04, 0.10440E+04, &
       0.10560E+04, 0.10680E+04, 0.10800E+04, 0.10920E+04, &
       0.11040E+04, 0.11160E+04 /)

  real(kind=8) :: qboundval(nqbound) = (/ &
       0.39000E+01, 0.39000E+01, 0.49800E+01, 0.57400E+01, &
       0.72000E+01, 0.67800E+01, 0.87600E+01, 0.10700E+02, &
       0.12400E+02, 0.13000E+02, 0.14500E+02, 0.17500E+02, &
       0.16700E+02, 0.20900E+02, 0.18400E+02, 0.25200E+02, &
       0.20900E+02, 0.24700E+02, 0.22300E+02, 0.25500E+02, &
       0.24900E+02, 0.33100E+02, 0.28200E+02, 0.30000E+02, &
       0.25500E+02, 0.29100E+02, 0.27400E+02, 0.16100E+02, &
       0.12700E+02, 0.94000E+01, 0.76200E+01, 0.53500E+01, &
       0.49800E+01, 0.41800E+01, 0.41800E+01, 0.55800E+01, &
       0.66200E+01, 0.96800E+01, 0.10500E+02, 0.14600E+02, &
       0.13400E+02, 0.19000E+02, 0.15100E+02, 0.13600E+02, &
       0.10600E+02, 0.76200E+01, 0.59800E+01, 0.41100E+01, &
       0.37600E+01, 0.32300E+01, 0.36900E+01, 0.67000E+01, &
       0.65400E+01, 0.94000E+01, 0.87600E+01, 0.15200E+02, &
       0.11800E+02, 0.19000E+02, 0.12800E+02, 0.15400E+02, &
       0.98500E+01, 0.94000E+01, 0.71200E+01, 0.74500E+01, &
       0.52000E+01, 0.62900E+01, 0.48300E+01, 0.36900E+01, &
       0.30400E+01, 0.25400E+01, 0.20200E+01, 0.19700E+01, &
       0.16500E+01, 0.26700E+01, 0.20200E+01, 0.27300E+01, &
       0.20800E+01, 0.27900E+01, 0.20800E+01, 0.29800E+01, &
       0.22500E+01, 0.26700E+01, 0.21900E+01, 0.32300E+01, &
       0.21900E+01, 0.23700E+01, 0.19100E+01, 0.24300E+01, &
       0.18600E+01, 0.23700E+01, 0.19100E+01, 0.17000E+01, &
       0.16500E+01, 0.13900E+01 /)
  
  real(kind=8) :: zbx(nzb) = (/ &
       0.00000E+00, 0.88000E+02, 0.14800E+03, 0.20900E+03, &
       0.29500E+03, 0.38200E+03, 0.47800E+03, 0.56600E+03, &
       0.64800E+03, 0.74600E+03, 0.83400E+03, 0.90000E+03, &
       0.97800E+03, 0.10620E+04, 0.11740E+04, 0.12600E+04, &
       0.13550E+04, 0.14260E+04, 0.14900E+04, 0.15610E+04, &
       0.15940E+04, 0.17750E+04, 0.18600E+04, 0.19410E+04, &
       0.20150E+04, 0.21010E+04, 0.21790E+04, 0.22710E+04, &
       0.23580E+04, 0.24480E+04, 0.25480E+04, 0.26540E+04, &
       0.27400E+04, 0.28350E+04, 0.29550E+04, 0.30360E+04, &
       0.31190E+04, 0.31810E+04, 0.32130E+04 /)

  real(kind=8) :: zbval(nzb) = (/ &
       0.74600E+01, 0.69100E+01, 0.71300E+01, 0.72100E+01, &
       0.71400E+01, 0.74200E+01, 0.70400E+01, 0.68100E+01, &
       0.70100E+01, 0.72100E+01, 0.65600E+01, 0.70200E+01, &
       0.68300E+01, 0.66100E+01, 0.68000E+01, 0.64100E+01, &
       0.65600E+01, 0.66200E+01, 0.65000E+01, 0.64800E+01, &
       0.62500E+01, 0.62100E+01, 0.62000E+01, 0.61400E+01, &
       0.62500E+01, 0.63800E+01, 0.59900E+01, 0.62900E+01, &
       0.58500E+01, 0.58700E+01, 0.55700E+01, 0.54800E+01, &
       0.53900E+01, 0.56400E+01, 0.52700E+01, 0.53800E+01, &
       0.51500E+01, 0.49400E+01, 0.49400E+01 /)

  real(kind=8) :: widthx(nwidth) = (/ &
       0.00000E+00, 0.88000E+02, 0.14800E+03, 0.20900E+03, &
       0.29500E+03, 0.38200E+03, 0.47800E+03, 0.56600E+03, &
       0.64800E+03, 0.74600E+03, 0.83400E+03, 0.90000E+03, &
       0.97800E+03, 0.10620E+04, 0.11740E+04, 0.12600E+04, &
       0.13550E+04, 0.14260E+04, 0.14900E+04, 0.15610E+04, &
       0.15940E+04, 0.17750E+04, 0.18600E+04, 0.19410E+04, &
       0.20150E+04, 0.21010E+04, 0.21790E+04, 0.22710E+04, &
       0.23580E+04, 0.24480E+04, 0.25480E+04, 0.26540E+04, &
       0.27400E+04, 0.28350E+04, 0.29550E+04, 0.30360E+04, &
       0.31190E+04, 0.31810E+04, 0.32130E+04 /)

  real(kind=8) :: widthval(nwidth) = (/ &
       0.16500E+02, 0.11500E+02, 0.14500E+02, 0.17500E+02, &
       0.15500E+02, 0.25500E+02, 0.21000E+02, 0.15500E+02, &
       0.20500E+02, 0.19000E+02, 0.20500E+02, 0.20000E+02, &
       0.24000E+02, 0.15500E+02, 0.20000E+02, 0.15000E+02, &
       0.14500E+02, 0.22500E+02, 0.19000E+02, 0.17000E+02, &
       0.16500E+02, 0.26000E+02, 0.12000E+02, 0.22500E+02, &
       0.18000E+02, 0.22000E+02, 0.22000E+02, 0.30500E+02, &
       0.14000E+02, 0.18500E+02, 0.16000E+02, 0.19500E+02, &
       0.14500E+02, 0.23000E+02, 0.14000E+02, 0.16000E+02, &
       0.14500E+02, 0.15000E+02, 0.16000E+02 /)

  logical :: obs(nobs) = .true.

  real(kind=8), parameter :: obsx1 = 751.0d0, obsx2 = 3256.0d0
  
  real(kind=8) :: obsx(nobs) = (/ ( obsx1, i=1,60 ), ( obsx2, i=61,122 ) /)

  real(kind=8) :: obst(nobs) = 3600.0d0 * (/ &
       0.72000E+02, 0.84000E+02, 0.96000E+02, 0.10800E+03, &
       0.12000E+03, 0.13200E+03, 0.14400E+03, 0.15600E+03, &
       0.16800E+03, 0.18000E+03, 0.19200E+03, 0.20400E+03, &
       0.21600E+03, 0.22800E+03, 0.24000E+03, 0.25200E+03, &
       0.26400E+03, 0.27600E+03, 0.28800E+03, 0.30000E+03, &
       0.31200E+03, 0.32400E+03, 0.33600E+03, 0.34800E+03, &
       0.36000E+03, 0.37200E+03, 0.38400E+03, 0.39600E+03, &
       0.40800E+03, 0.42000E+03, 0.43200E+03, 0.44400E+03, &
       0.45600E+03, 0.46800E+03, 0.48000E+03, 0.49200E+03, &
       0.50400E+03, 0.51600E+03, 0.52800E+03, 0.54000E+03, &
       0.55200E+03, 0.56400E+03, 0.57600E+03, 0.58800E+03, &
       0.60000E+03, 0.61200E+03, 0.62400E+03, 0.63600E+03, &
       0.64800E+03, 0.66000E+03, 0.67200E+03, 0.68400E+03, &
       0.69600E+03, 0.70800E+03, 0.72000E+03, 0.73200E+03, & 
       0.74400E+03, 0.75600E+03, 0.76800E+03, 0.78000E+03, &
       0.72000E+02, 0.84000E+02, 0.96000E+02, 0.10800E+03, &
       0.12000E+03, 0.13200E+03, 0.14400E+03, 0.15600E+03, &
       0.16800E+03, 0.18000E+03, 0.19200E+03, 0.20400E+03, &
       0.21600E+03, 0.22800E+03, 0.24000E+03, 0.25200E+03, &
       0.26400E+03, 0.27600E+03, 0.28800E+03, 0.30000E+03, &
       0.31200E+03, 0.32400E+03, 0.33600E+03, 0.34800E+03, &
       0.36000E+03, 0.37200E+03, 0.38400E+03, 0.39600E+03, &
       0.40800E+03, 0.42000E+03, 0.43200E+03, 0.44400E+03, &
       0.45600E+03, 0.46800E+03, 0.48000E+03, 0.49200E+03, &
       0.50400E+03, 0.51600E+03, 0.52800E+03, 0.54000E+03, &
       0.55200E+03, 0.56400E+03, 0.57600E+03, 0.58800E+03, &
       0.60000E+03, 0.61200E+03, 0.62400E+03, 0.63600E+03, &
       0.64800E+03, 0.66000E+03, 0.67200E+03, 0.68400E+03, &
       0.69600E+03, 0.70800E+03, 0.72000E+03, 0.73200E+03, & 
       0.74400E+03, 0.75600E+03, 0.76800E+03, 0.78000E+03, &
       0.79200E+03, 0.80400E+03 /)
  
  real(kind=8) :: obsz(nobs) = (/ &
       0.77600E+01, 0.78600E+01, 0.79300E+01, 0.79600E+01, &
       0.80300E+01, 0.81500E+01, 0.81000E+01, 0.82100E+01, &
       0.81400E+01, 0.82600E+01, 0.81900E+01, 0.82500E+01, &
       0.82000E+01, 0.82700E+01, 0.82500E+01, 0.83800E+01, &
       0.83100E+01, 0.83300E+01, 0.82600E+01, 0.83000E+01, &
       0.82900E+01, 0.80900E+01, 0.79100E+01, 0.77100E+01, &
       0.76100E+01, 0.74600E+01, 0.74300E+01, 0.73800E+01, &
       0.73900E+01, 0.75200E+01, 0.75300E+01, 0.77700E+01, &
       0.77700E+01, 0.80000E+01, 0.79300E+01, 0.81200E+01, &
       0.80100E+01, 0.79200E+01, 0.77800E+01, 0.76200E+01, &
       0.75200E+01, 0.73700E+01, 0.73600E+01, 0.73600E+01, &
       0.73500E+01, 0.75800E+01, 0.75500E+01, 0.77400E+01, &
       0.76800E+01, 0.80000E+01, 0.78400E+01, 0.81100E+01, &
       0.79100E+01, 0.80300E+01, 0.77500E+01, 0.77300E+01, &
       0.75900E+01, 0.76300E+01, 0.75100E+01, 0.75800E+01, &
       0.58200E+01, 0.59700E+01, 0.59900E+01, 0.60900E+01, &
       0.60900E+01, 0.62800E+01, 0.61800E+01, 0.63600E+01, &
       0.62300E+01, 0.64300E+01, 0.63000E+01, 0.63900E+01, &
       0.63100E+01, 0.64100E+01, 0.63700E+01, 0.66100E+01, &
       0.64800E+01, 0.65300E+01, 0.63900E+01, 0.64600E+01, &
       0.64500E+01, 0.61700E+01, 0.60000E+01, 0.58200E+01, &
       0.57300E+01, 0.56300E+01, 0.55800E+01, 0.55500E+01, &
       0.55200E+01, 0.56400E+01, 0.56500E+01, 0.58400E+01, &
       0.58400E+01, 0.60700E+01, 0.59700E+01, 0.61800E+01, &
       0.60800E+01, 0.60200E+01, 0.58700E+01, 0.57200E+01, &
       0.56400E+01, 0.55200E+01, 0.55100E+01, 0.55100E+01, &
       0.55000E+01, 0.56900E+01, 0.56700E+01, 0.58300E+01, &
       0.57700E+01, 0.60600E+01, 0.59100E+01, 0.61800E+01, &
       0.60000E+01, 0.61100E+01, 0.58500E+01, 0.58300E+01, &
       0.56900E+01, 0.57400E+01, 0.56000E+01, 0.56900E+01, &
       0.56000E+01, 0.55300E+01 /)

contains

  ! ******************************************************************
  ! ******************************************************************

  subroutine iniprob()

    implicit none

    ! PARAMETERS
    integer, parameter :: npoints = 10000

    ! LOCAL SCALARS
    integer :: i
    real(kind=8) :: qb,t,w,x,zb
    
    call nspline(nqbound,qboundt,qboundval,qboundb,qboundc,qboundd)
  
    call nspline(nzb,zbx,zbval,zbb,zbc,zbd)

    call nspline(nwidth,widthx,widthval,widthb,widthc,widthd)

    ! Output files to generate plots
    
    open(unit=10,file='qboundgiven.txt')
    do i = 1,nqbound
       write(10,*) qboundt(i),qboundval(i)
    end do
    close(10)
    
    open(unit=10,file='qboundspline.txt')
    do i = 1,npoints
       t = qboundt(1) + ( qboundt(nqbound) - qboundt(1) ) * ( i - 1 ) / ( npoints - 1 )
       call qbound(t,qb)
       write(10,*) t,qb
    end do
    close(10)
    
    open(unit=10,file='wgiven.txt')
    do i = 1,nwidth
       write(10,*) widthx(i),widthval(i)
    end do
    close(10)
    
    open(unit=10,file='wspline.txt')
    do i = 1,npoints
       x = min( qboundx, widthx(1) ) + &
            ( max( maxval( obsx(1:nobs) ), widthx(nwidth) ) - min( qboundx, widthx(1) ) ) * ( i - 1 ) / ( npoints - 1 )
       call width(x,w)
       write(10,*) x,w
    end do
    close(10)
    
    open(unit=10,file='zbgiven.txt')
    do i = 1,nzb
       write(10,*) zbx(i),zbval(i)
    end do
    close(10)
    
    open(unit=10,file='zbspline.txt')
    do i = 1,npoints
       x = min( qboundx, zbx(1) ) + &
            ( max( maxval( obsx(1:nobs) ), zbx(nzb) ) - min( qboundx, zbx(1) ) ) * ( i - 1 ) / ( npoints - 1 )
       call bedel(x,zb)
       write(10,*) x,zb
    end do
    close(10)

    open(unit=10,file='observations-x751.txt')
    open(unit=20,file='observations-x3256.txt')

    do i = 1,nobs
       if ( obs(i) ) then
          if ( obsx(i) .eq. obsx1 ) then
             write(10,*) obsx(i),obst(i),obsz(i)
          else ! if ( obsx(i) .eq. obsx2 ) then
             write(20,*) obsx(i),obst(i),obsz(i)
          end if
       end if
    end do
    
    close(10)
    close(20)

  end subroutine iniprob
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine setnsv(DayE)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: DayE
    
    ! LOCAL SCALARS
    integer :: j
    real(kind=8) :: tmax,x,xend,xini
    
    xini = min( qboundx, minval( obsx(1:nobs), obs(1:nobs) ) )
    xend = maxval( obsx(1:nobs), obs(1:nobs) )

    tmax = maxval( obst(1:nobs), obs(1:nobs) )

    do j = 1,nx
       x = xini + ( xend - xini ) * ( j - 1 ) / ( 100 - 1 )

       svx(j) = x
       svt(j) = DayE * 24.0d0 * 3600.0d0
    end do

  end subroutine setnsv
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine endprob()

    implicit none

  end subroutine endprob
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine qbound(t,qb)

    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: t
    real(kind=8), intent(out) :: qb

    ! FUNCTIONS
    real(kind=8), external :: seval
    
    qb = seval(nqbound,t,qboundt,qboundval,qboundb,qboundc,qboundd)
  
  end subroutine qbound
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine dqbound(t,dqb)

    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: t
    real(kind=8), intent(out) :: dqb
    
    ! FUNCTIONS
    real(kind=8), external :: dseval
    
    dqb = dseval(nqbound,t,qboundt,qboundval,qboundb,qboundc,qboundd)
  
  end subroutine dqbound
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine ddqbound(t,ddqb)

    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: t
    real(kind=8), intent(out) :: ddqb
    
    ! FUNCTIONS
    real(kind=8), external :: ddseval
    
    ddqb = ddseval(nqbound,t,qboundt,qboundval,qboundb,qboundc,qboundd)
  
  end subroutine ddqbound
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine width(x,w)

    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: w
    
    ! FUNCTIONS
    real(kind=8), external :: seval
    
    w = seval(nwidth,x,widthx,widthval,widthb,widthc,widthd)
  
  end subroutine width
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine dwidth(x,dw)

    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: dw

    ! FUNCTIONS
    real(kind=8), external :: dseval
    
    dw = dseval(nwidth,x,widthx,widthval,widthb,widthc,widthd)
  
  end subroutine dwidth
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine bedel(x,zb)
    
    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: zb
    
    ! FUNCTIONS
    real(kind=8), external :: seval
    
    zb = seval(nzb,x,zbx,zbval,zbb,zbc,zbd)
  
  end subroutine bedel

  ! ******************************************************************
  ! ******************************************************************

  subroutine dbedel(x,dzb)

    implicit none
    
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: dzb
    
    ! FUNCTIONS
    real(kind=8), external :: dseval
    
    dzb = dseval(nzb,x,zbx,zbval,zbb,zbc,zbd)
  
  end subroutine dbedel

end module bmproblem
