source: trunk/SOURCES/out_horiz_mod.f90 @ 23

!> \file out_horiz_mod.f90
!! Module avec les routines d'ecriture et de lecture de fichiers horizontaux
!<

!> \namespace out_hz
!! This module gathers routines to read and write horizontal files
!! \author ...
!! \date ...
!! @note Used module
!! @note   - use geography
!! @note   - use runparam
!<

module out_hz

use geography     ! permet d'avoir nx et ny et geoplace
use runparam       ! permet d'avoir tbegin,tend,runname,dirout
 
implicit none

integer, parameter :: ncol=50           !< nombre maxi de colonnes de sortie
integer,parameter :: nvar=100           !< nombre maxi de variables dans LISTE-VAR-HZ.dat
integer  :: ndtsortie   !< nombre de dtsortie
integer  :: npredeft    !< nombre de temps de sortie predefinis
integer  :: iglob_hz=0
integer :: ivar         !< index representant le numer d'une variable dans  LISTE-VAR-HZ.dat 
integer :: npos=0       !< position dans xxx, initialise a 0


! tableaux dont l'indice est npos (position en numero de colonne)

real,dimension(nx,ny,ncol) ::  xxx         !< tableau qui va contenir les variables à sortir
integer,dimension(ncol) :: varnumber    !< numero de la variable en position npos
character(len=8),dimension(ncol) :: formatcol !< le format de sortie pour une colonne
character(len=6),dimension(ncol) :: colname   !< le nom de la variable pour une colonne


! Pour chaque variable (definies dans  LISTE-VAR-HZ.dat)
! tableaux "caracteristiques de sortie" attachées a la variable
! indice dans ces tableaux : ivar=varnumber(npos) 

integer,dimension(nvar) :: isortie=0    !< si isortie=0, pas de sortie du tout.
integer,dimension(nvar) :: isort_time   !< 1 si sortie au temps time
integer,dimension(nvar) :: interv       !< entier qui code quel dtsortie utiliser
real,dimension(nvar) :: dtsortvar       !< pas de temps de sortie de chaque variable

real,dimension(nvar) :: coefsortvar     !< coefficient multiplicateur des sorties
real,dimension(nvar) :: minvar          !< valeur minimu de la variable    !pour éviter des ***
real,dimension(nvar) :: maxvar          !< valeur maximum de la variable   ! dans le format f
character(len=6),dimension(nvar) :: varname   !< le nom de la variable 
character(len=8),dimension(nvar) :: formatvar !< le format de sortie



double precision, dimension(:),allocatable :: dtsortie_hz  !< tableau des dtsortie : dimension (ndtsortie)
real,dimension(:),allocatable :: predef_tsort  !< tableau des temps predefinis pour sorties :
                                               !< dimension (npredft)



character(len=10) :: comment  
character (len=6) :: varchar

         
contains

