 SUBROUTINE rk2(calc)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    NUM/rk2  *
!
!    NOM
!   rk2(calc)
!    DESCRIPTION  
!*
! Integration en temps, par la methode de Runge-Kutta d ordre 1 ou 2 des bilans des flux.
!
! Integration in time, by the Runge-Kutta method of order 1 or 2 of flux balances.
!*
! 
!    ENTREES / INPUT
!   calc : objet calcul / calculation object
!    SORTIES / OUTPUT
!   calc : objet calcul mis a jour / updated calculation object
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use num_typ
  use num ,only:  dtloc, synchro_mpi, u_penalisation,display,calc_flux, calc_source,synchro_w,sharp
  IMPLICIT NONE
!.2-----  Declaration
  type(calcul),intent(inout) :: calc !! Objet calcul / Calcul object
  
  integer                    :: i,nitmax
  real(kind=kind(0.d0))      :: pm,e,vol
  integer                    :: lev,ic1,ic2

!=========================== DEBUT DU CODE EXECUTABLE ==================

! level: niveau de raffinement a considerer. Ici lev=grand.
  lev=1000
  
! Temps courant  
  tc=calc%tmin
  
! Nbre de pas de temps max  
  nitmax=100000

! Boucle en temps
  DO WHILE (tc.lt.calc%tmax)

! Determination du pas de temps du calcul
      CALL dtloc(calc)

! Ajustement du pas de temps si necessaire  au dernier pas de temps
      IF (tc+calc%dt.gt.calc%tmax) calc%dt=calc%tmax-tc
 
! Synchronisation du pas de temps et des cellules aux frontieres du domaine
      CALL synchro_mpi(calc)

! Calcul des gradients si necessaire
      IF (calc%iordre_e.eq.2) then
          CALL gradprim(calc%msh)
          SELECT CASE (calc%ilimiteur)
              CASE (0)   ! barth
                  CALL gradlim(calc%msh)
              CASE (1)   ! cartesien
                  CALL gradlimxy(calc%msh)
          END SELECT
          CALL synchro_w(calc)  
      END IF
  
! Calcul de l entropie a l instant t 
      DO i=1,calc%msh%nb_cell
          CALL entropie(calc%msh%list_cell(i)%w,calc%msh%list_cell(i)%w%enti)
      END DO

!------- Initialisation du stockage des flux et sources
      DO i=1,calc%msh%nb_cell
          calc%msh%list_cell(i)%balance%v=0.d0
      END DO

!------- Calcul des flux aux faces (et eventuellement stockage)
      CALL calc_flux(calc,lev)

!------- Calcul des termes source (et eventuellement stockage)
      CALL calc_source(calc,lev)

      tc=tc+calc%dt
      calc%ipas=calc%ipas+1
          
      IF (calc%msh%numdom.eq.1) write( *,1000) calc%ipas,tc,calc%dt,calc%vmax 
1000  format(' Ipas: ',i7,' T=',e12.5,' dt=',e12.5,' vmax=',e12.5)

! Mise a jour des variables conservatives
      pm=0.d0     
      DO i=1,calc%msh%nb_cell
          vol=calc%msh%list_cell(i)%vol
          calc%msh%list_cell(i)%w%vbal%v=calc%msh%list_cell(i)%w%vbal%v   & 
          - calc%dt/calc%iordre_t/vol * calc%msh%list_cell(i)%balance%v   ! point milieu
          CALL bal2prim(calc%msh%list_cell(i)%w)
          IF (dabs(calc%dt/vol * calc%msh%list_cell(i)%balance%v(5)).gt.pm& !!!!!! PEUT-ON LE METTRE AILLEURS ?
           .and. (imodel_phy.eq.2)) pm=dabs(calc%dt/vol * calc%msh%list_cell(i)%balance%v(5))
      END DO

if(imodel_phy.eq.1)then
    calc%dt=calc%dt/calc%iordre_t
    call friction(calc)
    calc%dt=calc%dt*calc%iordre_t
endif

! Correction par le point milieu ou heun si ordre 2 en temps
      IF (calc%iordre_t.eq.2) THEN
          IF (calc%iordre_e.eq.2) THEN
              CALL synchro_w(calc)
              CALL gradprim(calc%msh)  
              SELECT CASE(calc%ilimiteur)
                  CASE(0)   ! barth
                      CALL gradlim(calc%msh)
                  CASE(1)   ! cartesien
                      CALL gradlimxy(calc%msh)
              END SELECT
              CALL synchro_w(calc)  
          END IF

