I am using Fortran to resolve a series of sediment transport equations to solve changes in bed elevation, water depth, velocity, and sediment transport in both time and space for an infinitely large prismatic channel. Below is my nested loop structure. I have no errors coming up but the variables update in time but not space. Also, if I don't run the initialization as a separate do loop then nothing populates at all in the output text file. I'm not sure why this isn't working. The variables from the input file for Q, z0, slope, mann, domain, and y0 are as follows

13.00 -11.00 0.00002 0.0143 400000.00 11

program my_project
    implicit none

    ! Variable Declarations
    real :: Q, z0, y0, slope, mann, domain, S_b, S_f, Fr, dx, dt, g, U, C, p, ps, ds, d, dz_dx
    integer :: i, n_steps, n_time_steps, t_index, x_index
    real, dimension(:), allocatable :: X, z_vals, y_vals, U_vals, q_t, dz_vals

    ! Read Input File
    open(11, file='proj1pt1data.txt')
    read(11, *) Q, z0, slope, mann, domain, y0
    close(11)

    ! Constants
    g = 9.81
    dx = 100.0
    dt = 1.0  
    n_steps = int(domain / dx)  
    n_time_steps = 100  

    ! Morphodynamic Constants
    p = 0.4
    ps = 2650
    ds = 0.0001
    d = (ps - 1000) / 1000

    ! Allocate Arrays
    allocate(X(n_steps), z_vals(n_steps), y_vals(n_steps), U_vals(n_steps), q_t(n_steps), dz_vals(n_steps))

    ! Open Output File
    open(12, file='mydata.dat')
    write(12,*) 'Time (s)', 'X (m)', 'z (m)', 'Depth (m)', 'U (m/s)', 'qt (m2/s)'

    ! ***********************
    ! Initialization (Runs Once)
    ! ***********************
    do i = 1, n_steps
        X(i) = (i - 1) * dx
        z_vals(i) = z0 + slope * X(i)  
        y_vals(i) = y0  
        U_vals(i) = Q / y0   
    end do




    ! ***********************
    ! Start Time Loop 
    ! ***********************
    do t_index = 1, n_time_steps  
        print *, "t=", t_index, "Before update: y(1)=", y_vals(1), "y(2)=", y_vals(2)
        ! Start Spatial Loop
        do x_index = 2, n_steps-1

            ! Compute Froude Number
            Fr = U_vals(x_index) / sqrt(g * y_vals(x_index))

            ! Compute energy dissipation (S_f)
            S_f = mann**2 * U_vals(x_index)**2 / y_vals(x_index)**(4.0/3.0)
            
            ! Set known bed slope as S_b 
            S_b = slope

            ! Compute Change in Depth 
            y_vals(x_index) = y_vals(x_index) + dx * ((S_b - S_f) / (1.0 - Fr**2))
       
            ! Apply physical constraints
            if (y_vals(x_index) < 0.0) y_vals(x_index) = 0.0

            ! Compute velocity based on new depth
            U_vals(x_index) = Q / y_vals(x_index)
          
            ! Compute Sediment Transport Rate
            C = 1.0 / mann * sqrt(g / y_vals(x_index))
            q_t(x_index) = 0.05 * (U_vals(x_index)**5) / (C**3 * sqrt(g) * d**2 * ds)
            
            ! Compute Bed Elevation Change (Now Coupled)
            dz_dx = - (q_t(x_index+1) - q_t(x_index)) / ((1.0 - p) * 2 * dx)
           
            z_vals(x_index) = z_vals(x_index) + dz_dx * dt
            
! Write to output file
            write(12,*) t_index * dt, X(x_index), z_vals(x_index), y_vals(x_index), U_vals(x_index), q_t(x_index)
      
        end do ! End Spatial Loop


        

    end do ! End Time Loop  

    print *, "t=", t_index, "Before update: y(1)=", y_vals(1), "y(2)=", y_vals(2)

      
    ! Close Output File
    close(12)

    ! Deallocate Arrays
    deallocate(X, z_vals, y_vals, U_vals, q_t, dz_vals)

end program my_project

