solver_2d.f90

Go to the documentation of this file.

!*******************************************************************************
!
!
!
!********************************************************************************
MODULE solver_2d

  ! external variables

  USE constitutive_2d, ONLY : implicit_flag

  USE geometry_2d, ONLY : comp_cells_x,comp_cells_y
  USE geometry_2d, ONLY : comp_interfaces_x,comp_interfaces_y

  USE geometry_2d, ONLY : b_cent , b_prime_x , b_prime_y , b_stag_x , b_stag_y
  USE geometry_2d, ONLY : grav_surf

  USE parameters_2d, ONLY : n_eqns , n_vars , n_nh
  USE parameters_2d, ONLY : n_rk
  USE parameters_2d, ONLY : verbose_level

  USE parameters_2d, ONLY : bcw , bce , bcs , bcn

  IMPLICIT none

  REAL*8, ALLOCATABLE :: q(:,:,:)        
  REAL*8, ALLOCATABLE :: q0(:,:,:)        
  REAL*8, ALLOCATABLE :: q_fv(:,:,:)     
  REAL*8, ALLOCATABLE :: q_interfacel(:,:,:)        
  REAL*8, ALLOCATABLE :: q_interfacer(:,:,:)
  REAL*8, ALLOCATABLE :: q_interfaceb(:,:,:)        
  REAL*8, ALLOCATABLE :: q_interfacet(:,:,:)
  REAL*8, ALLOCATABLE :: a_interface_xneg(:,:,:)
  REAL*8, ALLOCATABLE :: a_interface_xpos(:,:,:)
  REAL*8, ALLOCATABLE :: a_interface_yneg(:,:,:)
  REAL*8, ALLOCATABLE :: a_interface_ypos(:,:,:)
  REAL*8, ALLOCATABLE :: h_interface_x(:,:,:)
  REAL*8, ALLOCATABLE :: h_interface_y(:,:,:)
  REAL*8, ALLOCATABLE :: qp(:,:,:)

  REAL*8, ALLOCATABLE :: source_xy(:,:)

  LOGICAL, ALLOCATABLE :: solve_mask(:,:) , solve_mask0(:,:)

  REAL*8 :: dt

  LOGICAL, ALLOCATABLE :: mask22(:,:) , mask21(:,:) , mask11(:,:) , mask12(:,:)

  INTEGER :: i_rk

  REAL*8, ALLOCATABLE :: a_tilde_ij(:,:)
  REAL*8, ALLOCATABLE :: a_dirk_ij(:,:)

  REAL*8, ALLOCATABLE :: omega_tilde(:)
  REAL*8, ALLOCATABLE :: omega(:)

  REAL*8, ALLOCATABLE :: a_tilde(:)
  REAL*8, ALLOCATABLE :: a_dirk(:)
  REAL*8 :: a_diag

  REAL*8, ALLOCATABLE :: q_rk(:,:,:,:)
  REAL*8, ALLOCATABLE :: divflux(:,:,:,:)
  REAL*8, ALLOCATABLE :: nh(:,:,:,:)
  REAL*8, ALLOCATABLE :: si_nh(:,:,:,:)

  REAL*8, ALLOCATABLE :: expl_terms(:,:,:,:)

  REAL*8, ALLOCATABLE :: divfluxj(:,:)
  REAL*8, ALLOCATABLE :: nhj(:,:)
  REAL*8, ALLOCATABLE :: expl_terms_j(:,:)
  REAL*8, ALLOCATABLE :: si_nhj(:,:)

  LOGICAL :: normalize_q

  LOGICAL :: normalize_f

  LOGICAL :: opt_search_nl

  REAL*8, ALLOCATABLE :: residual_term(:,:,:)


  REAL*8 :: t_imex1,t_imex2

