IMEX_Sflow_2d.f90

Go to the documentation of this file.

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

PROGRAM imex_sflow_2d

  USE constitutive_2d, ONLY : init_problem_param

  USE geometry_2d, ONLY : init_grid

  USE geometry_2d, ONLY : dx,dy,b_cent,cell_size

  USE init_2d, ONLY : riemann_problem, initial_conditions

  USE inpout_2d, ONLY : init_param
  USE inpout_2d, ONLY : read_param
  USE inpout_2d, ONLY : update_param
  USE inpout_2d, ONLY : output_solution
  USE inpout_2d, ONLY : output_runout
  USE inpout_2d, ONLY : read_solution
  USE inpout_2d, ONLY : close_units

  USE inpout_2d, ONLY : output_runout_flag


  USE solver_2d, ONLY : allocate_solver_variables
  USE solver_2d, ONLY : deallocate_solver_variables
  USE solver_2d, ONLY : imex_rk_solver
  USE solver_2d, ONLY : update_erosion_deposition_ver
  USE solver_2d, ONLY : timestep, timestep2
  ! USE solver_2d, ONLY : check_solve

  USE inpout_2d, ONLY : restart

  USE parameters_2d, ONLY : t_start
  USE parameters_2d, ONLY : t_end
  USE parameters_2d, ONLY : t_output
  USE parameters_2d, ONLY : t_runout
  USE parameters_2d, ONLY : t_steady
  USE parameters_2d, ONLY : dt0
  USE parameters_2d, ONLY : riemann_flag
  USE parameters_2d, ONLY : verbose_level
  USE parameters_2d, ONLY : solid_transport_flag


  USE solver_2d, ONLY : q , dt
  USE solver_2d, ONLY : q1max
  
  ! USE solver_2d, ONLY : solve_mask

  IMPLICIT NONE

  REAL*8 :: t
  REAL*8 :: t1 , t2
  REAL*8 :: dt_old , dt_old_old
  LOGICAL :: stop_flag
  LOGICAL :: stop_flag_old

  CALL cpu_time(t1)

  CALL init_param

  CALL read_param

  CALL init_grid

  CALL init_problem_param

  CALL allocate_solver_variables
  
  IF ( restart ) THEN

     CALL read_solution

  ELSE

     IF( riemann_flag .EQV. .true. )THEN

        ! riemann problem defined in file.inp
        CALL riemann_problem

     ELSE

        ! generic problem defined by initial conditions function (in init_2d.f90)
        CALL initial_conditions

     ENDIF

  END IF

  t = t_start
  
  WRITE(*,*) 
  WRITE(*,*) '******** START COMPUTATION *********'
  WRITE(*,*)

  IF ( verbose_level .GE. 1 ) THEN
     
     WRITE(*,*) 'Min q(1,:,:)=',minval(q(1,:,:))
     WRITE(*,*) 'Max q(1,:,:)=',maxval(q(1,:,:))

     WRITE(*,*) 'Min B(:,:)=',minval(b_cent(:,:))
     WRITE(*,*) 'Max B(:,:)=',maxval(b_cent(:,:))


     WRITE(*,*) 'size B_cent',size(b_cent,1),size(b_cent,2)

     WRITE(*,*) 'SUM(q(1,:,:)=',sum(q(1,:,:))
     WRITE(*,*) 'SUM(B_cent(:,:)=',sum(b_cent(:,:))
     
  END IF

  q1max(:,:) = q(1,:,:)
  
  dt_old = dt0
  dt_old_old = dt_old
  t_steady = t_end
  stop_flag = .false.


  IF ( output_runout_flag ) CALL output_runout(t,stop_flag)

  CALL output_solution(t)

  IF ( sum(q(1,:,:)) .EQ. 0.d0 ) t_steady = t_output
  
  IF ( solid_transport_flag ) THEN

     WRITE(*,fmt="(A3,F10.4,A5,F9.5,A15,ES12.3E3,A15,ES12.3E3,A15,ES12.3E3)")   &
          't =',t,'dt =',dt0,                                                   &
          ' total volume = ',dx*dy*sum(q(1,:,:)) ,                              &
          ' total area = ',dx*dy*count(q(1,:,:).GT.1.d-5) ,                     &
          ' total sediment fraction = ',dx*dy*sum(q(4,:,:)) 
          
  ELSE

     WRITE(*,fmt="(A3,F10.4,A5,F9.5,A15,ES12.3E3,A15,ES12.3E3)")                &
          't =',t,'dt =',dt0,                                                   &
          ' total volume = ',dx*dy*sum(q(1,:,:))  ,                             &
          ' total area = ',dx*dy*count(q(1,:,:).GT.1.d-5)

  END IF

  DO WHILE ( ( t .LT. t_end ) .AND. ( t .LT. t_steady ) )

     CALL update_param
     
     ! CALL check_solve
     ! WRITE(*,*) 'cells to solve:',COUNT(solve_mask)


     ! CALL timestep

     CALL timestep2

     IF ( t+dt .GT. t_end ) dt = t_end - t
     IF ( t+dt .GT. t_output ) dt = t_output - t
     
     IF ( output_runout_flag ) THEN

        IF ( t+dt .GT. t_runout ) dt = t_runout - t

     END IF

     dt = min(dt,1.1d0*0.5d0*(dt_old+dt_old_old))
     
     dt_old_old = dt_old
     dt_old = dt

     CALL imex_rk_solver

     IF ( solid_transport_flag ) THEN
                
        CALL update_erosion_deposition_ver(dt)
 
     END IF

     q1max(:,:) = max( q1max(:,:) , q(1,:,:) )
     
     t = t+dt

     IF ( solid_transport_flag ) THEN
        
        WRITE(*,fmt="(A3,F10.4,A5,F9.5,A15,ES12.3E3,A15,ES12.3E3,A15,ES12.3E3)")&
             't =',t,'dt =',dt,                                                 &
             ' total volume = ',dx*dy*sum(q(1,:,:)) ,                           &
             ' total area = ',dx*dy*count(q(1,:,:).GT.1.d-7) ,                  &
             ' total sediment fraction = ',dx*dy*sum(q(4,:,:)) 
        
     ELSE
        
        WRITE(*,fmt="(A3,F10.4,A5,F9.5,A15,ES12.3E3,A15,ES12.3E3)")             &
             't =',t,'dt =',dt,                                                 &
             ' total volume = ',dx*dy*sum(q(1,:,:))  ,                          &
             ' total area = ',dx*dy*count(q(1,:,:).GT.1.d-7)
        
     END IF
     
     IF ( output_runout_flag ) THEN

        IF ( ( t .GE. t_runout ) .OR. ( t .GE. t_steady ) ) THEN

           stop_flag_old = stop_flag
           
           IF ( output_runout_flag ) CALL output_runout(t,stop_flag)

           IF ( ( stop_flag ) .AND. (.NOT.stop_flag_old) ) THEN
              
              t_steady = min(t_end,t_output)
              t_runout = t_steady

           END IF

        END IF

     END IF

     IF ( ( t .GE. t_output ) .OR. ( t .GE. t_end ) ) THEN

        CALL output_solution(t)

     END IF

  END DO

  CALL deallocate_solver_variables

  CALL close_units

  CALL cpu_time(t2)

  WRITE(*,*) 'Time taken by the code was',t2-t1,'seconds'

END PROGRAM imex_sflow_2d