! extend_EOS.f
! Fortran program to extend EOS tables calculated using SESAME into the ideal gas range
! Original data are from Saumon & Guillot, ApJ (2004)
! T. Guillot
!
! USE: compile with gfortran, run
! INPUT: files data/EOSH/S_SESAME_waterP-q.dat and S_SESAME_drysandP-q.dat
! OUPUT: files data/EOSH/S_SESAME_ext_waterP-q.dat and S_SESAME_ext_drysandP-q.dat
! No options, code itself is to be modified as needed.

      program extend_EOS
!     Version: 16-Oct-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 fic_water_in,fic_water_out,fic_drysand_in,           &
     &     fic_drysand_out,ficin,ficout
!$$$      data fic_water_in/'data/EOSH/S_SESAME_waterP-q.dat'/,
!$$$     &     fic_water_out/'data/EOSH/S_SESAME_ext_waterP-q.dat'/,
!$$$     &     fic_drysand_in/'data/EOSH/S_SESAME_drysandP-q.dat'/,
!$$$     &     fic_drysand_out/'data/EOSH/S_SESAME_ext_drysandP-q.dat'/
      data fic_water_in/'data/EOS_NN/H2O-REOS_logTPDSgrid.dat'/,        &
     &     fic_water_out/'data/EOS_NN/H2O-REOS_ext_logTPDSgrid.dat'/


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     Added pressures and Q values                                     c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!$$$      naddx=5
!$$$      dlogp=0.2d0
!$$$      do i=1,naddx
!$$$       addlogp(i)=5.d0+dlogp*(i-1)
!$$$      enddo
!$$$      naddy=5
!$$$      dlogq=0.2
!$$$      do i=1,naddy
!$$$       addlogq(i)=6.16d0+dlogq*i
!$$$      enddo
      naddx=0
      naddy=5
      dlogq=0.2d0
      do i=1,naddy
       addlogq(i)=2d0+dlogq*(i-1)
      enddo

      do ii=2,1,-1              !first sand, then water
       if (ii.eq.1) then
        ficin=fic_water_in
        ficout=fic_water_out
       else
        ficin=fic_drysand_in
        ficout=fic_drysand_out
       endif

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

       nnx=nx+naddx
       nny=ny+naddy

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     Shifts table                                                     c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       do iq=1,ny
        newlogq(iq)=logq(iq)
        do ip=naddx+1,nnx
         newlogp(ip)=logp(ip-naddx)
         newlogro(ip+nnx*(iq-1))=logro(ip-naddx+nx*(iq-1))
         newlogs(ip+nnx*(iq-1))=logs(ip-naddx+nx*(iq-1))
        enddo
       enddo

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     Adds new points                                                  c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! Adds the high-Q lines
       do iq=ny+1,nny
        newlogq(iq)=addlogq(iq-ny)
        q0=10**(logq(ny))
        do ip=naddx+1,nnx
         p0=10**(newlogp(ip))
         ro0=10**(newlogro(ip+nnx*(ny-1)))
         newlogro(ip+nnx*(iq-1))=newlogro(ip+nnx*(ny-1))-               &
     &        1/4.d0*(newlogq(iq)-newlogq(ny))
         newlogs(ip+nnx*(iq-1))=log10(10**(newlogs(ip+nnx*(ny-1)))+     &
     &        p0**(3/4.d0)/ro0/q0**(1/4.)*                              &
     &        (5/8.d0*(newlogq(iq)-newlogq(ny))))
        enddo
       enddo

! Now adds the low-pressure columns 
       do ip=1,naddx
        newlogp(ip)=addlogp(ip)
        p0=10**(logp(1))
        do iq=1,nny
         q0=10**(logq(ip))
         ro0=10**(logro(1+nx*(iq-1)))
         newlogro(ip+nnx*(iq-1))=newlogro(naddx+1+nnx*(iq-1))+          &
     &        3/4.d0*(newlogp(ip)-newlogp(naddx+1))
         newlogs(ip+nnx*(iq-1))=log10(10**(newlogs(naddx+1+nx*(iq-1)))+ &
     &        p0**(3/4.d0)/ro0/q0**(1/4.)*                              &
     &        (-3/8.d0*(newlogp(ip)-newlogp(naddx+1))))
        enddo
       enddo

! 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)

      enddo                     !water / sand loop

! 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
        p0=10**(newlogp(i))
        q0=10**(newlogq(j))
        t0=10**((newlogq(j)+newlogp(i))/4)
        ro0=10**(newlogro(i+nnx*(j-1)))
        write(*,116)p0,q0,t0,ro0,                                       &
     &       10**(newlogs(i+nnx*(j-1))),                                &
     &       ro0*t0*8.3143d7/p0
 116    format(1p,"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