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,comp_cells_xy
  USE geometry_2d, ONLY : comp_interfaces_x,comp_interfaces_y

  USE geometry_2d, ONLY : source_cell

  USE variables, ONLY : wp , sp

  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

  ! external procedures
  USE geometry_2d, ONLY : limit

  USE omp_lib

  IMPLICIT none

  REAL(wp) :: t

  REAL(wp), ALLOCATABLE :: q(:,:,:)        
  REAL(wp), ALLOCATABLE :: q0(:,:,:)        
  REAL(wp), ALLOCATABLE :: q_fv(:,:,:)     

  REAL(wp), ALLOCATABLE :: q_interfacel(:,:,:)        
  REAL(wp), ALLOCATABLE :: q_interfacer(:,:,:)
  REAL(wp), ALLOCATABLE :: q_interfaceb(:,:,:)        
  REAL(wp), ALLOCATABLE :: q_interfacet(:,:,:)

  REAL(wp), ALLOCATABLE :: qp_interfacel(:,:,:)        
  REAL(wp), ALLOCATABLE :: qp_interfacer(:,:,:)
  REAL(wp), ALLOCATABLE :: qp_interfaceb(:,:,:)        
  REAL(wp), ALLOCATABLE :: qp_interfacet(:,:,:)

  REAL(wp), ALLOCATABLE :: q_cellnw(:,:,:)        
  REAL(wp), ALLOCATABLE :: q_cellne(:,:,:)
  REAL(wp), ALLOCATABLE :: q_cellsw(:,:,:)        
  REAL(wp), ALLOCATABLE :: q_cellse(:,:,:)

  REAL(wp), ALLOCATABLE :: qp_cellnw(:,:,:)        
  REAL(wp), ALLOCATABLE :: qp_cellne(:,:,:)
  REAL(wp), ALLOCATABLE :: qp_cellsw(:,:,:)        
  REAL(wp), ALLOCATABLE :: qp_cellse(:,:,:)


  REAL(wp), ALLOCATABLE :: q1max(:,:)

  REAL(wp), ALLOCATABLE :: a_interface_x_max(:,:,:)
  REAL(wp), ALLOCATABLE :: a_interface_y_max(:,:,:)


  REAL(wp), ALLOCATABLE :: a_interface_xneg(:,:,:)
  REAL(wp), ALLOCATABLE :: a_interface_xpos(:,:,:)
  REAL(wp), ALLOCATABLE :: a_interface_yneg(:,:,:)
  REAL(wp), ALLOCATABLE :: a_interface_ypos(:,:,:)
  REAL(wp), ALLOCATABLE :: h_interface_x(:,:,:)
  REAL(wp), ALLOCATABLE :: h_interface_y(:,:,:)
  REAL(wp), ALLOCATABLE :: qp(:,:,:)

  REAL(wp), ALLOCATABLE :: source_xy(:,:)

  LOGICAL, ALLOCATABLE :: solve_mask(:,:)
  LOGICAL, ALLOCATABLE :: solve_mask_x(:,:)
  LOGICAL, ALLOCATABLE :: solve_mask_y(:,:)

  INTEGER :: solve_cells
  INTEGER :: solve_interfaces_x
  INTEGER :: solve_interfaces_y

  REAL(wp) :: dt

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

  INTEGER :: i_rk

  REAL(wp), ALLOCATABLE :: a_tilde_ij(:,:)
  REAL(wp), ALLOCATABLE :: a_dirk_ij(:,:)

  REAL(wp), ALLOCATABLE :: omega_tilde(:)

  REAL(wp), ALLOCATABLE :: omega(:)

  REAL(wp), ALLOCATABLE :: a_tilde(:)

  REAL(wp), ALLOCATABLE :: a_dirk(:)

  REAL(wp) :: a_diag

  REAL(wp), ALLOCATABLE :: q_rk(:,:,:,:)

  REAL(wp), ALLOCATABLE :: qp_rk(:,:,:,:)

  REAL(wp), ALLOCATABLE :: divflux(:,:,:,:)

  REAL(wp), ALLOCATABLE :: nh(:,:,:,:)

  LOGICAL :: normalize_q

  LOGICAL :: normalize_f

  LOGICAL :: opt_search_nl

  REAL(wp), ALLOCATABLE :: residual_term(:,:,:)

  INTEGER, ALLOCATABLE :: j_cent(:)
  INTEGER, ALLOCATABLE :: k_cent(:)

  INTEGER, ALLOCATABLE :: j_stag_x(:)
  INTEGER, ALLOCATABLE :: k_stag_x(:)

  INTEGER, ALLOCATABLE :: j_stag_y(:)
  INTEGER, ALLOCATABLE :: k_stag_y(:)

  REAL(wp), ALLOCATABLE :: q_mg_old(:,:,:)
  REAL(wp), ALLOCATABLE :: q_mg_new(:,:,:)
  
  
