save_1D.f90

### TABLE OF CONTENTS

- 1. NUM/save_1D

## NUM/save_1D [ Modules ]

[ Top ] [ Modules ]

NOM

     save_1D(calc)

DESCRIPTION

     Posttraitement dans le cas ou le maillage est 1D
     Post-processing if the mesh is 1D

        ENTREES / INPUT
      calc: objet de type calcul / object of type calcul

        SORTIES / OUTPUT

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
      IMPLICIT NONE
    !.2-----  Declaration
      type(calcul), intent(inout)         :: calc !! Objet calcul / Calcul object

      real(kind=kind(0.d0))               :: r,u,p,psiL,psiR,psi,xi,alf,err1,err2,volt
      character*7                         :: fsol
      character*10                        :: cn
      integer                             :: i,ic,id,iu,ios
      character(len=256)                  :: iomsg
    !=========================== DEBUT DU CODE EXECUTABLE ==================
    !-------  Initialisations
      cn='0123456789'
      ic=int(calc%msh%numdom/100.)
      id=int((calc%msh%numdom-ic*100.)/10.)
      iu=calc%msh%numdom-id*10-ic*100

      IF (calc%msh%numdom.gt.999) call print_err('save_1D','Domain index > 999 is not supported')
      fsol(5:5)=cn(ic+1:ic+1)
      fsol(6:6)=cn(id+1:id+1)
      fsol(7:7)=cn(iu+1:iu+1)

    ! ouverture des fichiers
      SELECT CASE (imodel_phy)

        CASE(2)
          fsol(1:4)='n_r_'  ! densite
          OPEN(41,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='n_u_'  ! vitesse
          OPEN(42,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='n_p_'  ! pression
          OPEN(43,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='n_e_'  ! production numerique d entropie
          OPEN(49,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='n_c_'  ! fonction couleur
          OPEN(44,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='e_r_'
          OPEN(51,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='e_u_'
          OPEN(52,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='e_p_'
          OPEN(53,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='e_e_'
          open(54,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))

        CASE(1)
          fsol(1:4)='n_h_'  ! hauteur d eau
          OPEN(41,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='n_u_'  ! vitesse
          OPEN(42,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='n_e_'  ! production numerique d entropie
          OPEN(49,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='e_h_'
          OPEN(51,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          fsol(1:4)='e_u_'
          OPEN(52,file=fsol,iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open '//trim(fsol)//' : '//trim(iomsg))
          open(53,file='cerf_1d.dat',iostat=ios,iomsg=iomsg)
          if (ios /= 0) call print_err('save_1D','Cannot open cerf_1d.dat : '//trim(iomsg))
          write(53,*,iostat=ios,iomsg=iomsg)'TITLE="Solution Numerique CERF"'
          if (ios /= 0) call print_err('save_1D','Error writing cerf_1d header: '//trim(iomsg))
          write(53,*,iostat=ios,iomsg=iomsg)' VARIABLES="X","h","U","zb","eta"'
          if (ios /= 0) call print_err('save_1D','Error writing cerf_1d variables: '//trim(iomsg))
          write(53,111,iostat=ios,iomsg=iomsg)calc%msh%nb_cell
          if (ios /= 0) call print_err('save_1D','Error writing cerf_1d zone: '//trim(iomsg))
    111     format('ZONE T="", I=',i9,' F=POINT')

        CASE default
          call print_err('save_1D','Physical model not defined')

        END SELECT

      ! remplissage des fichiers de resultats num
      err1=0.d0
      err2=0.d0
      volt=0.d0
      if (calc%tmax == 0.d0) call print_err('save_1D','tmax is zero, xi cannot be computed')
      DO i=1,calc%msh%nb_cell
        xi=calc%msh%list_cell(i)%center%x/calc%tmax
        CALL Bal2Prim(calc%msh%list_cell(i)%w)
        volt=volt+calc%msh%list_cell(i)%vol

        SELECT CASE (imodel_phy)

          CASE(2)
            write(41,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(1)! densite
            if (ios /= 0) call print_err('save_1D','Error writing n_r_: '//trim(iomsg))
            write(42,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(2)! vitesse
            if (ios /= 0) call print_err('save_1D','Error writing n_u_: '//trim(iomsg))
            write(43,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(5)! pression
            if (ios /= 0) call print_err('save_1D','Error writing n_p_: '//trim(iomsg))
            write(49,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(6)
            if (ios /= 0) call print_err('save_1D','Error writing n_e_: '//trim(iomsg))
            write(44,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vbal%v(5)
            if (ios /= 0) call print_err('save_1D','Error writing n_c_: '//trim(iomsg))
            psiL=phy_state(1)%vbal%v(5)
            psiR=phy_state(2)%vbal%v(5)
            CALL riemisot(phy_state(1)%vprim%v(1),phy_state(1)%vprim%v(2),psiL,psiL, &
            phy_state(2)%vprim%v(1),phy_state(2)%vprim%v(2),psiR,psiR,xi,r,u,p,psi,alf)

            write(51,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,r
            if (ios /= 0) call print_err('save_1D','Error writing e_r_: '//trim(iomsg))
            write(52,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,u
            if (ios /= 0) call print_err('save_1D','Error writing e_u_: '//trim(iomsg))
            write(53,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,p
            if (ios /= 0) call print_err('save_1D','Error writing e_p_: '//trim(iomsg))
            write(54,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,psi
            if (ios /= 0) call print_err('save_1D','Error writing e_e_: '//trim(iomsg))

          CASE(1)
            write(41,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(1)+&
                                               calc%msh%list_cell(i)%w%vprim%v(5)! hauteur d eau
            if (ios /= 0) call print_err('save_1D','Error writing n_h_: '//trim(iomsg))
            write(42,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(2)! vitesse
            if (ios /= 0) call print_err('save_1D','Error writing n_u_: '//trim(iomsg))
            write(49,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(4)! production entropie
            if (ios /= 0) call print_err('save_1D','Error writing n_e_: '//trim(iomsg))
            CALL svriemann(phy_state(1)%vprim%v(1),phy_state(1)%vprim%v(2),&
                           phy_state(2)%vprim%v(1),phy_state(2)%vprim%v(2),xi,r,u)
            write(51,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,r
            if (ios /= 0) call print_err('save_1D','Error writing e_h_: '//trim(iomsg))
            write(52,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,u
            if (ios /= 0) call print_err('save_1D','Error writing e_u_: '//trim(iomsg))
            write(53,*,iostat=ios,iomsg=iomsg) calc%msh%list_cell(i)%center%x,calc%msh%list_cell(i)%w%vprim%v(1),&
                                               calc%msh%list_cell(i)%w%vprim%v(2), &
                                               calc%msh%list_cell(i)%w%vprim%v(5), &
                                               calc%msh%list_cell(i)%w%vprim%v(1)+calc%msh%list_cell(i)%w%vprim%v(5)
            if (ios /= 0) call print_err('save_1D','Error writing cerf_1d.dat row: '//trim(iomsg))

          CASE default
             call print_err('save_1D','Physical model not defined')

          END SELECT
          err1=err1+calc%msh%list_cell(i)%vol*(calc%msh%list_cell(i)%w%vprim%v(1)-r)**2
          err2=err2+calc%msh%list_cell(i)%vol*(calc%msh%list_cell(i)%w%vprim%v(2)-u)**2

      END DO
      if (volt <= 0.d0) call print_err('save_1D','Total cell volume is zero, cannot compute L2 errors')
      err1=sqrt(err1/volt)
      err2=sqrt(err2/volt)

    ! Fermeture des fichiers
      SELECT CASE (imodel_phy)

        CASE(2)
        write( *,1331)err1
        write( *,1332)err2
    1331 format('Density L2 error  :',e13.5)
    1332 format('Velocity L2 error :',e13.5)
          close(41)
          close(42)
          close(43)
          close(49)
          close(44)
          close(51)
          close(52)
          close(53)
          close(54)

        CASE(1)
        write( *,1001)err1
        write( *,1002)err2
    1001 format('Water level L2 error  :',e13.5)
    1002 format('Velocity L2 error     :',e13.5)
          close(41)
          close(42)
          close(49)
          close(51)
          close(52)
          close(53)

        CASE default
          call print_err('save_1D','Unexpected physical model case')

        END SELECT

    !===========================   FIN DE LA ROUTINE    ====================
    END SUBROUTINE save_1D