I have tried re-arranging this many ways and can get it to update in either time or space but never both and receive no error messages. The output should be 100 time points and 4000 spatial points. Every variable should update in both time and space. For example, along the directional axis (x) bed elevation (z), depth (y), velocity (u) and sediment transport (q_t) should differ. Based on these changes along the x-axis, each of these variables will change within their x designation at each time step until the system reaches equilibrium. I can get z to change in along the directional and time axis but y, u, and q_t only change as a function of time with each value being the same along x. Below are the equations being solved

steady state 1-d

dy/dx=(S_b-S_f)/(1-〖Fr〗^2 )

Euler's method:

 y^(n+1)=y^n+dx((S_b-S_f)/(1-〖Fr〗^2 ))_n

Engulend-Hansen Equation:

q^T=0.05 u^5/(C^3 〖√g ∆〗^2 d_s )

Exener equation:

〖(1-p)((dz_b)/dt)〗i=-(q_i^T (t)-q(i_UPWIND)^T (t))/∆x

I am using Fortran to resolve a series of sediment transport equations to solve changes in bed elevation, water depth, velocity, and sediment transport in both time and space for an infinitely large prismatic channel. Below is my nested loop structure. I have no errors coming up but the variables update in time but not space. Also, if I don't run the initialization as a separate do loop then nothing populates at all in the output text file. I'm not sure why this isn't working. The variables from the input file for Q, z0, slope, mann, domain, and y0 are as follows

13.00 -11.00 0.00002 0.0143 400000.00 11

program my_project
    implicit none

    ! Variable Declarations
    real :: Q, z0, y0, slope, mann, domain, S_b, S_f, Fr, dx, dt, g, U, C, p, ps, ds, d, dz_dx
    integer :: i, n_steps, n_time_steps, t_index, x_index
    real, dimension(:), allocatable :: X, z_vals, y_vals, U_vals, q_t, dz_vals

    ! Read Input File
    open(11, file='proj1pt1data.txt')
    read(11, *) Q, z0, slope, mann, domain, y0
    close(11)

    ! Constants
    g = 9.81
    dx = 100.0
    dt = 1.0  
    n_steps = int(domain / dx)  
    n_time_steps = 100  

    ! Morphodynamic Constants
    p = 0.4
    ps = 2650
    ds = 0.0001
    d = (ps - 1000) / 1000

    ! Allocate Arrays
    allocate(X(n_steps), z_vals(n_steps), y_vals(n_steps), U_vals(n_steps), q_t(n_steps), dz_vals(n_steps))

    ! Open Output File
    open(12, file='mydata.dat')
    write(12,*) 'Time (s)', 'X (m)', 'z (m)', 'Depth (m)', 'U (m/s)', 'qt (m2/s)'

    ! ***********************
    ! Initialization (Runs Once)
    ! ***********************
    do i = 1, n_steps
        X(i) = (i - 1) * dx
        z_vals(i) = z0 + slope * X(i)  
        y_vals(i) = y0  
        U_vals(i) = Q / y0   
    end do




    ! ***********************
    ! Start Time Loop 
    ! ***********************
    do t_index = 1, n_time_steps  
        print *, "t=", t_index, "Before update: y(1)=", y_vals(1), "y(2)=", y_vals(2)
        ! Start Spatial Loop
        do x_index = 2, n_steps-1

            ! Compute Froude Number
            Fr = U_vals(x_index) / sqrt(g * y_vals(x_index))

            ! Compute energy dissipation (S_f)
            S_f = mann**2 * U_vals(x_index)**2 / y_vals(x_index)**(4.0/3.0)
            
            ! Set known bed slope as S_b 
            S_b = slope

            ! Compute Change in Depth 
            y_vals(x_index) = y_vals(x_index) + dx * ((S_b - S_f) / (1.0 - Fr**2))
       
            ! Apply physical constraints
            if (y_vals(x_index) < 0.0) y_vals(x_index) = 0.0

            ! Compute velocity based on new depth
            U_vals(x_index) = Q / y_vals(x_index)
          
            ! Compute Sediment Transport Rate
            C = 1.0 / mann * sqrt(g / y_vals(x_index))
            q_t(x_index) = 0.05 * (U_vals(x_index)**5) / (C**3 * sqrt(g) * d**2 * ds)
            
            ! Compute Bed Elevation Change (Now Coupled)
            dz_dx = - (q_t(x_index+1) - q_t(x_index)) / ((1.0 - p) * 2 * dx)
           
            z_vals(x_index) = z_vals(x_index) + dz_dx * dt
            