!__________________________________________________________________________
!> SUBROUTINE: init_out_hz
!! Initialise les tableaux pour les sorties horizontaux
!>
  subroutine init_out_hz

    implicit none
    integer :: err   !< recuperation d'erreur 
    integer :: ivar
    integer :: i2
    integer :: i3
    integer  :: i,j,k !< indices de travail

    integer :: num_dat = 21

    ! initialise les tableaux 
    !----------------------------
    ! dtsortie_hz, predef_tsort  
    ! isortie,interv,dtsortvar,coefsortvar,varname,formatvar

    ! lecture des pas de temps de sortie
    !------------------------------------
    ! open(num_dat,file='../'//trim(dirsource)//'/Fichiers-parametres/TEMPS-HZ.dat')
    open(num_dat,file=trim(dirsource)//'/Fichiers-parametres/TEMPS-HZ.dat')

    ! passe les commentaires qui se terminent par une ligne de ~~~
    comment1: do k=1,500
       read(num_dat,'(a10)') comment

       if (comment.eq.'~~~~~~~~~~') exit comment1
    end do comment1


    ! lecture frequences de sortie
    read(num_dat,*) ndtsortie

    if (.not.allocated(dtsortie_hz)) THEN
       allocate(dtsortie_hz(ndtsortie),stat=err)
       if (err/=0) then
          print *,"Erreur à l'allocation du tableau dtsortie_hz ",err
          stop 4
       end if
    end if

    do k=1,ndtsortie
       read(num_dat,*) dtsortie_hz(k)
    end do

    read(num_dat,*)  ! saute une ligne

    ! lecture pas de temps predefinis
    read(num_dat,*) npredeft


    if (.not.allocated(predef_tsort)) THEN
       allocate(predef_tsort(npredeft),stat=err)
       if (err/=0) then
          print *,"Erreur à l'allocation du tableau dt-out_hz ",err
          stop 4
       end if
    end if

    do k=1,npredeft
       read(num_dat,*) predef_tsort(k)
    end do
    close(num_dat)

    ! Lecture des variables et de leur frequence de sortie
    !-----------------------------------------------------------

    ! open(num_dat,file='../'//trim(dirsource)//'/Fichiers-parametres/LISTE-VAR-HZ.dat')
    open(num_dat,file=trim(dirsource)//'/Fichiers-parametres/LISTE-VAR-HZ.dat')

    !saute les commentaires
    comment2: do k=1,500
       read(num_dat,'(a10)') comment
       if (comment.eq.'~~~~~~~~~~') exit comment2
    end do comment2


    do k=1,100
       read(num_dat,'(a6)',end=500,err=500) varchar
       read(num_dat,*,end=500,err=500) ivar,i2,i3

       varname(ivar)=' '//varchar
       isortie(ivar)=i2
       interv(ivar)=i3

       if ((i3.gt.0).and.(i3.le.ndtsortie)) then
          dtsortvar(ivar)=dtsortie_hz(i3)
       else 
          dtsortvar(ivar)=1.e10
       endif

       read(num_dat,'(a8)',end=500,err=500) formatvar(ivar)
       !  call minmax_format(minvar(ivar),maxvar(ivar),formatvar(ivar))

       read(num_dat,*,end=500,err=500) coefsortvar(ivar)
       !  print*,'k=',k,ivar,varname(ivar),interv(ivar),dtsortvar(ivar),formatvar(ivar),coefsortvar(ivar)
       read(num_dat,*,end=500,err=500)
    end do

    goto 510
500 continue
    !    write(6,*) 'nombre de variables dans liste_var',k
510 continue


    close (num_dat)
    return
  end subroutine init_out_hz

!> SUBROUTINE: testsort_time
!! Teste variable par variable si la sortie hz est faite à ce temps là
!! \param tsortie temps de sortie
!>

subroutine testsort_time(tsortie)

implicit none
double precision :: dbltime
real :: tsortie
real :: difftime  !< difference  tsortie-predef_tsort(npr)
real :: debtime   !< difference abs(tsortie-tbegin)
real :: fintime   !< difference abs(tsortie-tend)
integer   :: ipredef
integer   :: ideb
integer   :: ifin
integer   :: npr
integer  :: i,j,k ! indices de travail


!
!           exemple.  si dt_out_hz=(1000,5000,10000) 
!                    interv=2        la sortie se fera tous les 5000 ans
!                    interv=0        la sortie se fera seulement sur les pas de temps predefinis
!                    interv=-1       la sortie ne se fait qu'aux premier  pas de temps
!                    interv=-2       la sortie ne se fait qu'au premier et au dernier pas de temps


        isort_time(:)=0
        dbltime=dble(tsortie)
! recherche si ce pas de temps est un pas de temps predefini
        ipredef=0
        ideb=0
        ifin=0
        

predef:  do npr=1,npredeft
           difftime=abs(tsortie-predef_tsort(npr))
           if (difftime.lt.dttest) then 
              ipredef=1
              exit predef
           end if
           debtime=abs(tsortie-tbegin)
           fintime=abs(tsortie-tend)

           if ((debtime.lt.dttest).or.(nt.eq.1)) ideb=1
           if (fintime.lt.dttest) ifin=1

        end do predef

! boucle sur les numeros de variables 
boucle_var: do i=1,nvar

                if (isortie(i).eq.0) then  ! variables non attribuees et
                                           ! variables ou isortie est explicitement 0
                   isort_time(i)=0

                else                       ! variables dont on veut la sortie

                   if ((interv(i).ge.0).and.(ipredef.eq.1)) then   ! pas de temps predefini
                     isort_time(i)=1

                   else if ((interv(i).le.-1).and.(ideb.eq.1)) then ! sortie a Tbegin
                      isort_time(i)=1
                      
                   else if ((interv(i).eq.-2).and.(ifin.eq.1)) then ! sortie a Tend
                      isort_time(i)=1

                   else if (mod(abs(dbltime),dtsortvar(i)).lt.dble(dttest)) then
                      isort_time(i)=1

                                            ! le test est en dble car quand le temps est tres
                                            ! grand, on peut avoir des problemes d'arrondi
                   endif

                endif

             end do  boucle_var

! initialise npos et iglob_hz
     npos=0
     iglob_hz=maxval(isort_time)

return

end subroutine testsort_time
!--------------------------------------------------------------------------
!> SUBROUTINE: rempli_xxx
!! Rempli la colonne npos des tableaux xxx,varnumber,formatcol,colname
!! \param numvar  Le numero de la variable
!! \param Var     Nom de variable Var dans la liste LISTE-VAR-HZ.dat
!>
subroutine rempli_xxx(numvar,Var)
! rempli la colonne npos des tableaux xxx,varnumber,formatcol,colname
! numvar est le numero de la variable Var dans LISTE-VAR-HZ.dat
! 

implicit none

integer :: numvar
real,dimension(nx,ny) :: var
real :: mincol
real :: maxcol
real :: coef
 integer  :: i,j,k ! indices de travail

npos=npos+1
coef=coefsortvar(numvar) ! coefficient multiplicateur
mincol=minvar(numvar)
maxcol=maxvar(numvar)
varnumber(npos)=numvar
formatcol(npos)=formatvar(numvar)
colname(npos)=varname(numvar)


xxx(:,:,npos)=var(:,:)*coef 

! applique les minmax
!do j=1,ny
!   do i=1,nx
!      xxx(i,j,npos)=min(xxx(i,j,npos),maxcol)
!      xxx(i,j,npos)=max(xxx(i,j,npos),mincol)
!   end do
!end do   

return
end subroutine rempli_xxx
!---------------------------------------------------------------------------
!> SUBROUTINE: hz_output()
!! Remplace plotoutput. Ecrit le tableau xxx dans un fichier avec le nom runname//snapname//.hz
!! \param tsortie  temps de sortie 
!>

subroutine hz_output(tsortie)

 implicit none
 real tsortie

 character(len=1000) :: filout
 character(len=4) :: sep      !< pour le format de sortie
 character(len=1) :: fin      !< pour le format de sortie
 character(len=1) :: deb      !< pour le format de sortie
 character(len=1000) :: fmtxxx       !< pour le format de sortie xxx
 character(len=1000) :: fmtcolname   !< pour le format colname
 character(len=1000) :: fmtvarnumber !< pour le format varnumber
 character(len=3) :: charncol   
 character(len=30) :: snapname
 integer  :: i,j,k ! indices de travail
      integer :: num_forc = 20

!write(6,*) 'hz_output  time=', tsortie,'npos=',npos

if (npos.eq.0) goto 900

! nom du fichier
call snaptime(tsortie,snapname)
      filout =trim(runname)//trim(snapname)//'.hz'
      filout = TRIM(DIRNAMEOUT)//TRIM(filout)

!write(6,*) 'sortie hz pour time=',tsortie,'nb colonnes=',npos

open(num_forc,file=filout)

! ecriture de la ligne format pour xxx
 sep=',1x,'
 deb='('
 fin=')'
 fmtxxx=deb

 do k=1,npos-1
  fmtxxx=trim(fmtxxx)//trim(formatcol(k))//sep 
 end do

 fmtxxx=trim(fmtxxx)//trim(formatcol(npos))//fin

! met npos dans un character pour faire le format

 write(charncol,fmt='(i3)') npos
 charncol=adjustl( charncol) ! justifie a gauche

! format pour varname
 fmtcolname=deb//trim(charncol)//deb//'a6'//',1x),1x'//fin

! format pour varnumber
 fmtvarnumber=deb//trim(charncol)//deb//'i3'//',1x),1x'//fin

 
! ecriture dans le fichier sortie
  write(num_forc,*) tsortie, geoplace, '       time, geoplace'
  write(num_forc,'(10(i0,1x),a46)') nx*ny,npos,nx,ny,nint(dx/1000.),nint(seasea),xmin,xmax,ymin,ymax,& 
            'nx*ny,ncol,nx,ny,dx,sealev,xmin,xmax,ymin,ymax'
  write(num_forc,fmt=trim(fmtvarnumber)) (varnumber(k),k=1,npos)
  write(num_forc,fmt=trim(fmtcolname)) (colname(k),k=1,npos)

 do j=1,ny
  do i=1,nx
     write(num_forc,fmt=trim(fmtxxx)) (xxx(i,j,k),k=1,npos)
  end do
 end do
 close(num_forc)

900 continue
 return

end subroutine  hz_output
 

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


end module out_hz