fit_tankmodel.f90 Source File

program main
    use TankModelClass
    implicit none

    type(IO_) :: fw,fp
    type(TankModel_) :: model
    real(real64),allocatable :: rainfall(:),waterheads(:,:),params(:)
    ! time increment
    real(real64) :: dt
    ! tank model fitting
    ! unit: mm, s
    ! 表面排水=0.1 mm/s/mm, 1mmH2Oで0.1mm/s排水
    ! 表層30cmとして，-20cmに表面排水，-40cmに耕盤，その下に耕盤下タンク
    ! rain_profile.txt is generated by generate_rain_profile.f90

    rainfall = to_vector(fp%to_array("rain_profile.txt",column=[2],header=0))
    
    call model%init()
    call model%add(to_tank( &
            coeff=[0.0205d0/60.0d0/60.0d0,0.030d0/60.0d0/60.0d0],&
            height=[200.0d0,0.0d0], &
            GL_bottom = -300.0d0, &
            name = "Surface Layer" &
            ) &
        )
    call model%add(to_tank(&
            coeff=[0.0010d0/60.0d0/60.0d0], &
            height=[0.0d0], &
            GL_bottom=-600.0d0 , &
            name = "Bottom Layer" &
            ) &
        )
    call model%connect(tankIdx=[1,2],coeff=0.0090d0/60.0d0/60.0d0,height=0.0d0)
    
    dt = 60.0d0

    waterheads = this%simulate(rainfall=rainfall,GLs=[-5.0d0,-50.0d0,-100.0d0,-200.0d0],dt=dt)

    call fw%open("waterhead_true.txt","w")
    call fw%write(linspace([0.0d0,dt*size(rainfall)],size(rainfall))/60.0d0/60.0d0/24.0d0 .h. waterheads)
    call fw%close()
    ! 正解値はここまで．    
    
    ! 試行は以下


    call model%init()
    call model%add(to_tank( &
            coeff=[1.0d0/60.0d0/60.0d0,1.0d0/60.0d0/60.0d0],&
            height=[200.0d0,0.0d0], &
            GL_bottom = -300.0d0 , &
            name = "Surface Layer" &
            ) &
        )
    call model%add(to_tank(&
            coeff=[1.0d0/60.0d0/60.0d0], &
            height=[0.0d0], &
            GL_bottom=-600.0d0 , &
            name = "Bottom Layer" &
            ) &
        )
    call model%connect(tankIdx=[1,2],coeff=1.0d0/60.0d0/60.0d0,height=0.0d0)

    waterheads = fp%to_array("waterhead_true.txt",column=[2,3,4,5],header=0)
    
    params = model%optimize(rainfall=rainfall,waterheads=waterheads,GLs=[-5.0d0,-50.0d0,-100.0d0,-200.0d0],dt=dt)
    
    call model%setParams(params)

    waterheads = model%simulate(rainfall=rainfall,GLs=[-5.0d0,-50.0d0,-100.0d0,-200.0d0],dt=dt)
    
    call fw%open("waterhead_sim.txt","w")
    call fw%write(linspace([0.0d0,dt*size(rainfall)],size(rainfall))/60.0d0/60.0d0/24.0d0 .h. waterheads)
    call fw%close()
    
    print *, "True ans:"
    print *, [0.0205d0/60.0d0/60.0d0,0.030d0/60.0d0/60.0d0] // [0.0090d0/60.0d0/60.0d0] // [0.0010d0/60.0d0/60.0d0] 

end program main