! Write to output file
            write(12,*) t_index * dt, X(x_index), z_vals(x_index), y_vals(x_index), U_vals(x_index), q_t(x_index)
      
        end do ! End Spatial Loop


        

    end do ! End Time Loop  

    print *, "t=", t_index, "Before update: y(1)=", y_vals(1), "y(2)=", y_vals(2)

      
    ! Close Output File
    close(12)

    ! Deallocate Arrays
    deallocate(X, z_vals, y_vals, U_vals, q_t, dz_vals)

end program my_project

I have tried re-arranging this many ways and can get it to update in either time or space but never both and receive no error messages. The output should be 100 time points and 4000 spatial points. Every variable should update in both time and space. For example, along the directional axis (x) bed elevation (z), depth (y), velocity (u) and sediment transport (q_t) should differ. Based on these changes along the x-axis, each of these variables will change within their x designation at each time step until the system reaches equilibrium. I can get z to change in along the directional and time axis but y, u, and q_t only change as a function of time with each value being the same along x. Below are the equations being solved

steady state 1-d

dy/dx=(S_b-S_f)/(1-〖Fr〗^2 )

Euler's method:

 y^(n+1)=y^n+dx((S_b-S_f)/(1-〖Fr〗^2 ))_n

Engulend-Hansen Equation:

q^T=0.05 u^5/(C^3 〖√g ∆〗^2 d_s )

Exener equation:

〖(1-p)((dz_b)/dt)〗i=-(q_i^T (t)-q(i_UPWIND)^T (t))/∆x

• fortran

Share Improve this question Follow edited Mar 23 at 16:07 user30026192 asked Mar 22 at 20:21 user30026192user30026192 1133 bronze badges 24

• Welcome, I suggest to read the tour and How to Ask. – Vladimir F Героям слава Commented Mar 22 at 20:39
• You have written your question for somebody knowledgable in your field. There will be few such people here, but many experienced programmers - I am a computational chemist and really can't tell much from your description. Tell us exactly what you expect, with values, what you get, with values, and why that is wrong. Best provide a minimal, reproducible example – Ian Bush Commented Mar 22 at 21:05
• "ridiculous comments with asterisks. Makes it too hard to read" - I beg to differ @steve . – Ian Bush Commented Mar 22 at 21:10
• Thank you both for your input, I've updated the code. I agree with Ian that the comments make it easier to follow along but removed them before reading his comment. I hope this updated question will help. I'm sorry if the question is not concise as expected, I am beginning to learn this language and am struggling to apply it even though the math is making sense to me. Thank you all for your help! – user30026192 Commented Mar 22 at 21:16
• @steve thank you for pointing this issue out and I apologize for the oversight on my part. I think everything should be there for you now, please let me know if there is anything else I can add. – user30026192 Commented Mar 22 at 22:29

 |  Show 19 more comments

1 Answer 1

Sorted by: Reset to default 3

Well, you can have a play with the code below. Note that it is set for 10 years, with a timestep of 1 day. However, with your current parameters (particularly slope) the sediment transport is negligible.

Equations to be solved are:

(1) depth (y) - from the GVF (gradually-varied-flow) equation. Note that, for your subcritical conditions, the flow is controlled from DOWNSTREAM (i.e. sea level if it is a tidal estuary), so the boundary condition is at xmax:

(2) bed level (z) - from conservation of sediment:

(3) sediment total load (I think you have coded this wrongly for the Engelund-Hansen model):

Sample Code

program my_project
    implicit none

    real :: Q, z0, y0, slope, mann, domain, S_b, S_f, Fr, dx, dt, g, p, ps, ds, d
    integer :: i, n_steps, n_time_steps, t_index, x_index
    real dydx, dzdt
    real, dimension(:), allocatable :: X, z_vals, y_vals, U_vals, q_t, dz_vals