CONTAINS

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

  SUBROUTINE allocate_solver_variables

    IMPLICIT NONE

    REAL(wp) :: 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( qp( n_vars+2 , comp_cells_x , comp_cells_y ) )

    q(1:n_vars,1:comp_cells_x,1:comp_cells_y) = 0.0_wp
    qp(1:n_vars+2,1:comp_cells_x,1:comp_cells_y) = 0.0_wp

    ALLOCATE( q1max( 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 ) )

    a_interface_xneg = 0.0_wp
    a_interface_xpos = 0.0_wp
    
    ALLOCATE( h_interface_x( n_eqns , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( h_interface_y( n_eqns , comp_cells_x, comp_interfaces_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 ) )

    a_interface_yneg = 0.0_wp
    a_interface_ypos = 0.0_wp

    
    ALLOCATE( qp_interfacel( n_vars+2 , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( qp_interfacer( n_vars+2 , comp_interfaces_x, comp_cells_y ) )
    ALLOCATE( qp_interfaceb( n_vars+2 , comp_cells_x, comp_interfaces_y ) )
    ALLOCATE( qp_interfacet( n_vars+2 , comp_cells_x, comp_interfaces_y ) )


    ALLOCATE( q_cellnw( n_vars , comp_cells_x , comp_cells_y ) )
    ALLOCATE( q_cellne( n_vars , comp_cells_x , comp_cells_y ) )
    ALLOCATE( q_cellsw( n_vars , comp_cells_x , comp_cells_y ) )
    ALLOCATE( q_cellse( n_vars , comp_cells_x , comp_cells_y ) )

    ALLOCATE( qp_cellnw( n_vars+2 , comp_cells_x , comp_cells_y ) )
    ALLOCATE( qp_cellne( n_vars+2 , comp_cells_x , comp_cells_y ) )
    ALLOCATE( qp_cellsw( n_vars+2 , comp_cells_x , comp_cells_y ) )
    ALLOCATE( qp_cellse( n_vars+2 , comp_cells_x , comp_cells_y ) )


    ALLOCATE ( a_interface_x_max(n_eqns,comp_interfaces_x,comp_cells_y) )
    ALLOCATE ( a_interface_y_max(n_eqns,comp_cells_x,comp_interfaces_y) )

    ALLOCATE( solve_mask( comp_cells_x , comp_cells_y ) )

    solve_mask(1,1:comp_cells_y) = .true.
    solve_mask(comp_cells_x,1:comp_cells_y) = .true.
    solve_mask(1:comp_cells_x,1) = .true.
    solve_mask(1:comp_cells_x,comp_cells_y) = .true.


    ALLOCATE( solve_mask_x( comp_interfaces_x , comp_cells_y ) )
    ALLOCATE( solve_mask_y( comp_cells_x , comp_interfaces_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.0_wp

    ! Weight coefficients of the explicit part in the final assemblage
    omega_tilde = 0.0_wp

    ! Implicit part coefficients (a_dirk_ij=0 for j>i)
    a_dirk_ij = 0.0_wp

    ! Weight coefficients of the explicit part in the final assemblage
    omega = 0.0_wp

    gamma = 1.0_wp - 1.0_wp / sqrt(2.0_wp)
    delta = 1.0_wp - 1.0_wp / ( 2.0_wp * gamma )

    IF ( n_rk .EQ. 1 ) THEN

       a_tilde_ij(1,1) = 1.0_wp

       omega_tilde(1) = 1.0_wp

       a_dirk_ij(1,1) = 0.0_wp

       omega(1) = 0.0_wp

    ELSEIF ( n_rk .EQ. 2 ) THEN

       a_tilde_ij(2,1) = 1.0_wp

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

       a_dirk_ij(2,2) = 1.0_wp

       omega(1) = 0.0_wp
       omega(2) = 1.0_wp

    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.5_wp
       a_tilde_ij(3,1) = 0.5_wp
       a_tilde_ij(3,2) = 0.5_wp

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

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

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

    ELSEIF ( n_rk .EQ. 4 ) THEN

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

       a_tilde_ij(2,1) = 0.5_wp
       a_tilde_ij(3,1) = 1.0_wp / 3.0_wp
       a_tilde_ij(4,2) = 1.0_wp

       omega_tilde(1) = 0.0_wp
       omega_tilde(2) = 1.0_wp
       omega_tilde(3) = 0.0_wp
       omega_tilde(4) = 0.0_wp

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

       omega(1) = 0.0_wp
       omega(2) = 0.0_wp
       omega(3) = 0.75_wp
       omega(4) = 0.25_wp

    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( qp_rk( n_vars+2 , 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( residual_term( n_vars , comp_cells_x , comp_cells_y ) )

    comp_cells_xy = comp_cells_x * comp_cells_y

    ALLOCATE( j_cent( comp_cells_xy ) )
    ALLOCATE( k_cent( comp_cells_xy ) )

    ALLOCATE( j_stag_x( comp_interfaces_x * comp_cells_y ) )
    ALLOCATE( k_stag_x( comp_interfaces_x * comp_cells_y ) )

    ALLOCATE( j_stag_y( comp_cells_x * comp_interfaces_y ) )
    ALLOCATE( k_stag_y( comp_cells_x * comp_interfaces_y ) )

    RETURN
    
  END SUBROUTINE allocate_solver_variables

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

  SUBROUTINE deallocate_solver_variables

    DEALLOCATE( q , q0 )

    DEALLOCATE( qp )

    DEALLOCATE( q1max )

    DEALLOCATE( q_fv )

    DEALLOCATE( q_interfacel )
    DEALLOCATE( q_interfacer )
    DEALLOCATE( a_interface_xneg )
    DEALLOCATE( a_interface_xpos )

    DEALLOCATE( h_interface_x )
    DEALLOCATE( h_interface_y )


    DEALLOCATE( q_interfaceb )
    DEALLOCATE( q_interfacet )
    DEALLOCATE( a_interface_yneg )
    DEALLOCATE( a_interface_ypos )
    
    DEALLOCATE( qp_interfacel )
    DEALLOCATE( qp_interfacer )
    DEALLOCATE( qp_interfaceb )
    DEALLOCATE( qp_interfacet )


    DEALLOCATE( q_cellnw )
    DEALLOCATE( q_cellne )
    DEALLOCATE( q_cellsw )
    DEALLOCATE( q_cellse )

    DEALLOCATE( qp_cellnw )
    DEALLOCATE( qp_cellne )
    DEALLOCATE( qp_cellsw )
    DEALLOCATE( qp_cellse )


    DEALLOCATE( a_interface_x_max )
    DEALLOCATE( a_interface_y_max )

    DEALLOCATE( solve_mask )

    DEALLOCATE( solve_mask_x )
    DEALLOCATE( solve_mask_y )

    DEALLOCATE( source_xy )


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


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


    DEALLOCATE( a_tilde )
    DEALLOCATE( a_dirk )

    DEALLOCATE( q_rk )
    DEALLOCATE( qp_rk )
    DEALLOCATE( divflux )
    DEALLOCATE( nh )
    DEALLOCATE( residual_term )

    DEALLOCATE( j_cent )
    DEALLOCATE( k_cent )

    DEALLOCATE( j_stag_x )
    DEALLOCATE( k_stag_x )

    DEALLOCATE( j_stag_y )
    DEALLOCATE( k_stag_y )


    RETURN
    
  END SUBROUTINE deallocate_solver_variables


  SUBROUTINE allocate_multigrid

    ALLOCATE( q_mg_old(n_vars , comp_cells_x , comp_cells_y) )
    ALLOCATE( q_mg_new(n_vars , 2*comp_cells_x , 2*comp_cells_y) )
  
  END SUBROUTINE allocate_multigrid

  SUBROUTINE deallocate_multigrid

    DEALLOCATE( q_mg_old )
    DEALLOCATE( q_mg_new )
  
  END SUBROUTINE deallocate_multigrid

  SUBROUTINE remap_solution

  USE geometry_2d, ONLY : x0 , y0 , cell_size
  USE geometry_2d, ONLY : x_stag , y_stag
  
  USE geometry_2d, ONLY : regrid_scalar
    
    IMPLICIT NONE

    INTEGER :: i ,j ,k
    
    REAL(wp) :: xl , xr
    REAL(wp) :: yl , yr

    WRITE(*,*) 'INT(q(1,:,:)=',sum(q(1,:,:))*cell_size*cell_size
        
    DO k = 1,2*comp_cells_y

       yl = y0 + (k-1)*0.5_wp*cell_size
       yr = y0 + (k)*0.5_wp*cell_size
       
       DO j = 1,2*comp_cells_x

          xl = x0 + (j-1)*0.5_wp*cell_size
          xr = x0 + (j)*0.5_wp*cell_size
          
          DO i=1,n_vars
             
             CALL regrid_scalar(x_stag, y_stag, q(i,:,:), xl, xr , yl, yr, q_mg_new(i,j,k) )

          END DO

       END DO

    END DO

    WRITE(*,*) 'size q_mg_new',SIZE(q_mg_new)
    WRITE(*,*) 'SUM(q_mg_new(1,:,:)=',sum(q_mg_new(1,:,:))*0.25_wp*cell_size*cell_size
    READ(*,*)
    
    RETURN
    
  END SUBROUTINE remap_solution
  
  !******************************************************************************
  !
  !
  !
  !******************************************************************************

  SUBROUTINE check_solve

    IMPLICIT NONE

    INTEGER :: i,j,k

    !$OMP PARALLEL
    
    !$OMP WORKSHARE
    WHERE ( q(1,2:comp_cells_x-1,2:comp_cells_y-1) .GT. 0.0_wp )  &
         solve_mask(2:comp_cells_x,2:comp_cells_y) = .true.
    
    !$OMP END WORKSHARE
        
    !$OMP SINGLE

    DO i = 1,n_rk

       ! solution domain is extended to neighbours of positive-mass cells
       solve_mask(2:comp_cells_x-1,2:comp_cells_y-1) =                          &
            solve_mask(2:comp_cells_x-1,2:comp_cells_y-1) .OR.                  &
            solve_mask(1:comp_cells_x-2,2:comp_cells_y-1) .OR.                  &
            solve_mask(3:comp_cells_x,2:comp_cells_y-1) .OR.                    &
            solve_mask(2:comp_cells_x-1,1:comp_cells_y-2) .OR.                  &
            solve_mask(2:comp_cells_x-1,3:comp_cells_y) 

    END DO

    !$OMP END SINGLE

    !$OMP WORKSHARE

    solve_mask_x(1:comp_interfaces_x,1:comp_cells_y) = .false.
    solve_mask_y(1:comp_cells_x,1:comp_interfaces_y) = .false.

    !$OMP END WORKSHARE

    !$OMP END PARALLEL

    !----- check for cells where computation is needed
    i = 0
        
    DO k = 1,comp_cells_y

       DO j = 1,comp_cells_x

          IF ( solve_mask(j,k) ) THEN

             i = i+1
             j_cent(i) = j
             k_cent(i) = k

             solve_mask_x(j,k) = .true.
             solve_mask_x(j+1,k) = .true.
             solve_mask_y(j,k) = .true.
             solve_mask_y(j,k+1) = .true.

          END IF

       END DO

    END DO

    solve_cells = i

    !----- check for y-interfaces where computation is needed
    i = 0
        
    DO k = 1,comp_cells_y

       DO j = 1,comp_interfaces_x

          IF ( solve_mask_x(j,k) ) THEN

             i = i+1
             j_stag_x(i) = j
             k_stag_x(i) = k

          END IF

       END DO

    END DO

    solve_interfaces_x = i
   
    !----- check for y-interfaces where computation is needed

    i = 0
        
    DO k = 1,comp_interfaces_y

       DO j = 1,comp_cells_x

          IF ( solve_mask_y(j,k) ) THEN

             i = i+1
             j_stag_y(i) = j
             k_stag_y(i) = k

          END IF

       END DO

    END DO

    solve_interfaces_y = i

    RETURN

    WRITE(*,*) solve_mask_y
    READ(*,*)


  END SUBROUTINE check_solve

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

  SUBROUTINE timestep

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

    USE constitutive_2d, ONLY : qc_to_qp

    IMPLICIT none

    REAL(wp) :: dt_cfl

    REAL(wp) :: dt_interface_x, dt_interface_y

    INTEGER :: i,j,k,l

    REAL(wp) :: max_a

    dt = max_dt

    IF ( cfl .NE. -1.0_wp ) THEN

       !$OMP PARALLEL DO private(j,k)

       DO l = 1,solve_cells

          j = j_cent(l)
          k = k_cent(l)

          IF ( q(1,j,k) .GT. 0.0_wp ) THEN

             CALL qc_to_qp( q(1:n_vars,j,k) , qp(1:n_vars+2,j,k) )

          ELSE

             qp(1:n_vars+2,j,k) = 0.0_wp

          END IF

       END DO

       !$OMP END PARALLEL DO

       ! Compute the physical and conservative variables at the interfaces
       CALL reconstruction( q , qp )

       ! Compute the max/min eigenvalues at the interfaces
       CALL eval_speeds


       !$OMP PARALLEL DO private(j,k)
       DO l = 1,solve_interfaces_x

          j = j_stag_x(l)
          k = k_stag_x(l)

          DO i=1,n_vars

             a_interface_x_max(i,j,k) =                                         &
                  max( a_interface_xpos(i,j,k) , -a_interface_xneg(i,j,k) )
 
          END DO

       END DO
       !$OMP END PARALLEL DO
    

       !$OMP PARALLEL DO private(j,k)
       DO l = 1,solve_interfaces_y

          j = j_stag_y(l)
          k = k_stag_y(l)

          DO i=1,n_vars

             a_interface_y_max(i,j,k) =                                         &
                  max( a_interface_ypos(i,j,k) , -a_interface_yneg(i,j,k) )

 
          END DO

       END DO
       !$OMP END PARALLEL DO
       
       DO l = 1,solve_cells

          j = j_cent(l)
          k = k_cent(l)

          max_a =  max( maxval(a_interface_x_max(:,j,k)) ,                      &
               maxval(a_interface_x_max(:,j+1,k)) )

          IF ( max_a .GT. 0.0_wp ) THEN

             dt_interface_x = cfl * dx / max_a

          ELSE

             dt_interface_x = dt

          END IF

          max_a =  max( maxval(a_interface_y_max(:,j,k)) ,                      &
               maxval(a_interface_y_max(:,j,k+1)) )

          IF ( max_a .GT. 0.0_wp ) THEN

             dt_interface_y = cfl * dy / max_a

          ELSE

             dt_interface_y = dt

          END IF

          dt_cfl = min( dt_interface_x , dt_interface_y )

          dt = min(dt,dt_cfl)

       END DO

    END IF

    RETURN

  END SUBROUTINE timestep

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

  SUBROUTINE imex_rk_solver

    USE constitutive_2d, ONLY : eval_nonhyperbolic_terms

    USE constitutive_2d, ONLY : qc_to_qp

    USE geometry_2d, ONLY : cell_fract , dx_rel , dy_rel

    IMPLICIT NONE

    REAL(wp) :: q_si(n_vars)
    REAL(wp) :: q_guess(n_vars)
    INTEGER :: j,k,l
    REAL(wp) :: Rj_not_impl(n_eqns)

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

    !$OMP PARALLEL DO private(j,k,q_guess,q_si,Rj_not_impl)
    init_cells_loop:DO l = 1,solve_cells

       j = j_cent(l)
       k = k_cent(l)

       ! Initialization of the solution guess
       q0( 1:n_vars , j , k ) =                          &
            q( 1:n_vars , j , k )
       
       ! Initialization of the variables for the Runge-Kutta scheme
       q_rk( 1:n_vars , j , k , 1:n_rk ) = 0.0_wp
       qp_rk( 1:n_vars+2 , j , k , 1:n_rk ) = 0.0_wp
       divflux(1:n_eqns , j , k , 1:n_rk ) = 0.0_wp
       nh( 1:n_eqns, j , k , 1:n_rk ) = 0.0_wp
       
    END DO init_cells_loop
    !$OMP END PARALLEL DO

    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.0_wp
       a_dirk = 0.0_wp

       ! 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)

       !$OMP PARALLEL 
       !$OMP DO private(j,k,q_guess,q_si,Rj_not_impl)

       solve_cells_loop:DO l = 1,solve_cells

          j = j_cent(l)
          k = k_cent(l)

          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

          ! 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( divflux(1:n_eqns,j,k,1:i_rk) , a_tilde(1:i_rk) ) &
               - matmul( nh(1:n_eqns,j,k,1:i_rk) , a_dirk(1:i_rk) ) )

          IF ( verbose_level .GE. 2 ) THEN

             WRITE(*,*) 'q_guess',q_guess
             IF ( q_guess(1) .GT. 0.0_wp  ) THEN 

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

             END IF

          END IF

          adiag_pos:IF ( a_diag .NE. 0.0_wp ) THEN

             pos_thick:IF ( ( q_fv(1,j,k) .GT.  0.0_wp ) .OR. ( cell_fract(j,k) .GT.  0.0_wp ) )   THEN
             !pos_thick:IF ( q_fv(1,j,k) .GT.  0.0_wp )   THEN

                ! Initialize the guess for the NR solver
                q_guess(1:n_vars) = q_fv(1:n_vars,j,k)

                rj_not_impl =  matmul( divflux(1:n_eqns,j,k,1:i_rk-1) , a_tilde(1:i_rk-1) ) &
                     - matmul( nh(1:n_eqns,j,k,1:i_rk-1) , a_dirk(1:i_rk-1) )
                
                ! Solve the implicit system to find the solution at the 
                ! i_RK step of the IMEX RK procedure
                CALL solve_rk_step( cell_fract(j,k) , dx_rel(j,k) , dy_rel(j,k),&
                     q_guess(1:n_vars) , q0(1:n_vars,j,k ) ,     &
                     a_tilde , a_dirk , a_diag , rj_not_impl ,                  &
                     divflux( 1:n_eqns , j , k , 1:n_rk ) ,                     &
                     nh( 1:n_eqns , j , k , 1:n_rk ) )

                IF ( comp_cells_y .EQ. 1 ) THEN

                   q_guess(3) = 0.0_wp

                END IF

                IF ( comp_cells_x .EQ. 1 ) THEN

                   q_guess(2) = 0.0_wp

                END IF

                ! Eval and store the implicit term at the i_RK step
                CALL eval_nonhyperbolic_terms( cell_fract(j,k), dx_rel(j,k) ,   &
                     dy_rel(j,k) , r_qj = q_guess ,                             &
                     r_nh_term_impl = nh(1:n_eqns,j,k,i_rk) )

             ELSE

                ! If h=0 nothing has to be changed 
                q_guess(1:n_vars) = q_fv( 1:n_vars , j , k ) 
                nh(1:n_eqns,j,k,i_rk) = 0.0_wp

             END IF pos_thick

          END IF adiag_pos

          IF ( a_diag .NE. 0.0_wp ) THEN

             ! Update the implicit term with correction on the new velocity
             nh(1:n_vars,j,k,i_rk) = ( q_guess(1:n_vars) - q_fv(1:n_vars,j,k) ) &
                  / ( 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

             IF ( q_guess(1) .GT. 0.0_wp ) THEN

                CALL qc_to_qp( q_guess , qp(1:n_vars+2,j,k) )
                WRITE(*,*) 'imex_RK_solver: qp',qp(1:n_vars+2,j,k)

             END IF
             
             READ(*,*)

          END IF

          IF ( omega_tilde(i_rk) .GT. 0.0_wp ) THEN
          
             IF ( q_rk(1,j,k,i_rk) .GT. 0.0_wp ) THEN

                CALL qc_to_qp( q_rk(1:n_vars,j,k,i_rk) ,                        &
                     qp_rk(1:n_vars+2,j,k,i_rk) )

             ELSE

                qp_rk(1:n_vars+2,j,k,i_rk) = 0.0_wp

             END IF
             
          END IF

       END DO solve_cells_loop

       !$OMP END DO

       !$OMP END PARALLEL 

       IF ( omega_tilde(i_rk) .GT. 0.0_wp ) THEN

          ! 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) ,              &
               qp_rk(1:n_vars+2,1:comp_cells_x,1:comp_cells_y,i_rk) ,           &
               divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )

       END IF

    END DO runge_kutta

    !$OMP PARALLEL DO private(j,k)

    assemble_sol:DO l = 1,solve_cells

       j = j_cent(l)
       k = k_cent(l)

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

       IF ( verbose_level .GE. 1 ) THEN

          WRITE(*,*) 'cell jk =',j,k
          WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
          IF ( q0(1,j,k) .GT. 0.0_wp ) THEN

             CALL qc_to_qp(q0(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
             WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars+2,j,k)

          END IF

       END IF

       IF ( ( sum(abs( omega_tilde(:)-a_tilde_ij(n_rk,:))) .EQ. 0.0_wp  )       &
            .AND. ( sum(abs(omega(:)-a_dirk_ij(n_rk,:))) .EQ. 0.0_wp ) ) 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

       negative_thickness_check:IF ( q(1,j,k) .LT. 0.0_wp ) THEN

          IF ( q(1,j,k) .GT. -1.e-7_wp ) THEN

             q(1,j,k) = 0.0_wp
             q(2:n_vars,j,k) = 0.0_wp

          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)
             IF ( q0(1,j,k) .GT. 0.0_wp ) THEN

                CALL qc_to_qp(q0(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
                WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars+2,j,k)

             END IF
             WRITE(*,*) 'after imex_RK_solver: qc',q(1:n_vars,j,k)

             WRITE(*,*) 'divFlux(1,j,k,1:n_RK)',divflux(1,j,k,1:n_rk) 

             WRITE(*,*) h_interface_x(1,j+1,k), h_interface_x(1,j,k)
             WRITE(*,*) qp_interfacer(1:n_vars,j,k)
             WRITE(*,*) qp(1:n_vars,j,k)
             WRITE(*,*) qp_interfacel(1:n_vars,j+1,k)

             WRITE(*,*) 'NH(1,j,k,1:n_RK)',nh(1,j,k,1:n_rk) 

             READ(*,*)

          END IF

       END IF negative_thickness_check

    END DO assemble_sol

    !$OMP END PARALLEL DO

    RETURN

  END SUBROUTINE imex_rk_solver

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

  SUBROUTINE solve_rk_step( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj, qj_old, &
       a_tilde, a_dirk, a_diag, rj_not_impl, divfluxj , nhj )

    USE parameters_2d, ONLY : max_nl_iter , tol_rel , tol_abs

    USE constitutive_2d, ONLY : qc_to_qp

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: cell_fract_jk
    REAL(wp), INTENT(IN) :: dx_rel_jk
    REAL(wp), INTENT(IN) :: dy_rel_jk

    REAL(wp), INTENT(INOUT) :: qj(n_vars)
    REAL(wp), INTENT(IN) :: qj_old(n_vars)
    REAL(wp), INTENT(IN) :: a_tilde(n_rk)
    REAL(wp), INTENT(IN) :: a_dirk(n_rk)
    REAL(wp), INTENT(IN) :: a_diag
    REAL(wp), INTENT(IN) :: Rj_not_impl(n_eqns)
    REAL(wp), INTENT(IN) :: divFluxj(n_eqns,n_rk)
    REAL(wp), INTENT(IN) :: NHj(n_eqns,n_rk)

    REAL(wp) :: qj_init(n_vars)

    REAL(wp) :: qj_org(n_vars) , qj_rel(n_vars)

    REAL(wp) :: left_matrix(n_eqns,n_vars)
    REAL(wp) :: right_term(n_eqns)

    REAL(wp) :: scal_f

    REAL(wp) :: coeff_f(n_eqns)

    REAL(wp) :: qj_rel_NR_old(n_vars)
    REAL(wp) :: scal_f_old
    REAL(wp) :: desc_dir(n_vars)
    REAL(wp) :: grad_f(n_vars)

    INTEGER :: pivot(n_vars)

    REAL(wp) :: left_matrix_small22(n_nh,n_nh)
    REAL(wp) :: left_matrix_small21(n_eqns-n_nh,n_nh)
    REAL(wp) :: left_matrix_small11(n_eqns-n_nh,n_vars-n_nh)
    ! REAL(wp) :: left_matrix_small12(n_nh,n_vars-n_nh)

    REAL(wp) :: desc_dir_small2(n_nh)
    INTEGER :: pivot_small2(n_nh)

    REAL(wp) :: desc_dir_small1(n_vars-n_nh)

    INTEGER :: ok

    INTEGER :: i 
    INTEGER :: nl_iter

    REAL(wp), PARAMETER :: STPMX=100.0_wp
    REAL(wp) :: stpmax
    LOGICAL :: check

    REAL(wp), PARAMETER :: TOLF=1.e-10_wp , tolmin=1.e-6_wp
    REAL(wp) :: TOLX

    ! REAL(wp) :: qpj(n_vars+2)

    REAL(wp) :: desc_dir2(n_vars)

    REAL(wp) :: desc_dir_temp(n_vars)

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

    coeff_f(1:n_eqns) = 1.0_wp

    grad_f(1:n_eqns) = 0.0_wp

    qj_init = qj

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

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

       CALL eval_f( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj , qj_old ,       &
            a_diag , coeff_f , rj_not_impl , 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.0_wp ) coeff_f(i) = 1.0_wp/right_term(i)

       END DO

       right_term = coeff_f * right_term

       scal_f = 0.5_wp * 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.e-3_wp )

    ELSE 

       qj_org(1:n_vars) = 1.0_wp

    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( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj , qj_old ,       &
            a_diag , coeff_f , rj_not_impl , 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.e-6_wp ) ) 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( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj_rel ,  &
               qj_org, coeff_f , left_matrix )

          desc_dir_temp = - right_term

          IF ( wp .EQ. sp ) THEN

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

          ELSE

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

          desc_dir = desc_dir_temp

       ELSE

          CALL eval_jacobian( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj_rel ,  &
               qj_org,coeff_f , left_matrix )

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

          ! not needed for computation
          !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 )

          IF ( wp .EQ. sp ) THEN

             CALL sgesv(n_nh,1, left_matrix_small22 , n_nh , pivot_small2 ,     &
                  desc_dir_small2 , n_nh, ok)
             
          ELSE

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

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

       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( cell_fract_jk , dx_rel_jk , dy_rel_jk , 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 , rj_not_impl )

       ELSE

          qj_rel = qj_rel_nr_old + desc_dir

          qj = qj_rel * qj_org

          CALL eval_f( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj , qj_old ,    &
               a_diag , coeff_f , rj_not_impl , 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

       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.0_wp ) /   &
               max( scal_f , 0.5_wp * 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.0_wp ) ) < tolx ) THEN

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

       END IF

    END DO newton_raphson_loop

    RETURN
    
  END SUBROUTINE solve_rk_step

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

  SUBROUTINE lnsrch( cell_fract_jk , dx_rel_jk , dy_rel_jk , 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 , rj_not_impl )

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: cell_fract_jk
    REAL(wp), INTENT(IN) :: dx_rel_jk
    REAL(wp), INTENT(IN) :: dy_rel_jk

    REAL(wp), DIMENSION(:), INTENT(IN) :: qj_rel_NR_old

    REAL(wp), DIMENSION(:), INTENT(IN) :: qj_org

    REAL(wp), DIMENSION(:), INTENT(IN) :: qj_old

    REAL(wp), DIMENSION(:), INTENT(IN) :: grad_f

    REAL(wp), INTENT(IN) :: scal_f_old

    REAL(wp), DIMENSION(:), INTENT(INOUT) :: desc_dir

    REAL(wp), INTENT(IN) :: stpmax

    REAL(wp), DIMENSION(:), INTENT(IN) :: coeff_f

    REAL(wp), DIMENSION(:), INTENT(OUT) :: qj_rel

    REAL(wp), INTENT(OUT) :: scal_f

    REAL(wp), INTENT(OUT) :: right_term(n_eqns)

    LOGICAL, INTENT(OUT) :: check

    REAL(wp), INTENT(IN) :: Rj_not_impl(n_eqns)

    REAL(wp), PARAMETER :: TOLX=epsilon(qj_rel)

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

    REAL(wp) :: scal_f_min , alam_min

    REAL(wp) :: qj(n_vars)

    alf = 1.0e-4_wp

    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.0_wp ))

    IF ( alamin .EQ. 0.0_wp ) THEN

       qj_rel(:) = qj_rel_nr_old(:)

       RETURN

    END IF

    alam = 1.0_wp

    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( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj , qj_old ,       &
            a_diag , coeff_f , rj_not_impl , 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.0_wp ) THEN

             tmplam = - slope / ( 2.0_wp * ( 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 - rhs2/alam2**2 ) / ( alam - alam2 )
             b = ( -alam2*rhs1/alam**2 + alam*rhs2/alam2**2 ) / ( alam - alam2 )

             IF ( a .EQ. 0.0_wp ) THEN

                tmplam = - slope / ( 2.0_wp * b )

             ELSE

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

                IF ( disc .LT. 0.0_wp ) THEN

                   tmplam = 0.5_wp * alam

                ELSE IF ( b .LE. 0.0_wp ) THEN

                   tmplam = ( - b + sqrt(disc) ) / ( 3.0_wp * a )

                ELSE

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

                ENDIF

             END IF

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

          END IF

       END IF

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

    END DO optimal_step_search

    RETURN
    
  END SUBROUTINE lnsrch

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

  SUBROUTINE eval_f( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj , qj_old ,      &
       a_diag , coeff_f , rj_not_impl , f_nl , scal_f )

    USE constitutive_2d, ONLY : eval_nonhyperbolic_terms

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: cell_fract_jk 
    REAL(wp), INTENT(IN) :: dx_rel_jk
    REAL(wp), INTENT(IN) :: dy_rel_jk

    REAL(wp), INTENT(IN) :: qj(n_vars)
    REAL(wp), INTENT(IN) :: qj_old(n_vars)
    REAL(wp), INTENT(IN) :: a_diag
    REAL(wp), INTENT(IN) :: coeff_f(n_eqns)
    REAL(wp), INTENT(IN) :: Rj_not_impl(n_eqns)

    REAL(wp), INTENT(OUT) :: f_nl(n_eqns)
    REAL(wp), INTENT(OUT) :: scal_f

    REAL(wp) :: nh_term_impl(n_eqns)
    REAL(wp) :: Rj(n_eqns)

    CALL eval_nonhyperbolic_terms( cell_fract_jk, dx_rel_jk , dy_rel_jk ,       &
         r_qj = qj , r_nh_term_impl=nh_term_impl ) 

    rj = rj_not_impl - a_diag * nh_term_impl

    f_nl = qj - qj_old + dt * rj

    f_nl = coeff_f * f_nl

    scal_f = 0.5_wp * dot_product( f_nl , f_nl )

    RETURN
    
  END SUBROUTINE eval_f

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

  SUBROUTINE eval_jacobian( cell_fract_jk , dx_rel_jk , dy_rel_jk , qj_rel ,    &
       qj_org , coeff_f, left_matrix)

    USE constitutive_2d, ONLY : eval_nonhyperbolic_terms

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: cell_fract_jk 
    REAL(wp), INTENT(IN) :: dx_rel_jk
    REAL(wp), INTENT(IN) :: dy_rel_jk

    REAL(wp), INTENT(IN) :: qj_rel(n_vars)
    REAL(wp), INTENT(IN) :: qj_org(n_vars)
    REAL(wp), INTENT(IN) :: coeff_f(n_eqns)

    REAL(wp), INTENT(OUT) :: left_matrix(n_eqns,n_vars)

    REAL(wp) :: Jacob_relax(n_eqns,n_vars)
    COMPLEX(wp) :: nh_terms_cmplx_impl(n_eqns)
    COMPLEX(wp) :: qj_cmplx(n_vars) , qj_rel_cmplx(n_vars)
    COMPLEX(wp) :: qj_rel_cmplx_init(n_vars)

    REAL(wp) :: h

    INTEGER :: i

    h = n_vars * epsilon(1.0_wp)

    ! initialize the matrix of the linearized system and the Jacobian

    left_matrix(1:n_eqns,1:n_vars) = 0.0_wp
    jacob_relax(1:n_eqns,1:n_vars) = 0.0_wp

    ! evaluate the jacobian of the non-hyperbolic terms

    DO i=1,n_vars

       qj_rel_cmplx_init(i) = cmplx(qj_rel(i),0.0_wp,wp)

    END DO

    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_cmplx_init(1:n_vars)
          qj_rel_cmplx(i) = cmplx(qj_rel(i), h,wp)

          qj_cmplx = qj_rel_cmplx * qj_org

          CALL eval_nonhyperbolic_terms( cell_fract_jk , dx_rel_jk , dy_rel_jk ,&
               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

    RETURN
    
  END SUBROUTINE eval_jacobian

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

  SUBROUTINE eval_hyperbolic_terms( q_expl , qp_expl , divFlux_iRK )

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

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: q_expl(n_vars,comp_cells_x,comp_cells_y)
    REAL(wp), INTENT(IN) :: qp_expl(n_vars+2,comp_cells_x,comp_cells_y)
    REAL(wp), INTENT(OUT) :: divFlux_iRK(n_eqns,comp_cells_x,comp_cells_y)

    INTEGER :: l , i, j, k

    ! Linear reconstruction of the physical variables at the interfaces
    CALL reconstruction(q_expl,qp_expl)

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

    ! 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

    CASE ("UP")

       CALL eval_flux_up

    END SELECT

    !WRITE(*,*) 'H_interface_x(1,2,1)',H_interface_x(1,2,1)
    !WRITE(*,*) 'q_interfaceR(1,1,1)',q_interfaceR(1,1,1)
    !WRITE(*,*) 'qp_interfaceR(1,1,1)',qp_interfaceR(1,1,1)
    !WRITE(*,*) 'H_interface_x(1,1,1)',H_interface_x(1,1,1)
    !WRITE(*,*) 'H_interface_y(1,1,2)',H_interface_y(1,1,2)
    !WRITE(*,*) 'H_interface_y(1,1,1)',H_interface_y(1,1,1)

    !$OMP PARALLEL DO private(l,j,k,i)

    ! Advance in time the solution
    cells_loop:DO l = 1,solve_cells

       j = j_cent(l)
       k = k_cent(l)

       DO i=1,n_eqns

          divflux_irk(i,j,k) = 0.0_wp

          IF ( comp_cells_x .GT. 1 ) THEN

             divflux_irk(i,j,k) = divflux_irk(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_irk(i,j,k) = divflux_irk(i,j,k) +                          &
                  ( h_interface_y(i,j,k+1) - h_interface_y(i,j,k) ) / dy

          END IF

       END DO

    END DO cells_loop

    !$OMP END PARALLEL DO

    RETURN

  END SUBROUTINE eval_hyperbolic_terms

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE eval_flux_up

    ! External procedures
    USE constitutive_2d, ONLY : eval_fluxes

    IMPLICIT NONE

    REAL(wp) :: fluxL(n_eqns)
    REAL(wp) :: fluxR(n_eqns)
    REAL(wp) :: fluxB(n_eqns)
    REAL(wp) :: fluxT(n_eqns)

    INTEGER :: j,k,l

    h_interface_x = 0.0_wp
    h_interface_y = 0.0_wp

    IF ( comp_cells_x .GT. 1 ) THEN

       DO l = 1,solve_interfaces_x

          j = j_stag_x(l)
          k = k_stag_x(l)

          CALL eval_fluxes( q_interfacel(1:n_vars,j,k) ,                        &
               qp_interfacel(1:n_vars+2,j,k) , 1 , fluxl)

          CALL eval_fluxes( q_interfacer(1:n_vars,j,k) ,                        &
               qp_interfacer(1:n_vars+2,j,k) , 1 , fluxr)

          IF ( ( q_interfacel(2,j,k) .GT. 0.0_wp ) .AND.                        &
               ( q_interfacer(2,j,k) .GE. 0.0_wp ) ) THEN

             h_interface_x(:,j,k) = fluxl

          ELSEIF ( ( q_interfacel(2,j,k) .LE. 0.0_wp ) .AND.                    &
               ( q_interfacer(2,j,k) .LT. 0.0_wp ) ) THEN

             h_interface_x(:,j,k) = fluxr

          ELSE

             h_interface_x(:,j,k) = 0.5_wp * ( fluxl + fluxr )

          END IF

          IF ( (  q_interfacel(2,j,k) .EQ. 0.0_wp ) .AND.                       &
               (  q_interfacer(2,j,k) .EQ. 0.0_wp ) ) THEN

             h_interface_x(1,j,k) = 0.0_wp
             h_interface_x(4:n_vars,j,k) = 0.0_wp

          END IF

       END DO

    END IF


    IF ( comp_cells_y .GT. 1 ) THEN

       DO l = 1,solve_interfaces_y

          j = j_stag_y(l)
          k = k_stag_y(l)

          CALL eval_fluxes( q_interfaceb(1:n_vars,j,k) ,                        &
               qp_interfaceb(1:n_vars+2,j,k) , 2 , fluxb)

          CALL eval_fluxes( q_interfacet(1:n_vars,j,k) ,                        &
               qp_interfacet(1:n_vars+2,j,k) ,2 , fluxt)

          IF ( ( q_interfaceb(3,j,k) .GT. 0.0_wp ) .AND.                        &
               ( q_interfacet(3,j,k) .GE. 0.0_wp ) ) THEN

             h_interface_y(:,j,k) = fluxb

          ELSEIF ( ( q_interfaceb(3,j,k) .LE. 0.0_wp ) .AND.                    &
               ( q_interfacet(3,j,k) .LT. 0.0_wp ) ) THEN

             h_interface_y(:,j,k) = fluxt

          ELSE

             h_interface_y(:,j,k) = 0.5_wp * ( fluxb + fluxt )

          END IF

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

             h_interface_y(1,j,k) = 0.0_wp
             h_interface_y(4:n_vars,j,k) = 0.0_wp

          END IF

       END DO

    END IF

    RETURN

  END SUBROUTINE eval_flux_up


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

  SUBROUTINE eval_flux_kt

    ! External procedures
    USE constitutive_2d, ONLY : eval_fluxes

    IMPLICIT NONE

    REAL(wp) :: fluxL(n_eqns)
    REAL(wp) :: fluxR(n_eqns)
    REAL(wp) :: fluxB(n_eqns)
    REAL(wp) :: fluxT(n_eqns)

    REAL(wp) :: flux_avg_x(n_eqns)   
    REAL(wp) :: flux_avg_y(n_eqns)   

    INTEGER :: i,j,k,l

    ! WRITE(*,*) 'eval_flux_KT: qp_interfaceR(1,1,1)',qp_interfaceR(1,1,1)


    !H_interface_x = 0.0_wp
    !H_interface_y = 0.0_wp

    !$OMP PARALLEL

    IF ( comp_cells_x .GT. 1 ) THEN

       !$OMP DO private(j,k,i,fluxL,fluxR,flux_avg_x)

       interfaces_x_loop:DO l = 1,solve_interfaces_x

          j = j_stag_x(l)
          k = k_stag_x(l)

          CALL eval_fluxes( q_interfacel(1:n_vars,j,k) ,                        &
               qp_interfacel(1:n_vars+2,j,k) , 1 , fluxl)

          CALL eval_fluxes( q_interfacer(1:n_vars,j,k) ,                        &
               qp_interfacer(1:n_vars+2,j,k) , 1 , fluxr)

          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.0_wp

             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,alphas) if the 
          ! velocities at the interfaces are null, then the flux is null
          IF ( (  qp_interfacel(2,j,k) .EQ. 0.0_wp ) .AND.                      &
               (  qp_interfacer(2,j,k) .EQ. 0.0_wp ) ) THEN

             h_interface_x(1,j,k) = 0.0_wp
             h_interface_x(4:n_vars,j,k) = 0.0_wp

          END IF

       END DO interfaces_x_loop

       !$OMP END DO NOWAIT

    END IF


    IF ( comp_cells_y .GT. 1 ) THEN

       !$OMP DO private(j,k,i,fluxB,fluxT,flux_avg_y)
       
       interfaces_y_loop:DO l = 1,solve_interfaces_y

          j = j_stag_y(l)
          k = k_stag_y(l)

          CALL eval_fluxes( q_interfaceb(1:n_vars,j,k) ,                        &
               qp_interfaceb(1:n_vars+2,j,k) , 2 , fluxb)    

          CALL eval_fluxes( q_interfacet(1:n_vars,j,k) ,                        &
               qp_interfacet(1:n_vars+2,j,k) , 2 , fluxt)

          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.0_wp

             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,alphas) if the 
          ! velocities at the interfaces are null, then the flux is null
          IF ( (  q_interfaceb(3,j,k) .EQ. 0.0_wp ) .AND.                       &
               (  q_interfacet(3,j,k) .EQ. 0.0_wp ) ) THEN

             h_interface_y(1,j,k) = 0.0_wp
             h_interface_y(4:n_vars,j,k) = 0.0_wp

          END IF

       END DO interfaces_y_loop

       !$OMP END DO

    END IF

    !$OMP END PARALLEL

    RETURN
    
  END SUBROUTINE eval_flux_kt

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

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

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: a1(:) , a2(:)
    REAL(wp), INTENT(IN) :: w1(:) , w2(:)
    REAL(wp), 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.5_wp * ( w1(i) + w2(i) )
          w_avg(i) = 0.0_wp

       ELSE

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

       END IF

    END DO

    RETURN
    
  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(q_expl,qp_expl)

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

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

    USE geometry_2d, ONLY : sourcew , sourcee , sourcen , sources
    USE geometry_2d, ONLY : sourcew_vect_x , sourcew_vect_y
    USE geometry_2d, ONLY : sourcee_vect_x , sourcee_vect_y
    USE geometry_2d, ONLY : sourcen_vect_x , sourcen_vect_y
    USE geometry_2d, ONLY : sources_vect_x , sources_vect_y

    USE parameters_2d, ONLY : reconstr_coeff

    USE geometry_2d, ONLY : minmod

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: q_expl(:,:,:)
    REAL(wp), INTENT(IN) :: qp_expl(:,:,:)

    REAL(wp) :: qrecW(n_vars+2)
    REAL(wp) :: qrecE(n_vars+2)
    REAL(wp) :: qrecS(n_vars+2)
    REAL(wp) :: qrecN(n_vars+2)

    REAL(wp) :: source_bdry(n_vars+2)
    REAL(wp) :: qrec_prime_x(n_vars+2)
    REAL(wp) :: qrec_prime_y(n_vars+2)

    REAL(wp) :: qp2recW(3) , qp2recE(3)
    REAL(wp) :: qp2recS(3) , qp2recN(3) 

    REAL(wp) :: qrec_stencil(3)
    REAL(wp) :: x_stencil(3)
    REAL(wp) :: y_stencil(3)

    INTEGER :: l,j,k
    INTEGER :: i

    REAL(wp) :: dq

    !$OMP PARALLEL DO private(j,k,i,qrecW,qrecE,qrecS,qrecN,x_stencil,y_stencil,&
    !$OMP & qrec_stencil,qrec_prime_x,qrec_prime_y,qp2recW,qp2recE,qp2recS,     &
    !$OMP & qp2recN,source_bdry,dq)

    DO l = 1,solve_cells

       j = j_cent(l)
       k = k_cent(l)

       qrecw(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
       qrece(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
       qrecs(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
       qrecn(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)

       x_stencil(2) = x_comp(j)
       y_stencil(2) = y_comp(k)

       vars_loop:DO i=1,n_vars

          qrec_stencil(2) = qp_expl(i,j,k)

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

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

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

                   x_stencil(1) = x_stag(1)
                   x_stencil(3) = x_comp(j+1)

                   qrec_stencil(1) = bcw(i)%value
                   qrec_stencil(3) = qp_expl(i,j+1,k)

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

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

                   qrec_prime_x(i) = bcw(i)%value

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

                   qrec_prime_x(i) = ( qp_expl(i,2,k) - qp_expl(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)= qp_expl(i,j-1,k)

                   x_stencil(3) = x_stag(comp_interfaces_x)
                   x_stencil(1) = x_comp(j-1)

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

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

                   qrec_prime_x(i) = bce(i)%value

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

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

                END IF

                ! internal x cells
             ELSE

                x_stencil(1) = x_comp(j-1)
                x_stencil(3) = x_comp(j+1)

                qrec_stencil(1) = qp_expl(i,j-1,k)
                qrec_stencil(3) = qp_expl(i,j+1,k)

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

             ENDIF check_x_boundary

             dq = reconstr_coeff* dx2 * qrec_prime_x(i) 

             qrecw(i) = qrec_stencil(2) - dq
             qrece(i) = qrec_stencil(2) + dq

          END IF check_comp_cells_x

          ! 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

                   y_stencil(1) = y_stag(1)
                   y_stencil(3) = y_comp(k+1)

                   qrec_stencil(1) = bcs(i)%value
                   qrec_stencil(3) = qp_expl(i,j,k+1)

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

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

                   qrec_prime_y(i) = bcs(i)%value

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

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

                END IF

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

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

                   y_stencil(3) = y_stag(comp_interfaces_y)
                   y_stencil(1) = y_comp(k-1)

                   qrec_stencil(3) = bcn(i)%value
                   qrec_stencil(1)= qp_expl(i,j,k-1)

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

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

                   qrec_prime_y(i) = bcn(i)%value

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

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

                END IF

                ! Internal y cells
             ELSE

                y_stencil(1) = y_comp(k-1)
                y_stencil(3) = y_comp(k+1)

                qrec_stencil(1) = qp_expl(i,j,k-1)
                qrec_stencil(3) = qp_expl(i,j,k+1)

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

             ENDIF check_y_boundary

             dq = reconstr_coeff* dy2 * qrec_prime_y(i)

             qrecs(i) = qrec_stencil(2) - dq
             qrecn(i) = qrec_stencil(2) + dq

          ENDIF check_comp_cells_y

       ENDDO vars_loop

       add_vars_loop:DO i=n_vars+1,n_vars+2
          ! reconstruction on u and v with same limiters of hu,hv

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

             IF ( (j .GT. 1) .AND. (j .LT. comp_cells_x) ) THEN

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

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

                dq = reconstr_coeff*dx2*qrec_prime_x(i)

                qrecw(i) = qrec_stencil(2) - dq
                qrece(i) = qrec_stencil(2) + dq

             END IF

          END IF check_comp_cells_x2

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

             IF ( ( k .GT. 1 ) .AND. ( k .LT. comp_cells_y ) ) THEN

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

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

                dq = reconstr_coeff*dy2*qrec_prime_y(i) 

                qrecs(i) = qrec_stencil(2) - dq
                qrecn(i) = qrec_stencil(2) + dq

             ENDIF

          ENDIF check_comp_cells_y2

       ENDDO add_vars_loop


       IF ( comp_cells_x .GT. 1 ) THEN

          IF ( ( j .GT. 1 ) .AND. ( j .LT. comp_cells_x ) ) THEN

             IF ( q_expl(1,j,k) .EQ. 0.0_wp ) THEN

                ! In the internal cell, if thickness h is 0 at the center
                ! of the cell, then all the variables are 0 at the center
                ! and at the interfaces (no conversion back is needed from
                ! reconstructed to conservative)
                q_interfacer(:,j,k) = 0.0_wp
                q_interfacel(:,j+1,k) = 0.0_wp

                qp_interfacer(1:3,j,k) = 0.0_wp
                qp_interfacer(4:n_vars,j,k) = qrecw(4:n_vars)
                qp_interfacer(n_vars+1:n_vars+2,j,k) = 0.0_wp

                qp_interfacel(1:3,j+1,k) = 0.0_wp
                qp_interfacel(4:n_vars,j+1,k) = qrece(4:n_vars)
                qp_interfacel(n_vars+1:n_vars+2,j+1,k) = 0.0_wp

             END IF

          END IF

          IF ( j.EQ.1 ) THEN

             ! Dirichelet boundary condition at the west of the domain
             DO i=1,n_vars

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

                   qrecw(i) = bcw(i)%value 

                END IF

             ENDDO

             DO i=n_vars+1,n_vars+2

                CALL qp_to_qp2( qrecw(1:n_vars) , qp2recw ) 
                qrecw(i) = qp2recw(i-n_vars+1)

                CALL qp_to_qp2( qrece(1:n_vars) , qp2rece ) 
                qrece(i) = qp2rece(i-n_vars+1)

             END DO

          ELSEIF ( j.EQ.comp_cells_x ) THEN

             ! Dirichelet boundary condition at the east of the domain
             DO i=1,n_vars

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

                   qrece(i) = bce(i)%value 

                END IF

             ENDDO

             DO i=n_vars+1,n_vars+2

                CALL qp_to_qp2( qrecw(1:n_vars) , qp2recw ) 
                qrecw(i) = qp2recw(i-n_vars+1)

                CALL qp_to_qp2( qrece(1:n_vars) , qp2rece ) 
                qrece(i) = qp2rece(i-n_vars+1)

             END DO

          ELSE

          END IF

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

          qp_interfacer(1:n_vars+2,j,k) = qrecw(1:n_vars+2)
          qp_interfacel(1:n_vars+2,j+1,k) = qrece(1:n_vars+2)

          IF ( j.EQ.1 ) THEN

             ! Interface value at the left of first x-interface (external)
             q_interfacel(:,j,k) = q_interfacer(:,j,k)
             qp_interfacel(:,j,k) = qp_interfacer(:,j,k)

             !WRITE(*,*) 'q_interfaceL(:,j,k)',q_interfaceL(:,j,k)
             !READ(*,*)

          ELSEIF ( j.EQ.comp_cells_x ) THEN

             ! Interface value at the right of last x-interface (external)
             q_interfacer(:,j+1,k) = q_interfacel(:,j+1,k)
             qp_interfacer(:,j+1,k) = qp_interfacel(:,j+1,k)

          ELSE

          END IF

       ELSE

          ! for case comp_cells_x = 1 
          q_interfacer(1:n_vars,j,k) = q_expl(1:n_vars,j,k)
          q_interfacel(1:n_vars,j+1,k) = q_expl(1:n_vars,j,k)

          qp_interfacer(1:n_vars+2,j,k) = qp_expl(1:n_vars+2,j,k)
          qp_interfacel(1:n_vars+2,j+1,k) = qp_expl(1:n_vars+2,j,k)

       END IF

       IF ( comp_cells_y .GT. 1 ) THEN

          IF ( ( k .GT. 1 ) .AND. ( k .LT. comp_cells_y ) ) THEN

             IF ( q_expl(1,j,k) .EQ. 0.0_wp ) THEN

                ! In the internal cell, if thickness h is 0 at the center
                ! of the cell, then all the variables are 0 at the center
                ! and at the interfaces (no conversion back is needed from
                ! reconstructed to conservative)

                q_interfacet(:,j,k) = 0.0_wp
                q_interfaceb(:,j,k+1) = 0.0_wp

                qp_interfacet(1:3,j,k) = 0.0_wp
                qp_interfacet(4:n_vars,j,k) = qrecs(4:n_vars)
                qp_interfacet(n_vars+1:n_vars+2,j,k) = 0.0_wp

                qp_interfaceb(1:3,j,k+1) = 0.0_wp
                qp_interfaceb(4:n_vars,j,k+1) = qrecn(4:n_vars)
                qp_interfaceb(n_vars+1:n_vars+2,j,k+1) = 0.0_wp

             END IF

          END IF

          IF ( k .EQ. 1 ) THEN

             ! Dirichelet boundary condition at the south of the domain
             DO i=1,n_vars

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

                   qrecs(i) = bcs(i)%value 

                END IF

             ENDDO

             DO i=n_vars+1,n_vars+2

                CALL qp_to_qp2( qrecs(1:n_vars) , qp2recs ) 
                qrecs(i) = qp2recs(i-n_vars+1)

                CALL qp_to_qp2( qrecn(1:n_vars) , qp2recn ) 
                qrecn(i) = qp2recn(i-n_vars+1)

             END DO

          ELSEIF ( k .EQ. comp_cells_y ) THEN

             ! Dirichelet boundary condition at the north of the domain
             DO i=1,n_vars

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

                   qrecn(i) = bcn(i)%value 

                END IF

             ENDDO

             DO i=n_vars+1,n_vars+2

                CALL qp_to_qp2( qrecs(1:n_vars) , qp2recs ) 
                qrecs(i) = qp2recs(i-n_vars+1)

                CALL qp_to_qp2( qrecn(1:n_vars) , qp2recn ) 
                qrecn(i) = qp2recn(i-n_vars+1)

             END DO

          ELSE

          END IF

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

          qp_interfacet(1:n_vars+2,j,k) = qrecs(1:n_vars+2)
          qp_interfaceb(1:n_vars+2,j,k+1) = qrecn(1:n_vars+2)

          IF ( k .EQ. 1 ) THEN

             ! Interface value at the bottom of first y-interface (external)
             q_interfaceb(:,j,k) = q_interfacet(:,j,k)
             qp_interfaceb(:,j,k) = qp_interfacet(:,j,k)

          ELSEIF ( k .EQ. comp_cells_y ) THEN

             ! Interface value at the top of last y-interface (external)
             q_interfacet(:,j,k+1) = q_interfaceb(:,j,k+1)
             qp_interfacet(:,j,k+1) = qp_interfaceb(:,j,k+1)

          ELSE

          END IF

       ELSE

          ! case comp_cells_y = 1

          q_interfaceb(:,j,k) = q_expl(:,j,k)
          q_interfacet(:,j,k) = q_expl(:,j,k)
          q_interfaceb(:,j,k+1) = q_expl(:,j,k)
          q_interfacet(:,j,k+1) = q_expl(:,j,k)

          qp_interfaceb(:,j,k) = qp_expl(:,j,k)
          qp_interfacet(:,j,k) = qp_expl(:,j,k)
          qp_interfaceb(:,j,k+1) = qp_expl(:,j,k)
          qp_interfacet(:,j,k+1) = qp_expl(:,j,k)

       END IF

    END DO

    !$OMP END PARALLEL DO

    RETURN

  END SUBROUTINE reconstruction


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

  SUBROUTINE eval_speeds

    ! External procedures
    USE constitutive_2d, ONLY : eval_local_speeds_x, eval_local_speeds_y 

    IMPLICIT NONE

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

    INTEGER :: j,k,l

    !$OMP PARALLEL

    IF ( comp_cells_x .GT. 1 ) THEN

       !$OMP DO private(j , k , abslambdaL_min , abslambdaL_max ,               &
       !$OMP & abslambdaR_min , abslambdaR_max , min_r , max_r )

       x_interfaces_loop:DO l = 1,solve_interfaces_x

          j = j_stag_x(l)
          k = k_stag_x(l)

          CALL eval_local_speeds_x( qp_interfacel(:,j,k) , abslambdal_min ,     &
               abslambdal_max )

          CALL eval_local_speeds_x( qp_interfacer(:,j,k) , abslambdar_min ,     &
               abslambdar_max )

          min_r = min(abslambdal_min , abslambdar_min , 0.0_wp)
          max_r = max(abslambdal_max , abslambdar_max , 0.0_wp)

          a_interface_xneg(:,j,k) = min_r
          a_interface_xpos(:,j,k) = max_r

       END DO x_interfaces_loop

       !$OMP END DO NOWAIT

    END IF

    IF ( comp_cells_y .GT. 1 ) THEN

       !$OMP DO private(j , k , abslambdaB_min , abslambdaB_max ,               &
       !$OMP & abslambdaT_min , abslambdaT_max , min_r , max_r )

       y_interfaces_loop:DO l = 1,solve_interfaces_y

          j = j_stag_y(l)
          k = k_stag_y(l)

          CALL eval_local_speeds_y( qp_interfaceb(:,j,k) , abslambdab_min ,     &
               abslambdab_max )

          CALL eval_local_speeds_y( qp_interfacet(:,j,k) , abslambdat_min ,     &
               abslambdat_max )

          min_r = min(abslambdab_min , abslambdat_min , 0.0_wp)
          max_r = max(abslambdab_max , abslambdat_max , 0.0_wp)

          a_interface_yneg(:,j,k) = min_r
          a_interface_ypos(:,j,k) = max_r

       END DO y_interfaces_loop

       !$OMP END DO

    END IF

    !$OMP END PARALLEL

    RETURN
    
  END SUBROUTINE eval_speeds

END MODULE solver_2d