SW_VAR_DENS_MODEL.f90

Go to the documentation of this file.

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

PROGRAM sw_var_dens_model

  USE constitutive_2d, ONLY : init_problem_param

  USE geometry_2d, ONLY : init_grid
  USE geometry_2d, ONLY : init_source
  USE geometry_2d, ONLY : topography_reconstruction

  USE geometry_2d, ONLY : dx,dy,b_cent
  USE geometry_2d, ONLY : comp_cells_x,comp_cells_y

  USE init_2d, ONLY : riemann_problem
  USE init_2d, ONLY : collapsing_volume

  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_cell
  USE solver_2d, ONLY : timestep
  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 : topo_change_flag
  USE parameters_2d, ONLY : verbose_level
  USE parameters_2d, ONLY : n_solid
  USE parameters_2d, ONLY : n_vars
  USE parameters_2d, ONLY : radial_source_flag
  USE parameters_2d, ONLY : collapsing_volume_flag

  USE solver_2d, ONLY : q , qp , t, dt
  USE solver_2d, ONLY : q1max
  
  USE constitutive_2d, ONLY : qc_to_qp

  USE solver_2d, ONLY : solve_mask

  IMPLICIT NONE

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

  INTEGER j,k

  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 ( collapsing_volume_flag ) THEN

        CALL collapsing_volume

     END IF

     IF( riemann_flag ) THEN

        ! riemann problem defined in file.inp
        CALL riemann_problem

     ENDIF

  END IF

  IF ( radial_source_flag ) CALL init_source

  CALL check_solve
  
  IF ( topo_change_flag ) CALL topography_reconstruction

  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.

  DO k = 1,comp_cells_y
     
     DO j = 1,comp_cells_x
        
        CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars,j,k) )
        
     END DO
     
  END DO

  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
  
  WRITE(*,fmt="(A3,F10.4,A5,F9.5,A9,ES11.3E3,A11,ES11.3E3,A9,ES11.3E3,A15,ES11.3E3)")   &
       't =',t,'dt =',dt0,                                                      &
       ' mass = ',dx*dy*sum(q(1,:,:)) ,                                         &
       ' volume = ',dx*dy*sum(qp(1,:,:)) ,                                      &
       ' area = ',dx*dy*count(q(1,:,:).GT.1.d-5) ,                              &
       ' solid mass = ',dx*dy*sum(q(5:4+n_solid,:,:)) 

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

     CALL update_param
     
     CALL check_solve

     IF ( verbose_level .GE. 1 ) THEN

        WRITE(*,*) 'cells to solve and reconstruct:' , count(solve_mask)

     END IF

     CALL timestep

     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

     CALL update_erosion_deposition_cell(dt)
 
     IF ( topo_change_flag ) CALL topography_reconstruction

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

     DO k = 1,comp_cells_y
        
        DO j = 1,comp_cells_x
           
           CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars,j,k) )
           
        END DO

     END DO
  

     WRITE(*,fmt="(A3,F10.4,A5,F9.5,A9,ES11.3E3,A11,ES11.3E3,A9,ES11.3E3,A15,ES11.3E3)")&
          't =',t,'dt =',dt,                                                    &
          ' mass = ',dx*dy*sum(q(1,:,:)) ,                                      &
          ' volume = ',dx*dy*sum(qp(1,:,:)) ,                                   &
          ' area = ',dx*dy*count(q(1,:,:).GT.1.d-7) ,                           &
          ' solid mass = ',dx*dy*sum(q(5:4+n_solid,:,:)) 
     
     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

     t_steady = t_end

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

        CALL output_solution(t)

        CALL cpu_time(t2)

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

     END IF

  END DO

  CALL deallocate_solver_variables

  CALL close_units

  CALL cpu_time(t2)

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

END PROGRAM sw_var_dens_model