CONTAINS

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

  SUBROUTINE allocate_solver_variables

    IMPLICIT NONE

    REAL*8 :: gamma , delta

    INTEGER :: i,j

    ALLOCATE( q( n_vars , comp_cells_x , comp_cells_y ) , q0( n_vars ,          &
         comp_cells_x , comp_cells_y ) )

    ALLOCATE( q_fv( n_vars , comp_cells_x , comp_cells_y ) )

    ALLOCATE( q_interfacel( n_vars , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( q_interfacer( n_vars , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( a_interface_xneg( n_eqns , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( a_interface_xpos( n_eqns , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( h_interface_x( n_eqns , comp_interfaces_x, comp_cells_y ) )


    ALLOCATE( q_interfaceb( n_vars , comp_cells_x, comp_interfaces_y ) )
    ALLOCATE( q_interfacet( n_vars , comp_cells_x, comp_interfaces_y ) )
    ALLOCATE( a_interface_yneg( n_eqns , comp_cells_x, comp_interfaces_y ) )
    ALLOCATE( a_interface_ypos( n_eqns , comp_cells_x, comp_interfaces_y ) )


    ALLOCATE( h_interface_y( n_eqns , comp_cells_x, comp_interfaces_y ) )

    ALLOCATE( solve_mask( comp_cells_x , comp_cells_y ) )
    ALLOCATE( solve_mask0( comp_cells_x , comp_cells_y ) )

    ALLOCATE( qp( n_vars , comp_cells_x , comp_cells_y ) )

    ALLOCATE( source_xy( comp_cells_x , comp_cells_y ) )

    
    ALLOCATE( a_tilde_ij(n_rk,n_rk) )
    ALLOCATE( a_dirk_ij(n_rk,n_rk) )
    ALLOCATE( omega_tilde(n_rk) )
    ALLOCATE( omega(n_rk) )


    ! Allocate the logical arrays defining the implicit parts of the system
    ALLOCATE( mask22(n_eqns,n_eqns) )
    ALLOCATE( mask21(n_eqns,n_eqns) )
    ALLOCATE( mask11(n_eqns,n_eqns) )
    ALLOCATE( mask12(n_eqns,n_eqns) )

    ! Initialize the logical arrays with all false (everything is implicit)
    mask11(1:n_eqns,1:n_eqns) = .false.
    mask12(1:n_eqns,1:n_eqns) = .false.
    mask22(1:n_eqns,1:n_eqns) = .false.
    mask21(1:n_eqns,1:n_eqns) = .false.

    ! Set to .TRUE. the elements not corresponding to equations and variables to 
    ! be solved implicitly
    DO i = 1,n_eqns

       DO j = 1,n_eqns

          IF ( .NOT.implicit_flag(i) .AND. .NOT.implicit_flag(j) )              &
               mask11(j,i) = .true.
          IF ( implicit_flag(i) .AND. .NOT.implicit_flag(j) )                   &
               mask12(j,i) = .true.
          IF ( implicit_flag(i) .AND. implicit_flag(j) )                        &
               mask22(j,i) = .true.
          IF ( .NOT.implicit_flag(i) .AND. implicit_flag(j) )                   &
               mask21(j,i) = .true.

       END DO

    END DO

    ! Initialize the coefficients for the IMEX Runge-Kutta scheme
    ! Please note that with respect to the schemes described in Pareschi & Russo 
    ! (2000) we do not have the coefficient vectors c_tilde and c, because the 
    ! explicit and implicit terms do not depend explicitly on time.

    ! Explicit part coefficients (a_tilde_ij=0 for j>=i)
    a_tilde_ij = 0.d0

    ! Weight coefficients of the explicit part in the final assemblage
    omega_tilde = 0.d0

    ! Implicit part coefficients (a_dirk_ij=0 for j>i)
    a_dirk_ij = 0.d0

    ! Weight coefficients of the explicit part in the final assemblage
    omega = 0.d0

    gamma = 1.d0 - 1.d0 / sqrt(2.d0)
    delta = 1.d0 - 1.d0 / ( 2.d0 * gamma )

    IF ( n_rk .EQ. 1 ) THEN

       a_tilde_ij(1,1) = 1.d0

       omega_tilde(1) = 1.d0

       a_dirk_ij(1,1) = 0.d0

       omega(1) = 0.d0

    ELSEIF ( n_rk .EQ. 2 ) THEN

       a_tilde_ij(2,1) = 1.0d0

       omega_tilde(1) = 1.0d0
       omega_tilde(2) = 0.0d0

       a_dirk_ij(2,2) = 1.0d0

       omega(1) = 0.d0
       omega(2) = 1.d0

      
    ELSEIF ( n_rk .EQ. 3 ) THEN

       ! Tableau for the IMEX-SSP(3,3,2) Stiffly Accurate Scheme
       ! from Pareschi & Russo (2005), Table IV

       a_tilde_ij(2,1) = 0.5d0
       a_tilde_ij(3,1) = 0.5d0
       a_tilde_ij(3,2) = 0.5d0

       omega_tilde(1) =  1.0d0 / 3.0d0
       omega_tilde(2) =  1.0d0 / 3.0d0
       omega_tilde(3) =  1.0d0 / 3.0d0

       a_dirk_ij(1,1) = 0.25d0
       a_dirk_ij(2,2) = 0.25d0
       a_dirk_ij(3,1) = 1.0d0 / 3.0d0
       a_dirk_ij(3,2) = 1.0d0 / 3.0d0
       a_dirk_ij(3,3) = 1.0d0 / 3.0d0

       omega(1) =  1.0d0 / 3.0d0
       omega(2) =  1.0d0 / 3.0d0
       omega(3) =  1.0d0 / 3.0d0


       
    ELSEIF ( n_rk .EQ. 4 ) THEN

       ! LRR(3,2,2) from Table 3 in Pareschi & Russo (2000)

       a_tilde_ij(2,1) = 0.5d0
       a_tilde_ij(3,1) = 1.d0 / 3.d0
       a_tilde_ij(4,2) = 1.0d0

       omega_tilde(1) = 0.d0
       omega_tilde(2) = 1.0d0
       omega_tilde(3) = 0.0d0
       omega_tilde(4) = 0.d0

       a_dirk_ij(2,2) = 0.5d0
       a_dirk_ij(3,3) = 1.0d0 / 3.0d0
       a_dirk_ij(4,3) = 0.75d0
       a_dirk_ij(4,4) = 0.25d0

       omega(1) = 0.d0
       omega(2) = 0.d0
       omega(3) = 0.75d0
       omega(4) = 0.25d0

    END IF

    ALLOCATE( a_tilde(n_rk) )
    ALLOCATE( a_dirk(n_rk) )

    ALLOCATE( q_rk( n_vars , comp_cells_x , comp_cells_y , n_rk ) )
    ALLOCATE( divflux( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
    ALLOCATE( nh( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
    ALLOCATE( si_nh( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )

    ALLOCATE( expl_terms( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )

    ALLOCATE( divfluxj(n_eqns,n_rk) )
    ALLOCATE( nhj(n_eqns,n_rk) )
    ALLOCATE( si_nhj(n_eqns,n_rk) )
    ALLOCATE( expl_terms_j(n_eqns,n_rk) )

    ALLOCATE( residual_term( n_vars , comp_cells_x , comp_cells_y ) )

  END SUBROUTINE allocate_solver_variables

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

  SUBROUTINE deallocate_solver_variables

    DEALLOCATE( q , q0 )

    DEALLOCATE( q_fv )

    DEALLOCATE( q_interfacel )
    DEALLOCATE( q_interfacer )
    DEALLOCATE( q_interfaceb )
    DEALLOCATE( q_interfacet )

    DEALLOCATE( a_interface_xneg )
    DEALLOCATE( a_interface_xpos )
    DEALLOCATE( a_interface_yneg )
    DEALLOCATE( a_interface_ypos )

    DEALLOCATE( h_interface_x )
    DEALLOCATE( h_interface_y )

    DEALLOCATE( solve_mask , solve_mask0 )

    Deallocate( qp )

    DEALLOCATE( source_xy )
    
    DEALLOCATE( a_tilde_ij )
    DEALLOCATE( a_dirk_ij )
    DEALLOCATE( omega_tilde )
    DEALLOCATE( omega )

    DEALLOCATE( implicit_flag )

    DEALLOCATE( a_tilde )
    DEALLOCATE( a_dirk )

    DEALLOCATE( q_rk )
    DEALLOCATE( divflux )
    DEALLOCATE( nh )
    DEALLOCATE( si_nh )
    DEALLOCATE( expl_terms )

    DEALLOCATE( divfluxj )
    DEALLOCATE( nhj )
    DEALLOCATE( si_nhj )
    DEALLOCATE( expl_terms_j )

    DEALLOCATE( mask22 , mask21 , mask11 , mask12 )

    DEALLOCATE( residual_term )

  END SUBROUTINE deallocate_solver_variables


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

  SUBROUTINE check_solve

    IMPLICIT NONE

    INTEGER :: i

    solve_mask0(1:comp_cells_x,1:comp_cells_y) = .false.

    solve_mask = solve_mask0

    WHERE ( q(1,:,:) - b_cent(:,:) .GT. 0.d0 ) solve_mask = .true.

    DO i = 1,n_rk

       solve_mask(1+i:comp_cells_x,:) =  solve_mask(1+i:comp_cells_x,:) .OR.    &
            solve_mask(1:comp_cells_x-i,:) 

       solve_mask(1:comp_cells_x-i,:) =  solve_mask(1:comp_cells_x-i,:) .OR.    &
            solve_mask(1+i:comp_cells_x,:) 

       solve_mask(:,1+i:comp_cells_y) =  solve_mask(:,1+i:comp_cells_y) .OR.    &
            solve_mask(:,1:comp_cells_y-i) 

       solve_mask(:,1:comp_cells_y-i) =  solve_mask(:,1:comp_cells_y-i) .OR.    &
            solve_mask(:,1+i:comp_cells_y) 

    END DO


  END SUBROUTINE check_solve

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

  SUBROUTINE timestep

    ! External variables
    USE geometry_2d, ONLY : dx,dy
    USE parameters_2d, ONLY : max_dt , cfl

    ! External procedures
    USE constitutive_2d, ONLY : eval_local_speeds_x, eval_local_speeds_y

    IMPLICIT none

    REAL*8 :: vel_max(n_vars)
    REAL*8 :: vel_min(n_vars)
    REAL*8 :: vel_j
    REAL*8 :: dt_cfl
    REAL*8 :: qj(n_vars)

    INTEGER :: j,k

    dt = max_dt

    IF ( cfl .NE. -1.d0 ) THEN

       DO j = 1,comp_cells_x

          DO k = 1,comp_cells_y

             qj = q( 1:n_vars , j , k )

             IF ( comp_cells_x .GT. 1 ) THEN

                ! x direction
                CALL eval_local_speeds_x( qj , b_cent(j,k) , vel_min , vel_max )
                
                vel_j = max( maxval(abs(vel_min)) , maxval(abs(vel_max)) )
                
                dt_cfl = cfl * dx / vel_j
                
                dt = min( dt , dt_cfl )

             END IF

             IF ( comp_cells_y .GT. 1 ) THEN
                
                ! y direction
                CALL eval_local_speeds_y( qj , b_cent(j,k) , vel_min , vel_max )
                
                vel_j = max( maxval(abs(vel_min)) , maxval(abs(vel_max)) )
                
                dt_cfl = cfl * dy / vel_j
                                
                dt = min( dt , dt_cfl )

             END IF
                
          ENDDO

       END DO

    END IF

  END SUBROUTINE timestep


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

  SUBROUTINE timestep2

    ! External variables
    USE geometry_2d, ONLY : dx,dy
    USE parameters_2d, ONLY : max_dt , cfl

    ! External procedures
    USE constitutive_2d, ONLY : eval_local_speeds2_x, eval_local_speeds2_y

    IMPLICIT none

    REAL*8 :: dt_cfl

    REAL*8 :: a_interface_x_max(n_eqns,comp_interfaces_x,comp_cells_y)
    REAL*8 :: a_interface_y_max(n_eqns,comp_cells_x,comp_interfaces_y)
    REAL*8 :: dt_interface_x, dt_interface_y

    INTEGER :: i,j,k

    dt = max_dt

    IF ( cfl .NE. -1.d0 ) THEN

       CALL reconstruction

       CALL eval_speeds

       DO i=1,n_vars

          a_interface_x_max(i,:,:) =                                            &
               max( a_interface_xpos(i,:,:) , -a_interface_xneg(i,:,:) )

          a_interface_y_max(i,:,:) =                                            &
               max( a_interface_ypos(i,:,:) , -a_interface_yneg(i,:,:) )

       END DO

       DO j = 1,comp_cells_x

          DO k = 1,comp_cells_y

             dt_interface_x = cfl * dx / max( maxval(a_interface_x_max(:,j,k))  &
                  , maxval(a_interface_x_max(:,j+1,k)) )

             dt_interface_y = cfl * dy / max( maxval(a_interface_y_max(:,j,k))  &
                  , maxval(a_interface_y_max(:,j,k+1)) )

             dt_cfl = min( dt_interface_x , dt_interface_y )

             dt = min(dt,dt_cfl)

          ENDDO

       END DO

    END IF

  END SUBROUTINE timestep2

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

  SUBROUTINE imex_rk_solver

    USE constitutive_2d, ONLY : eval_nonhyperbolic_terms

    USE constitutive_2d, ONLY : eval_nh_semi_impl_terms

    USE constitutive_2d, ONLY : qc_to_qp

    IMPLICIT NONE

    REAL*8 :: q_si(n_vars)
    REAL*8 :: q_guess(n_vars)
    INTEGER :: j,k

    REAL*8 :: h_new
    
    ! Initialization of the solution guess
    q0( 1:n_vars , 1:comp_cells_x , 1:comp_cells_y ) =                          &
         q( 1:n_vars , 1:comp_cells_x , 1:comp_cells_y )

    IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: beginning'

    ! Initialization of the variables for the Runge-Kutta scheme
    q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0

    divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0

    nh(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0

    si_nh(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0

    expl_terms(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0

    runge_kutta:DO i_rk = 1,n_rk

       IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: i_RK',i_rk

       ! define the explicits coefficients for the i-th step of the Runge-Kutta
       a_tilde = 0.d0
       a_dirk = 0.d0

       ! in the first step of the RK scheme all the coefficients remain to 0
       a_tilde(1:i_rk-1) = a_tilde_ij(i_rk,1:i_rk-1)
       a_dirk(1:i_rk-1) = a_dirk_ij(i_rk,1:i_rk-1)

       ! define the implicit coefficient for the i-th step of the Runge-Kutta
       a_diag = a_dirk_ij(i_rk,i_rk)

       CALL cpu_time(t_imex1)

       loop_over_ycells:DO k = 1,comp_cells_y
       
          loop_over_xcells:DO j = 1,comp_cells_x

             IF ( verbose_level .GE. 2 ) THEN

                WRITE(*,*) 'solver, imex_RK_solver: j',j,k

             END IF

             ! initialize the RK step
             IF ( i_rk .EQ. 1 ) THEN

                ! solution from the previous time step
                q_guess(1:n_vars) = q0( 1:n_vars , j , k) 

             ELSE

                ! solution from the previous RK step
                !q_guess(1:n_vars) = q_rk( 1:n_vars , j , k , MAX(1,i_RK-1) )

             END IF

             ! RK explicit hyperbolic terms for the volume (i,j)
             divfluxj(1:n_eqns,1:n_rk) = divflux( 1:n_eqns , j , k , 1:n_rk )

             ! RK implicit terms for the volume (i,j)
             nhj(1:n_eqns,1:n_rk) = nh( 1:n_eqns , j , k , 1:n_rk )

             ! RK semi-implicit terms for the volume (i,j)
             si_nhj(1:n_eqns,1:n_rk) = si_nh( 1:n_eqns , j , k , 1:n_rk )
             
             ! RK additional explicit terms for the volume (i,j)
             expl_terms_j(1:n_eqns,1:n_rk) = expl_terms( 1:n_eqns,j,k,1:n_rk )

             ! New solution at the i_RK step without the implicit  and
             ! semi-implicit term
             q_fv( 1:n_vars , j , k ) = q0( 1:n_vars , j , k )                  &
                  - dt * (matmul( divfluxj(1:n_eqns,1:i_rk)                     &
                  + expl_terms_j(1:n_eqns,1:i_rk) , a_tilde(1:i_rk) )           &
                  - matmul( nhj(1:n_eqns,1:i_rk) + si_nhj(1:n_eqns,1:i_rk) ,    &
                  a_dirk(1:i_rk) ) )
            
             IF ( verbose_level .GE. 2 ) THEN

                WRITE(*,*) 'q_guess',q_guess
                CALL qc_to_qp( q_guess , b_cent(j,k) , qp(1:n_vars,j,k) )
                WRITE(*,*) 'q_guess: qp',qp(1:n_vars,j,k)

             END IF

             adiag_pos:IF ( a_diag .NE. 0.d0 ) THEN

                pos_thick:IF ( q_fv(1,j,k) .GT.  b_cent(j,k) )  THEN

                   ! Eval the semi-implicit discontinuous terms
                   CALL eval_nh_semi_impl_terms( b_cent(j,k) , grav_surf(j,k) , &
                        r_qj = q_fv( 1:n_vars , j , k ) ,                       &
                        r_nh_semi_impl_term = si_nh(1:n_eqns,j,k,i_rk) ) 

                   si_nhj(1:n_eqns,i_rk) = si_nh( 1:n_eqns,j,k,i_rk )

                   ! Assemble the initial guess for the implicit solver
                   q_si(1:n_vars) = q_fv(1:n_vars,j,k ) + dt * a_diag *         &
                        si_nh(1:n_eqns,j,k,i_rk)
                   
                   IF ( q_fv(2,j,k)**2 + q_fv(3,j,k)**2 .EQ. 0.d0 ) THEN
                      
                      !Case 1: if the velocity was null, then it must stay null
                      q_si(2:3) = 0.d0 
                      
                   ELSEIF ( ( q_si(2)*q_fv(2,j,k) .LT. 0.d0 ) .OR.              &
                        ( q_si(3)*q_fv(3,j,k) .LT. 0.d0 ) ) THEN
                      
                      ! If the semi-impl. friction term changed the sign of the 
                      ! velocity then set it to zero
                      q_si(2:3) = 0.d0 

                   ELSE
                      
                      ! Align the velocity vector with previous one
                      q_si(2:3) = dsqrt( q_si(2)**2 + q_si(3)**2 ) *            &
                           q_fv(2:3,j,k) / dsqrt( q_fv(2,j,k)**2                &
                           + q_fv(3,j,k)**2 ) 
                      
                   END IF

                   ! Update the semi-implicit term accordingly with the
                   ! corrections above
                   si_nh(1:n_eqns,j,k,i_rk) = ( q_si(1:n_vars) -                &
                        q_fv(1:n_vars,j,k ) ) / ( dt*a_diag )

                   si_nhj(1:n_eqns,i_rk) = si_nh( 1:n_eqns,j,k,i_rk )

                   ! Initialize the guess for the NR solver
                   q_guess(1:n_vars) = q_si(1:n_vars)
                   
                   ! Solve the implicit system to find the solution at the 
                   ! i_RK step of the IMEX RK procedure
                   CALL solve_rk_step( b_cent(j,k) , b_prime_x(j,k) ,           &
                        b_prime_y(j,k), grav_surf(j,k) ,                        &
                        q_guess(1:n_vars) , q0(1:n_vars,j,k ) , a_tilde ,       &
                        a_dirk , a_diag )

                   IF ( comp_cells_y .EQ. 1 ) THEN

                      q_guess(3) = 0.d0

                   END IF

                   IF ( comp_cells_x .EQ. 1 ) THEN

                      q_guess(2) = 0.d0

                   END IF

                   ! Eval and store the implicit term at the i_RK step
                   CALL eval_nonhyperbolic_terms( b_cent(j,k) , b_prime_x(j,k) ,&
                        b_prime_y(j,k) , grav_surf(j,k) , r_qj = q_guess ,      &
                        r_nh_term_impl = nh(1:n_eqns,j,k,i_rk) )
                   
                   IF ( q_si(2)**2 + q_si(3)**2 .EQ. 0.d0 ) THEN
                      
                      q_guess(2:3) = 0.d0 
                      
                   ELSEIF ( ( q_guess(2)*q_si(2) .LE. 0.d0 ) .AND.              &
                        ( q_guess(3)*q_si(3) .LE. 0.d0 ) ) THEN
                      
                      ! If the impl. friction term changed the sign of the 
                      ! velocity then set it to zero
                      q_guess(2:3) = 0.d0 
                      
                   ELSE
                      
                      ! Align the velocity vector with previous one
                      q_guess(2:3) = dsqrt( q_guess(2)**2 + q_guess(3)**2 ) *   &
                           q_si(2:3) / dsqrt( q_si(2)**2 + q_si(3)**2 ) 
                                            
                   END IF
                                      
                ELSE

                   ! If h=0 nothing has to be changed 
                   q_guess(1:n_vars) = q_fv( 1:n_vars , j , k ) 
                   q_si(1:n_vars) = q_fv( 1:n_vars , j , k ) 
                   si_nh(1:n_eqns,j,k,i_rk) = 0.d0
                   nh(1:n_eqns,j,k,i_rk) = 0.d0

                END IF pos_thick

             END IF adiag_pos

             ! Check the sign of the flow thickness
             h_new = q_guess(1) - b_cent(j,k)

             IF ( j .EQ. 0 ) THEN
                
                WRITE(*,*) 'h',q_guess(1) -  b_cent(j,k)
                READ(*,*)

             END IF
             
             ! IF ( h_new .LT. 1.D-10 ) THEN
             !
             !   q_guess(1) = B_cent(j,k)
             !   q_guess(2:n_vars) = 0.D0
             !      
             ! END IF

             IF ( a_diag .NE. 0.d0 ) THEN
                
                ! Update the viscous term with the correction on the new velocity
                nh(1:n_vars,j,k,i_rk) = ( q_guess(1:n_vars) - q_si(1:n_vars))   &
                     / ( dt*a_diag ) 

             END IF

             ! Store the solution at the end of the i_RK step
             q_rk( 1:n_vars , j , k , i_rk ) = q_guess

             IF ( verbose_level .GE. 2 ) THEN

                WRITE(*,*) 'imex_RK_solver: qc',q_guess
                CALL qc_to_qp( q_guess, b_cent(j,k) , qp(1:n_vars,j,k) )
                WRITE(*,*) 'imex_RK_solver: qp',qp(1:n_vars,j,k)
                READ(*,*)

             END IF

          END DO loop_over_xcells

       ENDDO loop_over_ycells

       CALL cpu_time(t_imex2)
       ! WRITE(*,*) 'Time taken by implicit',t_imex2-t_imex1,'seconds'


       ! Eval and store the explicit hyperbolic (fluxes) terms
       CALL eval_hyperbolic_terms( q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y, &
            i_rk) , divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )

       CALL cpu_time(t_imex1)
       ! WRITE(*,*) 'Time taken by explicit',t_imex1-t_imex2,'seconds'

       ! Eval and store the other explicit terms (e.g. gravity or viscous forces)
       CALL eval_explicit_terms( q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y,   &
            i_rk) , expl_terms(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )


       CALL cpu_time(t_imex1)
       ! WRITE(*,*) 'Time taken by explicit',t_imex1-t_imex2,'seconds'

       IF ( verbose_level .GE. 1 ) THEN

          WRITE(*,*) 'div_flux(2),div_flux(3),expl_terms(2),expl_terms(3)'

          DO k = 1,comp_cells_y
          
             DO j = 1,comp_cells_x

                WRITE(*,*) divflux(2,j,k,i_rk) , divflux(3,j,k,i_rk) ,          &
                     expl_terms(2,j,k,i_rk) , expl_terms(3,j,k,i_rk)

             ENDDO

          END DO

          READ(*,*)

       END IF

    END DO runge_kutta

    DO k = 1,comp_cells_y
    
       DO j = 1,comp_cells_x

             residual_term(1:n_vars,j,k) = matmul( divflux(1:n_eqns,j,k,1:n_rk) &
                  + expl_terms(1:n_eqns,j,k,1:n_rk) , omega_tilde ) -           &
                  matmul( nh(1:n_eqns,j,k,1:n_rk) + si_nh(1:n_eqns,j,k,1:n_rk) ,&
                  omega )

       ENDDO

    END DO

    DO k = 1,comp_cells_y
    
       DO j = 1,comp_cells_x

          IF ( verbose_level .GE. 1 ) THEN

             WRITE(*,*) 'cell jk =',j,k
             WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
             CALL qc_to_qp(q0(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
             WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars,j,k)

          END IF

          IF ( ( sum(abs( omega_tilde(:)-a_tilde_ij(n_rk,:))) .EQ. 0.d0  )      &
               .AND. ( sum(abs(omega(:)-a_dirk_ij(n_rk,:))) .EQ. 0.d0 ) ) THEN

             ! The assembling coeffs are equal to the last step of the RK scheme
             q(1:n_vars,j,k) = q_rk(1:n_vars,j,k,n_rk)

          ELSE

             ! The assembling coeffs are different
             q(1:n_vars,j,k) = q0(1:n_vars,j,k) - dt*residual_term(1:n_vars,j,k)

          END IF


          IF ( q(1,j,k) .LT. b_cent(j,k) ) THEN

             IF ( dabs( q(1,j,k)-b_cent(j,k) ) .LT. 1.d-10 ) THEN
                
                q(1,j,k) = b_cent(j,k)
                q(2:n_vars,j,k) = 0.d0
                
             ELSE
                
                WRITE(*,*) 'j,k,n_RK',j,k,n_rk
                WRITE(*,*) 'dt',dt
                
                WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
                CALL qc_to_qp(q0(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
                WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars,j,k)
                WRITE(*,*) 'h old',q0(1,j,k) - b_cent(j,k)
               
                READ(*,*)
                
             END IF

          END IF

          IF ( verbose_level .GE. 1 ) THEN
          
             CALL qc_to_qp(q(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
             
             WRITE(*,*) 'after imex_RK_solver: qc',q(1:n_vars,j,k)
             WRITE(*,*) 'after imex_RK_solver: qp',qp(1:n_vars,j,k)
             WRITE(*,*) 'h new',q(1,j,k) - b_cent(j,k)
             READ(*,*)

          END IF


       ENDDO

    END DO

  END SUBROUTINE imex_rk_solver

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

  SUBROUTINE solve_rk_step( Bj, Bprimej_x, Bprimej_y, grav3_surf,               &
       qj, qj_old, a_tilde , a_dirk , a_diag )

    USE parameters_2d, ONLY : max_nl_iter , tol_rel , tol_abs

    USE constitutive_2d, ONLY : qc_to_qp

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bj
    REAL*8, INTENT(IN) :: Bprimej_x
    REAL*8, INTENT(IN) :: Bprimej_y
    REAL*8, INTENT(IN) :: grav3_surf
    REAL*8, INTENT(INOUT) :: qj(n_vars)
    REAL*8, INTENT(IN) :: qj_old(n_vars)
    REAL*8, INTENT(IN) :: a_tilde(n_rk)
    REAL*8, INTENT(IN) :: a_dirk(n_rk)
    REAL*8, INTENT(IN) :: a_diag

    REAL*8 :: qj_init(n_vars)

    REAL*8 :: qj_org(n_vars) , qj_rel(n_vars)

    REAL*8 :: left_matrix(n_eqns,n_vars)
    REAL*8 :: right_term(n_eqns)

    REAL*8 :: scal_f

    REAL*8 :: coeff_f(n_eqns)

    REAL*8 :: qj_rel_NR_old(n_vars)
    REAL*8 :: scal_f_old
    REAL*8 :: desc_dir(n_vars)
    REAL*8 :: grad_f(n_vars)
    REAL*8 :: mod_desc_dir

    INTEGER :: pivot(n_vars)

    REAL*8 :: left_matrix_small22(n_nh,n_nh)
    REAL*8 :: left_matrix_small21(n_eqns-n_nh,n_nh)
    REAL*8 :: left_matrix_small11(n_eqns-n_nh,n_vars-n_nh)
    REAL*8 :: left_matrix_small12(n_nh,n_vars-n_nh)

    REAL*8 :: desc_dir_small2(n_nh)
    INTEGER :: pivot_small2(n_nh)

    REAL*8 :: desc_dir_small1(n_vars-n_nh)

    INTEGER :: ok

    INTEGER :: i 
    INTEGER :: nl_iter

    REAL*8, PARAMETER :: STPMX=100.d0
    REAL*8 :: stpmax
    LOGICAL :: check

    REAL*8, PARAMETER :: TOLF=1.d-10 , tolmin=1.d-6
    REAL*8 :: TOLX

    REAL*8 :: qpj(n_vars)

    REAL*8 :: desc_dir2(n_vars)

    REAL*8 :: desc_dir_temp(n_vars)

    normalize_q = .true.
    normalize_f = .false.
    opt_search_nl = .true.

    coeff_f(1:n_eqns) = 1.d0

    qj_init = qj

    ! normalize the functions of the nonlinear system
    IF ( normalize_f ) THEN

       qj = qj_old - dt * ( matmul(divfluxj+ expl_terms_j,a_tilde)              &
            - matmul(nhj,a_dirk) )

       CALL eval_f( bj , bprimej_x , bprimej_y , grav3_surf ,                   &
            qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , right_term ,    &
            scal_f )

       IF ( verbose_level .GE. 3 ) THEN

          WRITE(*,*) 'solve_rk_step: non-normalized right_term'
          WRITE(*,*) right_term
          WRITE(*,*) 'scal_f',scal_f

       END IF

       DO i=1,n_eqns

          IF ( abs(right_term(i)) .GE. 1.d0 ) coeff_f(i) = 1.d0 / right_term(i)

       END DO

       right_term = coeff_f * right_term

       scal_f = 0.5d0 * dot_product( right_term , right_term )

       IF ( verbose_level .GE. 3 ) THEN                    
          WRITE(*,*) 'solve_rk_step: after normalization',scal_f
       END IF

    END IF

    !---- normalize the conservative variables ------

    IF ( normalize_q ) THEN

       qj_org = qj

       qj_org = max( abs(qj_org) , 1.d-3 )

    ELSE 

       qj_org(1:n_vars) = 1.d0

    END IF

    qj_rel = qj / qj_org

    ! -----------------------------------------------

    newton_raphson_loop:DO nl_iter=1,max_nl_iter

       tolx = epsilon(qj_rel)

       IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solve_rk_step: nl_iter',nl_iter

       CALL eval_f( bj , bprimej_x , bprimej_y , grav3_surf ,                   &
            qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , right_term ,    &
            scal_f )

       IF ( verbose_level .GE. 2 ) THEN

          WRITE(*,*) 'solve_rk_step: right_term',right_term

       END IF

       IF ( verbose_level .GE. 2 ) THEN

          WRITE(*,*) 'before_lnsrch: scal_f',scal_f

       END IF

       ! check the residual of the system

       IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN

          IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
          RETURN

       END IF

       IF ( ( normalize_f ) .AND. ( scal_f < 1.d-6 ) ) THEN

          IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check scal_f',check
          RETURN

       END IF

       ! ---- evaluate the descent direction ------------------------------------

       IF ( count( implicit_flag ) .EQ. n_eqns ) THEN

          CALL eval_jacobian(bj,bprimej_x,bprimej_y,grav3_surf,                 &
               qj_rel,qj_org,coeff_f,left_matrix)

          desc_dir_temp = - right_term

          CALL dgesv(n_eqns,1, left_matrix , n_eqns, pivot, desc_dir_temp ,     &
               n_eqns, ok)

          desc_dir = desc_dir_temp

       ELSE

          CALL eval_jacobian(bj,bprimej_x,bprimej_y,grav3_surf,                 &
               qj_rel,qj_org,coeff_f,left_matrix)

          left_matrix_small11 = reshape(pack(left_matrix, mask11),              &
               [n_eqns-n_nh,n_eqns-n_nh]) 

          left_matrix_small12 = reshape(pack(left_matrix, mask12),              &
               [n_nh,n_eqns-n_nh]) 

          left_matrix_small22 = reshape(pack(left_matrix, mask22),              &
               [n_nh,n_nh]) 

          left_matrix_small21 = reshape(pack(left_matrix, mask21),              &
               [n_eqns-n_nh,n_nh]) 


          desc_dir_small1 = pack( right_term, .NOT.implicit_flag )
          desc_dir_small2 = pack( right_term , implicit_flag )

          DO i=1,n_vars-n_nh

             desc_dir_small1(i) = desc_dir_small1(i) / left_matrix_small11(i,i)

          END DO

          desc_dir_small2 = desc_dir_small2 -                                   &
               matmul( desc_dir_small1 , left_matrix_small21 )

          CALL dgesv(n_nh,1, left_matrix_small22 , n_nh , pivot_small2 ,        &
               desc_dir_small2 , n_nh, ok)

          desc_dir = unpack( - desc_dir_small2 , implicit_flag , 0.0d0 )        &
               + unpack( - desc_dir_small1 , .NOT.implicit_flag , 0.0d0 )

       END IF

       mod_desc_dir = dsqrt( desc_dir(2)**2 + desc_dir(3)**2 )

       !IF (  qj(2)**2 + qj(3)**2 .GT. 0.D0 ) THEN 
       !
       !   desc_dir(2) = mod_desc_dir * qj(2) / ( qj(2)**2 + qj(3)**2 ) 
       !   desc_dir(3) = mod_desc_dir * qj(3) / ( qj(2)**2 + qj(3)**2 ) 
       !
       !END IF

       IF ( verbose_level .GE. 3 ) WRITE(*,*) 'desc_dir',desc_dir

       qj_rel_nr_old = qj_rel
       scal_f_old = scal_f

       IF ( ( opt_search_nl ) .AND. ( nl_iter .GT. 1 ) ) THEN
          ! Search for the step lambda giving a suffic. decrease in the solution 

          stpmax = stpmx * max( sqrt( dot_product(qj_rel,qj_rel) ) ,            &
               dble( SIZE(qj_rel) ) )

          grad_f = matmul( right_term , left_matrix )

          desc_dir2 = desc_dir

          CALL lnsrch( bj , bprimej_x , bprimej_y , grav3_surf ,                &
               qj_rel_nr_old , qj_org , qj_old , scal_f_old , grad_f ,          &
               desc_dir , coeff_f , qj_rel , scal_f , right_term , stpmax ,     &
               check )

       ELSE

          qj_rel = qj_rel_nr_old + desc_dir

          qj = qj_rel * qj_org

          CALL eval_f( bj , bprimej_x , bprimej_y , grav3_surf ,                &
               qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f ,              &
               right_term , scal_f )

       END IF

       IF ( verbose_level .GE. 2 ) WRITE(*,*) 'after_lnsrch: scal_f',scal_f

       qj = qj_rel * qj_org

       IF ( verbose_level .GE. 3 ) THEN

          WRITE(*,*) 'qj',qj
          CALL qc_to_qp( qj , bj , qpj)
          WRITE(*,*) 'qp',qpj

       END IF

       
       IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN

          IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
          check= .false.
          RETURN

       END IF

       IF (check) THEN

          check = ( maxval( abs(grad_f(:)) * max( abs( qj_rel(:) ),1.d0 ) /     &
               max( scal_f , 0.5d0 * SIZE(qj_rel) ) )  < tolmin )

          IF ( verbose_level .GE. 3 ) WRITE(*,*) '2: check',check
          !          RETURN

       END IF

       IF ( maxval( abs( qj_rel(:) - qj_rel_nr_old(:) ) / max( abs( qj_rel(:)) ,&
            1.d0 ) ) < tolx ) THEN

          IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check',check
          RETURN

       END IF

    END DO newton_raphson_loop

  END SUBROUTINE solve_rk_step

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

  SUBROUTINE lnsrch( Bj , Bprimej_x , Bprimej_y , grav3_surf ,                  &
       qj_rel_nr_old , qj_org , qj_old , scal_f_old , grad_f ,                  &
       desc_dir , coeff_f , qj_rel , scal_f , right_term , stpmax , check )

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bj

    REAL*8, INTENT(IN) :: Bprimej_x

    REAL*8, INTENT(IN) :: Bprimej_y

    REAL*8, INTENT(IN) :: grav3_surf

    REAL*8, DIMENSION(:), INTENT(IN) :: qj_rel_NR_old

    REAL*8, DIMENSION(:), INTENT(IN) :: qj_org

    REAL*8, DIMENSION(:), INTENT(IN) :: qj_old

    REAL*8, DIMENSION(:), INTENT(IN) :: grad_f

    REAL*8, INTENT(IN) :: scal_f_old

    REAL*8, DIMENSION(:), INTENT(INOUT) :: desc_dir

    REAL*8, INTENT(IN) :: stpmax

    REAL*8, DIMENSION(:), INTENT(IN) :: coeff_f

    REAL*8, DIMENSION(:), INTENT(OUT) :: qj_rel

    REAL*8, INTENT(OUT) :: scal_f

    REAL*8, INTENT(OUT) :: right_term(n_eqns)

    LOGICAL, INTENT(OUT) :: check

    REAL*8, PARAMETER :: TOLX=epsilon(qj_rel)

    INTEGER, DIMENSION(1) :: ndum
    REAL*8 :: ALF , a,alam,alam2,alamin,b,disc
    REAL*8 :: scal_f2
    REAL*8 :: desc_dir_abs
    REAL*8 :: rhs1 , rhs2 , slope, tmplam

    REAL*8 :: scal_f_min , alam_min

    REAL*8 :: qj(n_vars)

    alf = 1.0d-4

    IF ( size(grad_f) == size(desc_dir) .AND. size(grad_f) == size(qj_rel)      &
         .AND. size(qj_rel) == size(qj_rel_nr_old) ) THEN

       ndum = size(grad_f)

    ELSE

       WRITE(*,*) 'nrerror: an assert_eq failed with this tag:', 'lnsrch'
       stop 'program terminated by assert_eq4'

    END IF

    check = .false.

    desc_dir_abs = sqrt( dot_product(desc_dir,desc_dir) )

    IF ( desc_dir_abs > stpmax ) desc_dir(:) = desc_dir(:) * stpmax/desc_dir_abs  

    slope = dot_product(grad_f,desc_dir)

    alamin = tolx / maxval( abs( desc_dir(:))/max( abs(qj_rel_nr_old(:)),1.d0 ) )

    IF ( alamin .EQ. 0.d0) THEN

       qj_rel(:) = qj_rel_nr_old(:)

       RETURN

    END IF

    alam = 1.0d0

    scal_f_min = scal_f_old

    optimal_step_search: DO

       IF ( verbose_level .GE. 4 ) THEN

          WRITE(*,*) 'alam',alam

       END IF

       qj_rel = qj_rel_nr_old + alam * desc_dir

       qj = qj_rel * qj_org

       CALL eval_f( bj , bprimej_x , bprimej_y , grav3_surf , qj , qj_old ,     &
            a_tilde , a_dirk , a_diag , coeff_f , right_term , scal_f )

       IF ( verbose_level .GE. 4 ) THEN

          WRITE(*,*) 'lnsrch: effe_old,effe',scal_f_old,scal_f
          READ(*,*)

       END IF

       IF ( scal_f .LT. scal_f_min ) THEN

          scal_f_min = scal_f
          alam_min = alam

       END IF

       IF ( scal_f .LE. 0.9 * scal_f_old ) THEN   
          ! sufficient function decrease

          IF ( verbose_level .GE. 4 ) THEN

             WRITE(*,*) 'sufficient function decrease'

          END IF

          EXIT optimal_step_search   

       ELSE IF ( alam < alamin ) THEN   
          ! convergence on Delta_x

          IF ( verbose_level .GE. 4 ) THEN

             WRITE(*,*) ' convergence on Delta_x',alam,alamin

          END IF

          qj_rel(:) = qj_rel_nr_old(:)
          scal_f = scal_f_old
          check = .true.

          EXIT optimal_step_search

          !       ELSE IF ( scal_f .LE. scal_f_old + ALF * alam * slope ) THEN   
       ELSE  

          IF ( alam .EQ. 1.d0 ) THEN

             tmplam = - slope / ( 2.0d0 * ( scal_f - scal_f_old - slope ) )

          ELSE

             rhs1 = scal_f - scal_f_old - alam*slope
             rhs2 = scal_f2 - scal_f_old - alam2*slope

             a = ( rhs1/alam**2.d0 - rhs2/alam2**2.d0 ) / ( alam - alam2 )
             b = ( -alam2*rhs1/alam**2 + alam*rhs2/alam2**2 ) / ( alam - alam2 )

             IF ( a .EQ. 0.d0 ) THEN

                tmplam = - slope / ( 2.0d0 * b )

             ELSE

                disc = b*b - 3.0d0*a*slope

                IF ( disc .LT. 0.d0 ) THEN

                   tmplam = 0.5d0 * alam

                ELSE IF ( b .LE. 0.d0 ) THEN

                   tmplam = ( - b + sqrt(disc) ) / ( 3.d0 * a )

                ELSE

                   tmplam = - slope / ( b + sqrt(disc) )

                ENDIF

             END IF

             IF ( tmplam .GT. 0.5d0*alam ) tmplam = 0.5d0 * alam

          END IF

       END IF

       alam2 = alam
       scal_f2 = scal_f
       alam = max( tmplam , 0.5d0*alam )

    END DO optimal_step_search

  END SUBROUTINE lnsrch

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

  SUBROUTINE eval_f( Bj , Bprimej_x , Bprimej_y , grav3_surf , qj , qj_old ,    &
       a_tilde , a_dirk , a_diag , coeff_f , f_nl , scal_f )

    USE constitutive_2d, ONLY : eval_nonhyperbolic_terms

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bj
    REAL*8, INTENT(IN) :: Bprimej_x
    REAL*8, INTENT(IN) :: Bprimej_y
    REAL*8, INTENT(IN) :: grav3_surf
    REAL*8, INTENT(IN) :: qj(n_vars)
    REAL*8, INTENT(IN) :: qj_old(n_vars)
    REAL*8, INTENT(IN) :: a_tilde(n_rk)
    REAL*8, INTENT(IN) :: a_dirk(n_rk)
    REAL*8, INTENT(IN) :: a_diag
    REAL*8, INTENT(IN) :: coeff_f(n_eqns)
    REAL*8, INTENT(OUT) :: f_nl(n_eqns)
    REAL*8, INTENT(OUT) :: scal_f

    REAL*8 :: nh_term_impl(n_eqns)
    REAL*8 :: Rj(n_eqns)

    CALL eval_nonhyperbolic_terms( bj , bprimej_x , bprimej_y , grav3_surf ,    &
         r_qj = qj , r_nh_term_impl = nh_term_impl ) 

    rj = ( matmul( divfluxj(1:n_eqns,1:i_rk-1)+expl_terms_j(1:n_eqns,1:i_rk-1), &
         a_tilde(1:i_rk-1) ) - matmul( nhj(1:n_eqns,1:i_rk-1)                   &
         + si_nhj(1:n_eqns,1:i_rk-1) , a_dirk(1:i_rk-1) ) )                     &
         - a_diag * ( nh_term_impl + si_nhj(1:n_eqns,i_rk) )

    f_nl = qj - qj_old + dt * rj

    f_nl = coeff_f * f_nl

    scal_f = 0.5d0 * dot_product( f_nl , f_nl )

    !WRITE(*,*) 'i_RK',i_RK
    !WRITE(*,*) 'eval_f: qj(2) = ',qj(2)
    !WRITE(*,*) 'nh_term_impl(2) = ',nh_term_impl(2)
    !WRITE(*,*) ' SI_NHj(1:n_eqns,i_RK) = ', SI_NHj(1:n_eqns,i_RK)
    !WRITE(*,*) 'eval_f: new term = ',qj(2) + dt * ( - a_diag * nh_term_impl(2) )
    !WRITE(*,*) 'f_nl(2)=',f_nl(2)
    !READ(*,*)

  END SUBROUTINE eval_f

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

  SUBROUTINE eval_jacobian(Bj , Bprimej_x , Bprimej_y , grav3_surf ,            &
       qj_rel , qj_org , coeff_f, left_matrix)

    USE constitutive_2d, ONLY : eval_nonhyperbolic_terms

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bj
    REAL*8, INTENT(IN) :: Bprimej_x
    REAL*8, INTENT(IN) :: Bprimej_y
    REAL*8, INTENT(IN) :: grav3_surf
    REAL*8, INTENT(IN) :: qj_rel(n_vars)
    REAL*8, INTENT(IN) :: qj_org(n_vars)
    REAL*8, INTENT(IN) :: coeff_f(n_eqns)
    REAL*8, INTENT(OUT) :: left_matrix(n_eqns,n_vars)

    REAL*8 :: Jacob_relax(n_eqns,n_vars)
    COMPLEX*16 :: nh_terms_cmplx_impl(n_eqns)
    COMPLEX*16 :: qj_cmplx(n_vars) , qj_rel_cmplx(n_vars)

    REAL*8 :: h

    INTEGER :: i

    h = n_vars * epsilon(1.d0)

    ! initialize the matrix of the linearized system and the Jacobian

    left_matrix(1:n_eqns,1:n_vars) = 0.d0
    jacob_relax(1:n_eqns,1:n_vars) = 0.d0

    ! evaluate the jacobian of the non-hyperbolic terms

    DO i=1,n_vars

       left_matrix(i,i) = coeff_f(i) * qj_org(i)

       IF ( implicit_flag(i) ) THEN 

          qj_rel_cmplx(1:n_vars) = qj_rel(1:n_vars)
          qj_rel_cmplx(i) = dcmplx(qj_rel(i), h)

          qj_cmplx = qj_rel_cmplx * qj_org

          CALL eval_nonhyperbolic_terms( bj , bprimej_x , bprimej_y ,           &
               grav3_surf, c_qj = qj_cmplx ,                                    &
               c_nh_term_impl = nh_terms_cmplx_impl ) 

          jacob_relax(1:n_eqns,i) = coeff_f(i) *                                &
               aimag(nh_terms_cmplx_impl) / h

          left_matrix(1:n_eqns,i) = left_matrix(1:n_eqns,i) - dt * a_diag       &
               * jacob_relax(1:n_eqns,i)

       END IF

    END DO

  END SUBROUTINE eval_jacobian

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

  SUBROUTINE eval_explicit_terms( q_expl , expl_terms )

    USE constitutive_2d, ONLY : eval_expl_terms

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells_x,comp_cells_y)
    REAL*8, INTENT(OUT) :: expl_terms(n_eqns,comp_cells_x,comp_cells_y)

    REAL*8 :: qc(n_vars)
    REAL*8 :: expl_forces_term(n_eqns)

    INTEGER :: j,k

    DO k = 1,comp_cells_y
    
       DO j = 1,comp_cells_x

          qc = q_expl(1:n_vars,j,k)

          CALL eval_expl_terms( b_cent(j,k), b_prime_x(j,k),                    &
               b_prime_y(j,k), source_xy(j,k) , qc, expl_forces_term )

          expl_terms(1:n_eqns,j,k) =  expl_forces_term

       ENDDO

    END DO

  END SUBROUTINE eval_explicit_terms

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

  SUBROUTINE eval_hyperbolic_terms( q_expl , divFlux )

    ! External variables
    USE geometry_2d, ONLY : dx,dy
    USE parameters_2d, ONLY : solver_scheme

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells_x,comp_cells_y)
    REAL*8, INTENT(OUT) :: divFlux(n_eqns,comp_cells_x,comp_cells_y)

    REAL*8 :: q_old(n_vars,comp_cells_x,comp_cells_y)

    REAL*8 :: h_new

    REAL*8 :: tcpu0,tcpu1,tcpu2,tcpu3,tcpu4

    INTEGER :: i, j, k

    q_old = q

    q = q_expl

    CALL cpu_time(tcpu0)

    ! Linear reconstruction of the physical variables at the interfaces
    CALL reconstruction

    CALL cpu_time(tcpu1)
    !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu1-tcpu0,'seconds'


    ! Evaluation of the maximum local speeds at the interfaces
    CALL eval_speeds

    CALL cpu_time(tcpu2)
    !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu2-tcpu1,'seconds'

    ! Evaluation of the numerical fluxes
    SELECT CASE ( solver_scheme )

    CASE ("LxF")

       CALL eval_flux_lxf

    CASE ("GFORCE")

       CALL eval_flux_gforce

    CASE ("KT")

       CALL eval_flux_kt

    END SELECT

    CALL cpu_time(tcpu3)
    !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu3-tcpu2,'seconds'



    ! Advance in time the solution
    DO k = 1,comp_cells_y
    
       DO j = 1,comp_cells_x
  
          DO i=1,n_eqns
             
             divflux(i,j,k) = 0.d0
             
             IF ( comp_cells_x .GT. 1 ) THEN
                
                divflux(i,j,k) = divflux(i,j,k) +                               &
                     ( h_interface_x(i,j+1,k) - h_interface_x(i,j,k) ) / dx
                
             END IF

             IF ( comp_cells_y .GT. 1 ) THEN
                
                divflux(i,j,k) = divflux(i,j,k) +                               &
                     ( h_interface_y(i,j,k+1) - h_interface_y(i,j,k) ) / dy
                
             END IF

             
          END DO

          h_new = q_expl(1,j,k) - dt * divflux(1,j,k) - b_cent(j,k)

          ! IF ( h_new .NE. 0.D0 ) THEN
          IF ( j .EQ. 0 ) THEN

             WRITE(*,*) 'j,k,h,divF',j,k,h_new, divflux(1,j,k)
             WRITE(*,*) 'dt',dt

             WRITE(*,*) 'h_interface(j,k) ',q_interfacel(1,j,k)-b_stag_x(j,k) , &
                  q_interfacer(1,j,k) - b_stag_x(j,k)

             WRITE(*,*) 'hu_interface(j,k)' , q_interfacel(2,j,k) ,             &
                  q_interfacer(2,j,k)

             WRITE(*,*) 'h_interface(j+1,k) ' , q_interfacel(1,j+1,k)           &
                  - b_stag_y(j+1,k) , q_interfacer(1,j+1,k) - b_stag_y(j+1,k)

             WRITE(*,*) 'hu_interface(j+1,k)' , q_interfacel(2,j+1,k) ,         &
                  q_interfacer(2,j+1,k)

             WRITE(*,*) 'h flux x',h_interface_x(1,j,k) , h_interface_x(1,j+1,k)

             WRITE(*,*) 'h flux y',h_interface_y(i,j,k) , h_interface_y(1,j,k+1)

             WRITE(*,*) 'hu flux x',h_interface_x(2,j,k) , h_interface_x(2,j+1,k)

             WRITE(*,*) 'hu flux y',h_interface_y(2,j,k) , h_interface_y(2,j,k+1)

             READ(*,*) 

          END IF

       ENDDO

    END DO

    CALL cpu_time(tcpu4)
    !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu4-tcpu4,'seconds'

    q = q_old

  END SUBROUTINE eval_hyperbolic_terms


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

  SUBROUTINE eval_flux_kt

    ! External procedures
    USE constitutive_2d, ONLY : eval_fluxes

    IMPLICIT NONE

    REAL*8 :: fluxL(n_eqns)
    REAL*8 :: fluxR(n_eqns)
    REAL*8 :: fluxB(n_eqns)
    REAL*8 :: fluxT(n_eqns)

    REAL*8 :: flux_avg_x(n_eqns)   
    REAL*8 :: flux_avg_y(n_eqns)   

    INTEGER :: j,k
    INTEGER :: i

    IF ( comp_cells_x .GT. 1 ) THEN

       DO k = 1,comp_cells_y
       
          DO j = 1,comp_interfaces_x

             CALL eval_fluxes( b_stag_x(j,k) ,                                  &
                  r_qj = q_interfacel(1:n_vars,j,k) , r_flux=fluxl , dir=1 )

             CALL eval_fluxes( b_stag_x(j,k) ,                                  &
                  r_qj = q_interfacer(1:n_vars,j,k) , r_flux=fluxr , dir=1 )

             CALL average_kt( a_interface_xneg(:,j,k), a_interface_xpos(:,j,k) ,&
                  fluxl , fluxr , flux_avg_x )

             eqns_loop:DO i=1,n_eqns

                IF ( a_interface_xneg(i,j,k) .EQ. a_interface_xpos(i,j,k) ) THEN

                   h_interface_x(i,j,k) = 0.d0

                ELSE

                   h_interface_x(i,j,k) = flux_avg_x(i)                         &
                        + ( a_interface_xpos(i,j,k) * a_interface_xneg(i,j,k) ) &
                        / ( a_interface_xpos(i,j,k) - a_interface_xneg(i,j,k) ) &
                        * ( q_interfacer(i,j,k) - q_interfacel(i,j,k) )             

                END IF

             ENDDO eqns_loop

             ! In the equation for mass and for trasnport (T,xs) if the 
             ! velocities at the interfaces are null, then the flux is null
             IF ( (  q_interfacel(2,j,k) .EQ. 0.d0 ) .AND.                      &
                  (  q_interfacer(2,j,k) .EQ. 0.d0 ) ) THEN

                h_interface_x(1,j,k) = 0.d0
                h_interface_x(4:n_vars,j,k) = 0.d0

             END IF

             IF ( j .LE. 0 ) THEN 
                
                WRITE(*,*) 'eval_flux_KT'
                WRITE(*,*) 'j',j
                WRITE(*,*) a_interface_xpos(1,j,k), a_interface_xneg(1,j,k)
                WRITE(*,*) q_interfacer(1,j,k) , q_interfacel(1,j,k)
                READ(*,*)
                
             END IF

          END DO

       END DO

    END IF

    IF ( comp_cells_y .GT. 1 ) THEN

       DO k = 1,comp_interfaces_y
          
          DO j = 1,comp_cells_x

             CALL eval_fluxes( b_stag_y(j,k) ,                                  &
                  r_qj = q_interfaceb(1:n_vars,j,k) , r_flux=fluxb , dir=2 )

             CALL eval_fluxes( b_stag_y(j,k) ,                                  &
                  r_qj = q_interfacet(1:n_vars,j,k) , r_flux=fluxt , dir=2 )

             CALL average_kt( a_interface_yneg(:,j,k) ,                         &
                  a_interface_ypos(:,j,k) , fluxb , fluxt , flux_avg_y )

             DO i=1,n_eqns

                IF ( a_interface_yneg(i,j,k) .EQ. a_interface_ypos(i,j,k) ) THEN

                   h_interface_y(i,j,k) = 0.d0

                ELSE

                   h_interface_y(i,j,k) = flux_avg_y(i)                         &
                        + ( a_interface_ypos(i,j,k) * a_interface_yneg(i,j,k) ) &
                        / ( a_interface_ypos(i,j,k) - a_interface_yneg(i,j,k) ) &
                        * ( q_interfacet(i,j,k) - q_interfaceb(i,j,k) )             

                END IF

             END DO

             ! In the equation for mass and for trasnport (T,xs) if the 
             ! velocities at the interfaces are null, then the flux is null
             IF ( (  q_interfaceb(3,j,k) .EQ. 0.d0 ) .AND.                      &
                  (  q_interfacet(3,j,k) .EQ. 0.d0 ) ) THEN

                h_interface_y(1,j,k) = 0.d0
                h_interface_y(4:n_vars,j,k) = 0.d0

             END IF


          ENDDO

       END DO

    END IF

  END SUBROUTINE eval_flux_kt

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


  SUBROUTINE average_kt( a1 , a2 , w1 , w2 , w_avg )

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: a1(:) , a2(:)
    REAL*8, INTENT(IN) :: w1(:) , w2(:)
    REAL*8, INTENT(OUT) :: w_avg(:)

    INTEGER :: n
    INTEGER :: i 

    n = SIZE( a1 )

    DO i=1,n

       IF ( a1(i) .EQ. a2(i) ) THEN

          w_avg(i) = 0.5d0 * ( w1(i) + w2(i) )
          w_avg(i) = 0.d0

       ELSE

          w_avg(i) = ( a2(i) * w1(i) - a1(i) * w2(i) ) / ( a2(i) - a1(i) )  

       END IF

    END DO

  END SUBROUTINE average_kt

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

  SUBROUTINE eval_flux_gforce

    ! to be implemented
    WRITE(*,*) 'method not yet implemented in 2-d case'

  END SUBROUTINE eval_flux_gforce

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

  SUBROUTINE eval_flux_lxf

    ! to be implemented
    WRITE(*,*) 'method not yet implemented in 2-d case'

  END SUBROUTINE eval_flux_lxf


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

  SUBROUTINE reconstruction

    ! External procedures
    USE constitutive_2d, ONLY : qc_to_qp , qp_to_qc
    USE constitutive_2d, ONLY : qrec_to_qc
    USE parameters_2d, ONLY : limiter

    ! External variables
    USE geometry_2d, ONLY : x_comp , x_stag , y_comp , y_stag , dx , dx2 , dy , &
         dy2

    USE parameters_2d, ONLY : reconstr_variables , reconstr_coeff

    IMPLICIT NONE

    REAL*8 :: qc(n_vars)
    REAL*8 :: qrec(n_vars,comp_cells_x,comp_cells_y) 
    REAL*8 :: qrecW(n_vars)
    REAL*8 :: qrecE(n_vars)
    REAL*8 :: qrecS(n_vars)
    REAL*8 :: qrecN(n_vars)
    REAL*8 :: qrec_bdry(n_vars)

    REAL*8 :: qrec_stencil(3)
    REAL*8 :: x_stencil(3)
    REAL*8 :: y_stencil(3)
    REAL*8 :: qrec_prime_x
    REAL*8 :: qrec_prime_y

    INTEGER :: j,k
    INTEGER :: i

    ! Compute the variable to reconstruct (phys or cons)
    DO k = 1,comp_cells_y
    
       DO j = 1,comp_cells_x

          qc = q(1:n_vars,j,k)

          IF ( reconstr_variables .EQ. 'phys' ) THEN

             CALL qc_to_qp( qc , b_cent(j,k) , qrec(1:n_vars,j,k) )

          ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN

             qrec(1:n_vars,j,k) = qc

          END IF
          
       END DO

    ENDDO

    ! Linear reconstruction

    y_loop:DO k = 1,comp_cells_y

       x_loop:DO j = 1,comp_cells_x

          vars_loop:DO i=1,n_vars

             ! x direction
             IF ( comp_cells_x .GT. 1 ) THEN

                ! west boundary
                IF (j.EQ.1) THEN

                   IF ( bcw(i)%flag .EQ. 0 ) THEN

                      x_stencil(1) = x_stag(1)
                      x_stencil(2:3) = x_comp(1:2)

                      qrec_stencil(1) = bcw(i)%value
                      qrec_stencil(2:3) = qrec(i,1:2,k)

                      CALL limit( qrec_stencil , x_stencil , limiter(i) ,       &
                           qrec_prime_x ) 

                   ELSEIF ( bcw(i)%flag .EQ. 1 ) THEN

                      qrec_prime_x = bcw(i)%value

                   ELSEIF ( bcw(i)%flag .EQ. 2 ) THEN

                      qrec_prime_x = ( qrec(i,2,k) - qrec(i,1,k) ) / dx

                   END IF

                   !east boundary
                ELSEIF (j.EQ.comp_cells_x) THEN

                   IF ( bce(i)%flag .EQ. 0 ) THEN

                      qrec_stencil(3) = bce(i)%value
                      qrec_stencil(1:2) = qrec(i,comp_cells_x-1:comp_cells_x,k)

                      x_stencil(3) = x_stag(comp_interfaces_x)
                      x_stencil(1:2) = x_comp(comp_cells_x-1:comp_cells_x)

                      CALL limit( qrec_stencil , x_stencil , limiter(i) ,       &
                           qrec_prime_x ) 

                   ELSEIF ( bce(i)%flag .EQ. 1 ) THEN

                      qrec_prime_x = bce(i)%value

                   ELSEIF ( bce(i)%flag .EQ. 2 ) THEN

                      qrec_prime_x = ( qrec(i,comp_cells_x,k) -                 &
                           qrec(i,comp_cells_x-1,k) ) / dx

                   END IF

                   ! internal x cells
                ELSE

                   x_stencil(1:3) = x_comp(j-1:j+1)
                   qrec_stencil = qrec(i,j-1:j+1,k)

                   CALL limit( qrec_stencil , x_stencil , limiter(i) ,          &
                        qrec_prime_x )

                ENDIF

                qrecw(i) = qrec(i,j,k) - reconstr_coeff * dx2 * qrec_prime_x
                qrece(i) = qrec(i,j,k) + reconstr_coeff * dx2 * qrec_prime_x

                ! positivity preserving reconstruction for h (i=1, 1st variable)
                IF ( i .EQ. 1 ) THEN

                   IF ( qrece(i) .LT. b_stag_x(j+1,k) ) THEN

                      qrec_prime_x = ( b_stag_x(j+1,k) - qrec(i,j,k) ) / dx2

                      qrece(i) = qrec(i,j,k) + dx2 * qrec_prime_x

                      qrecw(i) = 2.d0 * qrec(i,j,k) - qrece(i) 

                   ENDIF

                   IF ( qrecw(i) .LT. b_stag_x(j,k) ) THEN

                      qrec_prime_x = ( qrec(i,j,k) - b_stag_x(j,k) ) / dx2

                      qrecw(i) = qrec(i,j,k) - dx2 * qrec_prime_x

                      qrece(i) = 2.d0 * qrec(i,j,k) - qrecw(i) 

                   ENDIF

                END IF

             END IF

             ! y-direction
             check_comp_cells_y:IF ( comp_cells_y .GT. 1 ) THEN

                ! South boundary
                check_y_boundary:IF (k.EQ.1) THEN

                   IF ( bcs(i)%flag .EQ. 0 ) THEN

                      qrec_stencil(1) = bcs(i)%value
                      qrec_stencil(2:3) = qrec(i,j,1:2)

                      y_stencil(1) = y_stag(1)
                      y_stencil(2:3) = y_comp(1:2)

                      CALL limit( qrec_stencil , y_stencil , limiter(i) ,       &
                           qrec_prime_y ) 

                   ELSEIF ( bcs(i)%flag .EQ. 1 ) THEN

                      qrec_prime_y = bcs(i)%value

                   ELSEIF ( bcs(i)%flag .EQ. 2 ) THEN

                      qrec_prime_y = ( qrec(i,j,2) - qrec(i,j,1) ) / dy 

                   END IF

                   ! North boundary
                ELSEIF ( k .EQ. comp_cells_y ) THEN

                   IF ( bcn(i)%flag .EQ. 0 ) THEN

                      qrec_stencil(3) = bcn(i)%value
                      qrec_stencil(1:2) = qrec(i,j,comp_cells_y-1:comp_cells_y)

                      y_stencil(3) = y_stag(comp_interfaces_y)
                      y_stencil(1:2) = y_comp(comp_cells_y-1:comp_cells_y)

                      CALL limit( qrec_stencil , y_stencil , limiter(i) ,       &
                           qrec_prime_y ) 

                   ELSEIF ( bcn(i)%flag .EQ. 1 ) THEN

                      qrec_prime_y = bcn(i)%value

                   ELSEIF ( bcn(i)%flag .EQ. 2 ) THEN

                      qrec_prime_y = ( qrec(i,j,comp_cells_y) -                 &
                           qrec(i,j,comp_cells_y-1) ) / dy 

                   END IF

                   ! Internal y cells
                ELSE

                   y_stencil(1:3) = y_comp(k-1:k+1)
                   qrec_stencil = qrec(i,j,k-1:k+1)

                   CALL limit( qrec_stencil , y_stencil , limiter(i) ,          &
                        qrec_prime_y )

                ENDIF check_y_boundary

                qrecs(i) = qrec(i,j,k) - reconstr_coeff * dy2 * qrec_prime_y
                qrecn(i) = qrec(i,j,k) + reconstr_coeff * dy2 * qrec_prime_y

                ! positivity preserving reconstruction for h
                IF ( i .EQ. 1 ) THEN

                   IF ( qrecn(i) .LT. b_stag_y(j,k+1) ) THEN

                      qrec_prime_y = ( b_stag_y(j,k+1) - qrec(i,j,k) ) / dy2

                      qrecn(i) = qrec(i,j,k) + dy2 * qrec_prime_y

                      qrecs(i) = 2.d0 * qrec(i,j,k) - qrecn(i) 

                   ENDIF

                   IF ( qrecs(i) .LT. b_stag_y(j,k) ) THEN

                      qrec_prime_y = ( qrec(i,j,k) - b_stag_y(j,k) ) / dy2

                      qrecs(i) = qrec(i,j,k) - dy2 * qrec_prime_y

                      qrecn(i) = 2.d0 * qrec(i,j,k) - qrecs(i)

                   ENDIF

                END IF

             ENDIF check_comp_cells_y

          ENDDO vars_loop


          IF ( reconstr_variables .EQ. 'phys' ) THEN

             ! Convert back from physical to conservative variables
             IF ( comp_cells_x .GT. 1 ) THEN

                CALL qp_to_qc( qrecw , b_stag_x(j,k) , q_interfacer(:,j,k) )
                CALL qp_to_qc( qrece , b_stag_x(j+1,k) , q_interfacel(:,j+1,k) )

             END IF

             IF ( comp_cells_y .GT. 1 ) THEN

                CALL qp_to_qc( qrecs , b_stag_y(j,k) , q_interfacet(:,j,k) )
                CALL qp_to_qc( qrecn , b_stag_y(j,k+1) , q_interfaceb(:,j,k+1) )

             END IF

          ELSE IF ( reconstr_variables .EQ. 'cons' ) THEN

             IF ( comp_cells_x .GT. 1 ) THEN

                CALL qrec_to_qc( qrecw, b_stag_x(j,k) , q_interfacer(:,j,k) )
                CALL qrec_to_qc( qrece, b_stag_x(j+1,k) , q_interfacel(:,j+1,k) )

             END IF

             IF ( comp_cells_y .GT. 1 ) THEN

                CALL qrec_to_qc( qrecs, b_stag_y(j,k) , q_interfacet(:,j,k) )
                CALL qrec_to_qc( qrecn, b_stag_y(j,k+1) , q_interfaceb(:,j,k+1) )

             END IF

          END IF

          IF ( j.EQ.0 ) THEN

             WRITE(*,*) 'qrecW',qrecw
             WRITE(*,*) 'qrecE',qrece

          END IF

          ! boundary conditions

          IF ( comp_cells_y .GT. 1 ) THEN

             ! South q_interfaceB(:,j,1)
             IF ( k .EQ. 1 ) THEN

                DO i=1,n_vars

                   IF ( bcs(i)%flag .EQ. 0 ) THEN

                      qrec_bdry(i) = bcs(i)%value 

                   ELSE

                      qrec_bdry(i) = qrecs(i)

                   END IF

                ENDDO

                IF ( reconstr_variables .EQ. 'phys' ) THEN

                   CALL qp_to_qc( qrec_bdry, b_stag_y(j,1), q_interfaceb(:,j,1) )

                ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN

                   CALL qrec_to_qc( qrec_bdry , b_stag_y(j,1) ,                 &
                        q_interfaceb(:,j,1) )

                END IF

             ENDIF

             ! North qT(i,j,comp_interfaces_y)
             IF(k.EQ.comp_cells_y)THEN

                DO i=1,n_vars

                   IF ( bcn(i)%flag .EQ. 0 ) THEN

                      qrec_bdry(i) = bcn(i)%value 

                   ELSE

                      qrec_bdry(i) = qrecn(i)

                   END IF

                ENDDO

                IF ( reconstr_variables .EQ. 'phys' ) THEN

                   CALL qp_to_qc( qrec_bdry , b_stag_y(j,comp_interfaces_y) ,   &
                        q_interfacet(:,j,comp_interfaces_y) )

                ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN

                   CALL qrec_to_qc( qrec_bdry , b_stag_y(j,comp_interfaces_y) , &
                        q_interfacet(:,j,comp_interfaces_y) )

                END IF

             ENDIF

          END IF

          IF ( comp_cells_x .GT. 1 ) THEN

             ! West q_interfaceL(:,1,k)
             IF ( j.EQ.1 ) THEN

                DO i=1,n_vars

                   IF ( bcw(i)%flag .EQ. 0 ) THEN

                      qrec_bdry(i) = bcw(i)%value 

                   ELSE

                      qrec_bdry(i) = qrecw(i)

                   END IF

                ENDDO

                IF ( reconstr_variables .EQ. 'phys' ) THEN

                   CALL qp_to_qc( qrec_bdry, b_stag_x(1,k), q_interfacel(:,1,k) )

                ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN

                   CALL qrec_to_qc( qrec_bdry , b_stag_x(1,k) ,                 &
                        q_interfacel(:,1,k) )

                END IF

                q_interfacer(:,1,k) = q_interfacel(:,1,k)

             ENDIF

             ! East q_interfaceR(:,comp_interfaces_x,k)
             IF ( j.EQ.comp_cells_x ) THEN

                DO i=1,n_vars

                   IF ( bce(i)%flag .EQ. 0 ) THEN

                      qrec_bdry(i) = bce(i)%value 

                   ELSE

                      qrec_bdry(i) = qrece(i)

                   END IF

                ENDDO

                IF ( reconstr_variables .EQ. 'phys' ) THEN

                   CALL qp_to_qc( qrec_bdry , b_stag_x(comp_interfaces_x,k) ,   &
                        q_interfacer(:,comp_interfaces_x,k) )

                ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN

                   CALL qrec_to_qc( qrec_bdry , b_stag_x(comp_interfaces_x,k) , &
                        q_interfacer(:,comp_interfaces_x,k) )

                END IF

                q_interfacel(:,comp_interfaces_x,k) =                           &
                     q_interfacer(:,comp_interfaces_x,k)

             ENDIF

          END IF

       END DO x_loop

    END DO y_loop

  END SUBROUTINE reconstruction


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

  SUBROUTINE eval_speeds

    ! External procedures
    USE constitutive_2d, ONLY : eval_local_speeds_x, eval_local_speeds_y 
    USE constitutive_2d, ONLY : eval_local_speeds2_x, eval_local_speeds2_y 

    IMPLICIT NONE

    REAL*8 :: abslambdaL_min(n_vars) , abslambdaL_max(n_vars)
    REAL*8 :: abslambdaR_min(n_vars) , abslambdaR_max(n_vars)
    REAL*8 :: abslambdaB_min(n_vars) , abslambdaB_max(n_vars)
    REAL*8 :: abslambdaT_min(n_vars) , abslambdaT_max(n_vars)
    REAL*8 :: min_r(n_vars) , max_r(n_vars)

    INTEGER :: j,k

    IF ( comp_cells_x .GT. 1 ) THEN
    
       DO j = 1,comp_interfaces_x
          
          DO k = 1, comp_cells_y
             
             CALL eval_local_speeds2_x( q_interfacel(:,j,k) , b_stag_x(j,k) ,   &
                  abslambdal_min , abslambdal_max )
             
             CALL eval_local_speeds2_x( q_interfacer(:,j,k) , b_stag_x(j,k) ,   &
                  abslambdar_min , abslambdar_max )
             
             min_r = min(abslambdal_min , abslambdar_min , 0.0d0)
             max_r = max(abslambdal_max , abslambdar_max , 0.0d0)
             
             a_interface_xneg(:,j,k) = min_r
             a_interface_xpos(:,j,k) = max_r
             
          ENDDO
          
       END DO

    END IF

    IF ( comp_cells_y .GT. 1 ) THEN
    
       DO j = 1,comp_cells_x
          
          DO k = 1,comp_interfaces_y
             
             CALL eval_local_speeds2_y( q_interfaceb(:,j,k) , b_stag_y(j,k) ,   &
                  abslambdab_min , abslambdab_max )
             
             CALL eval_local_speeds2_y( q_interfacet(:,j,k) , b_stag_y(j,k) ,   &
                  abslambdat_min , abslambdat_max )
             
             min_r = min(abslambdab_min , abslambdat_min , 0.0d0)
             max_r = max(abslambdab_max , abslambdat_max , 0.0d0)
             
             a_interface_yneg(:,j,k) = min_r
             a_interface_ypos(:,j,k) = max_r
             
          ENDDO
          
       END DO

    END IF
       
  END SUBROUTINE eval_speeds


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

  SUBROUTINE limit( v , z , limiter , slope_lim )

    USE parameters_2d, ONLY : theta

    IMPLICIT none

    REAL*8, INTENT(IN) :: v(3)
    REAL*8, INTENT(IN) :: z(3)
    INTEGER, INTENT(IN) :: limiter

    REAL*8, INTENT(OUT) :: slope_lim

    REAL*8 :: a , b , c

    REAL*8 :: sigma1 , sigma2

    a = ( v(3) - v(2) ) / ( z(3) - z(2) )
    b = ( v(2) - v(1) ) / ( z(2) - z(1) )
    c = ( v(3) - v(1) ) / ( z(3) - z(1) )

    SELECT CASE (limiter)

    CASE ( 0 )

       slope_lim = 0.d0

    CASE ( 1 )

       ! minmod

       slope_lim = minmod(a,b)

    CASE ( 2 )

       ! superbee

       sigma1 = minmod( a , 2.d0*b )
       sigma2 = minmod( 2.d0*a , b )
       slope_lim = maxmod( sigma1 , sigma2 )

    CASE ( 3 )

       ! generalized minmod

       slope_lim = minmod( c , theta * minmod( a , b ) )

    CASE ( 4 )

       ! monotonized central-difference (MC, LeVeque p.112)

       slope_lim = minmod( c , 2.0 * minmod( a , b ) )

    END SELECT

  END SUBROUTINE limit


  REAL*8 FUNCTION minmod(a,b)

    IMPLICIT none

    REAL*8 :: a , b , sa , sb 

    IF ( a*b .EQ. 0.d0 ) THEN

       minmod = 0.d0

    ELSE

       sa = a / abs(a)
       sb = b / abs(b)

       minmod = 0.5d0 * ( sa+sb ) * min( abs(a) , abs(b) )

    END IF

  END FUNCTION minmod

  REAL*8 function maxmod(a,b)

    IMPLICIT none

    REAL*8 :: a , b , sa , sb 

    IF ( a*b .EQ. 0.d0 ) THEN

       maxmod = 0.d0

    ELSE

       sa = a / abs(a)
       sb = b / abs(b)

       maxmod = 0.5d0 * ( sa+sb ) * max( abs(a) , abs(b) )

    END IF

  END function maxmod

END MODULE solver_2d