! On revient en arriere d un demi pas de temps (en variables conservatives)
          DO i=1,calc%msh%nb_cell
              calc%msh%list_cell(i)%balance%v=-0.5d0*calc%msh%list_cell(i)%balance%v 
              IF (ive_phy.ne.0) calc%msh%list_cell(i)%balance%v(ive_phy)=0.d0 
          END DO
              
!------- Calcul des flux aux faces (et eventuellement stockage)
          CALL calc_flux(calc,lev)
  
!------- Calcul des termes source (et eventuellement stockage)
          CALL calc_source(calc,lev)
   
          pm=0.d0     
          DO i=1,calc%msh%nb_cell
              vol=calc%msh%list_cell(i)%vol
              calc%msh%list_cell(i)%w%vbal%v=calc%msh%list_cell(i)%w%vbal%v & 
              - calc%dt/vol * calc%msh%list_cell(i)%balance%v
              CALL bal2prim(calc%msh%list_cell(i)%w)
              IF (dabs(calc%dt/vol * calc%msh%list_cell(i)%balance%v(5)).gt.pm&
              .and. imodel_phy.eq.2) pm=dabs(calc%dt/vol *&
              calc%msh%list_cell(i)%balance%v(5))
          END DO
        if(imodel_phy.eq.1)then
            calc%dt=calc%dt/calc%iordre_t
            call friction(calc)
            calc%dt=calc%dt*calc%iordre_t
        endif
      END IF

!------- Calcul de l entropie a l instant t+dt 
      DO i=1,calc%msh%nb_cell
          CALL entropie(calc%msh%list_cell(i)%w,e)
          calc%msh%list_cell(i)%w%entf=e
          IF (ive_phy.ne.0) THEN
              calc%msh%list_cell(i)%w%pe=(e-calc%msh%list_cell(i)%w%enti)/calc%dt +&
                                         calc%msh%list_cell(i)%balance%v(ive_phy)/calc%msh%list_cell(i)%vol
              calc%msh%list_cell(i)%w%vprim%v(ive_phy)=calc%msh%list_cell(i)%w%pe
          END IF
      END DO
! dans le cas de saint venant determination de la shoreline et calcul de la friction
if(imodel_phy.eq.1)then
      DO i=1,calc%msh%nb_cell
          calc%msh%list_cell(i)%w%chi=0.d0
      END DO
      DO i=1,calc%msh%nb_face
          ic1=calc%msh%list_face(i)%ic1
          ic2=calc%msh%list_face(i)%ic2
          if(ic2.gt.0)then
              if(calc%msh%list_cell(ic1)%w%vprim%v(1).gt.hzero_phy.and.&
                 calc%msh%list_cell(ic2)%w%vprim%v(1).lt.hzero_phy)calc%msh%list_cell(ic1)%w%chi=1.d0
              if(calc%msh%list_cell(ic1)%w%vprim%v(1).lt.hzero_phy.and.&
                 calc%msh%list_cell(ic2)%w%vprim%v(1).gt.hzero_phy)calc%msh%list_cell(ic2)%w%chi=1.d0
               if(calc%msh%list_cell(ic1)%w%vprim%v(1).gt.hzero_phy.and.&
                 calc%msh%list_cell(ic2)%w%vprim%v(1).lt.hzero_phy)calc%msh%list_cell(ic2)%w%chi=1.d0
              if(calc%msh%list_cell(ic1)%w%vprim%v(1).lt.hzero_phy.and.&
                 calc%msh%list_cell(ic2)%w%vprim%v(1).gt.hzero_phy)calc%msh%list_cell(ic1)%w%chi=1.d0
         endif
      END DO
      !call friction(calc)
 endif

! Compression de l interface
      IF (imodel_phy.eq.2.and.sharp_phy.ne.0.d0) CALL sharp(calc,pm)

! Penalisation si obstacle
      IF (calc%nb_sol .NE. 0) CALL u_penalisation(calc)
! Verification conservation
      !call  info_cons(calc)

  END DO
  !if(calc%msh%numdom.eq.1.and.calc%nb_sol .NE. 0)write(13,'(7e13.5)')tc,calc%list_sol(1)%xg%v
  CALL synchro_w(calc)
  !call display(calc%msh)

!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE RK2
!DOCFIN