!   open(11, file='proj1pt1data.txt')
!   read(11, *) Q, z0, slope, mann, domain, y0
!   close(11)
    Q = 13.00;   z0 = -11.00;   slope = 0.00002;   mann =  0.0143;   domain = 400000.00;   y0 = 11

    ! Constants
    g = 9.81                                     ! gravity (m/s^2)
    dx = 100.0                                   ! spatial step (m)
    dt = 60 * 60 * 24                            ! timestep (s) = one day
    n_steps = int(domain / dx)                   ! number of spatial steps
    n_time_steps = 365 * 10                      ! number of timesteps

    ! Morphodynamic Constants
    p = 0.4                                      ! porosity
    ps = 2650                                    ! sediment density (kg/m^3)
    ds = 0.0001                                  ! sediment diameter (m)
    d = (ps - 1000) / 1000                       ! relative density increase (relative to water)

    allocate(X(n_steps), z_vals(n_steps), y_vals(n_steps), U_vals(n_steps), q_t(n_steps), dz_vals(n_steps))

    open(12, file='mydata.dat')
    write(12,*) 'X (m)', 'z (m)', 'Depth (m)', 'U (m/s)', 'qt (m2/s)'

    !*****************
    ! Initialization *
    !*****************
    do i = 1, n_steps
        X(i) = (i - 1) * dx                      ! spatial location
        z_vals(i) = z0 - slope * X(i)            ! bed level (goes DOWNHILL!)
        y_vals(i) = y0                           ! OUTLET depth
        U_vals(i) = Q / y0                       ! velocity
    end do

    !***********************************************************************
    ! Solve the GVF (gradually-varied-flow) equation for depth and velocity
    ! Runs BEFORE the timestepping loop
    !***********************************************************************
    !*************************************************
    ! Timestepping loop - bed-load sediment transport
    !*************************************************
    do t_index = 1, n_time_steps  

       ! Calculate depth from GVF equation. NOTE: SUBCRITICAL FLOW IS CONTROLLED FROM DOWNSTREAM
       do x_index = n_steps - 1, 1, -1                                         ! downstream depth is fixed
          Fr = U_vals(x_index) / sqrt( g * y_vals(x_index) )                   ! Froude number
          S_f = mann ** 2 * U_vals(x_index)**2 / y_vals(x_index) ** (4.0/3.0)  ! friction slope
!         S_b = slope                                                          ! (downward) bed slope (fixed)
          S_b = -( z_vals(x_index+1) - z_vals(x_index) ) / dx                  ! (downward) bed slope (with sediment transport)
          dydx = ( S_b - S_f ) / ( 1.0 - Fr ** 2 )                             ! depth gradient (GVF equation)
          y_vals(x_index) = y_vals(x_index+1) - dx * dydx                      ! Euler update
       end do

       ! Set velocity and bed load
       do x_index = 1, n_steps
          U_vals(x_index) = Q / y_vals(x_index)                                ! velocity (from continuity)
          q_t(x_index) = 0.05 * U_vals(x_index) ** 5 * mann ** 3 / ( sqrt( g * y_vals(x_index) ) * d ** 2 * ds )
                                                                               ! sediment load (Engelund-Hansen)
       end do

       ! Compute bed elevation change (from (1-p) dz/dt = -dq/dx )
       do x_index = 2, n_steps-1
!         dzdt = -( q_t(x_index+1) - q_t(x_index-1) ) / ( 2 * dx ) / ( 1.0 - p )
          dzdt = -( q_t(x_index) - q_t(x_index-1) ) / dx / ( 1.0 - p )
          z_vals(x_index) = z_vals(x_index) + dzdt * dt
       end do

       ! Boundary conditions (extrapolate from interior)
       z_vals(1) = z_vals(2) + z_vals(2) - z_vals(3)
       z_vals(n_steps) = z_vals(n_steps-1) + z_vals(n_steps-1) - z_vals(n_steps-2)

    end do ! End Spatial Loop

    ! Final output
    do x_index = 1, n_steps
        write(12,*) X(x_index), z_vals(x_index), y_vals(x_index), U_vals(x_index), q_t(x_index)
    end do
    close(12)

end program my_project

Depth against distance (note that the boundary condition is DOWNSTREAM)

Bed sediment level (barely changed and nowhere near equilibrium):