! smooth_EOS.f
! Fortran program to smooth EOS tables calculated using SESAME (entropy only)
! Original data are from Saumon & Guillot, ApJ (2004)
! T. Guillot
!
! USE: compile with gfortran, run
! INPUT: files data/EOSH/S_SESAME_alt_waterP-q.dat
! OUPUT: files data/EOSH/S_SESAME_smo_waterP-q.dat
! No options, code itself is to be modified as needed.

      program extend_EOS
!     Version: 01-Jun-2009

      implicit none

      integer pnx,pny,pxt,pyt,dims
      parameter (pnx=101,pny=101)
      integer nx,ny,nnx,nny,naddx,naddy
      real*8 logp(pnx),logq(pny),logro(pnx*pny),logs(pnx*pny),          &
     &     newlogp(pnx),newlogq(pny),newlogro(pnx*pny),newlogs(pnx*pny),&
     &     addlogp(pnx),addlogq(pnx)
      
      integer ii,i,j,ip,iq
      real*8 q0,p0,t0,ro0,dlogp,dlogq

      character*80 ficin,ficout
!      data ficin/'data/EOSH/S_SESAME_alt_waterP-q.dat'/,
!     &     ficout/'data/EOSH/S_SESAME_smo_waterP-q.dat'/
      data ficin/'data/EOS_NN/H2O-REOS_logTPDSgrid.dat'/,               &
     &     ficout/'data/EOS_NN/H2O-REOS_logTPDSgrid_smo.dat'/


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     Reads the original file                                          c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      write(*,*)'On lit le fichier ',ficin
      open(unit=35,file=ficin,status='old',err=134)
      read(35,*) ny,nx
      write(*,*)'nx,ny:',nx,ny
      read(35,*) (logq(j),j=1,ny)
      read(35,*) (logp(i),i=1,nx)
      read(35,*) ((logro(i+nx*(j-1)),j=1,ny),i=1,nx)
      read(35,*) ((logs(i+nx*(j-1)),j=1,ny),i=1,nx)
 7    format(10f9.5)
 8    format(2i4)
      close(35)

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     Smooths table                                                     c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      do iq=1,ny
       newlogq(iq)=logq(iq)
       do ip=1,nx
        newlogp(ip)=logp(ip)
        newlogro(ip+nx*(iq-1))=logro(ip+nx*(iq-1))
        if ((iq.eq.1).or.(iq.eq.ny)) then
         if ((ip.eq.1).or.(ip.eq.nx)) then
          newlogs(ip+nx*(iq-1))=logs(ip+nx*(iq-1))
         else
          newlogs(ip+nx*(iq-1))=(2*logs(ip+nx*(iq-1))+                  &
     &         logs(ip-1+nx*(iq-1))+logs(ip+1+nx*(iq-1)))/4
         endif
        else
         if ((ip.eq.1).or.(ip.eq.nx)) then
          newlogs(ip+nx*(iq-1))=(2*logs(ip+nx*(iq-1))+                  &
     &         logs(ip+nx*(iq-2))+logs(ip+nx*(iq)))/4
         else
          newlogs(ip+nx*(iq-1))=(4*logs(ip+nx*(iq-1))+                  &
     &         logs(ip+nx*(iq-2))+logs(ip+nx*(iq))+                     &
     &         logs(ip-1+nx*(iq-1))+logs(ip+1+nx*(iq-1)))/8
         endif
        endif
        write(*,565)iq,ip,10**logq(iq),10**logp(ip),logro(ip+nx*(iq-1)),&
     &       newlogro(ip+nx*(iq-1)),logs(ip+nx*(iq-1)),                 &
     &       newlogs(ip+nx*(iq-1))
 565    format(i4,i4,'Q=',1p,d10.3,' P=',d10.3,' logRho=',d10.3,' => ', &
     &       d10.3,' logS=',d10.3," => ",d10.3)
       enddo
      enddo
      nnx=nx
      nny=ny

! Lets write the new table
      write(*,*)'On ecrit le fichier ',ficout
      open(unit=35,file=ficout,status='unknown')
      write(35,8) nny,nnx
      write(*,*)'nnx,nny:',nnx,nny
      write(35,7) (newlogq(j),j=1,nny)
      write(35,7) (newlogp(i),i=1,nnx)
      write(35,7) ((newlogro(i+nnx*(j-1)),j=1,nny),i=1,nnx)
      write(35,7) ((newlogs(i+nnx*(j-1)),j=1,nny),i=1,nnx)
      close(35)

! Tests
      i=1
      j=1
      do while ((i.gt.0).and.(j.gt.0))
       write(*,*)'Test: i (logP), j (logQ)?'
       read(*,*)i,j
       if ((i.gt.0).and.(j.gt.0)) then
! old table
        p0=10**(logp(i))
        q0=10**(logq(j))
        t0=10**((logq(j)+logp(i))/4)
        ro0=10**(logro(i+nnx*(j-1)))
        write(*,116)'old',p0,q0,t0,ro0,                                 &
     &       10**(logs(i+nnx*(j-1))),                                   &
     &       ro0*t0*8.3143d7/p0
! new table
        p0=10**(newlogp(i))
        q0=10**(newlogq(j))
        t0=10**((newlogq(j)+newlogp(i))/4)
        ro0=10**(newlogro(i+nnx*(j-1)))
        write(*,116)'new',p0,q0,t0,ro0,                                 &
     &       10**(newlogs(i+nnx*(j-1))),                                &
     &       ro0*t0*8.3143d7/p0
 116    format(1p,a3,"=> P=",d12.5," Q=",d12.5," T=",d12.5," Ro=",d12.5,&
     &       " S=",d12.5," mu=",d12.5)
       endif
      enddo

      return

 134  write(*,*)'Problem: file does not exist'

      end