rk2.f90

### TABLE OF CONTENTS

- 1. NUM/rk2

## NUM/rk2 [ Modules ]

[ Top ] [ Modules ]

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