program comp_rhop
!
! =============================================================================
! 
!  This program computes the potential density rhop based on the equation 
!       used in NEMO. Beware that the configuration is hard-coded. LX, LY
!       LZ, LT need to be changed for each new configuration and the program
!       needs to be recompiled with the makefile provided.
!
!  Usage ./comp_rhop with a namelist called namelist_rhop:
!              &density 
!                input_file = ''
!                input_Tvar = ''
!                input_Svar = ''
!                input_Dvar = ''
!                outfile    = ''
!              /
! 
!  History : Virginie Guemas - Initial version PhD tools       - February 2008
!                            - Namelist and remove hard-coding - April    2014
!  
! =============================================================================
!
      USE eosbn2        ! Routine from NEMO to compute rhop
      IMPLICIT NONE
      INCLUDE 'netcdf.inc'

      CHARACTER(80) :: &
        & input_file, & ! Input netcdf file name
        & input_Tvar, & ! Input netcdf potential temperature variable name
        & input_Svar, & ! Input netcdf salinity variable name
        & input_Dvar, & ! Input netcdf depth variable name
        & outfile       ! Ouput netcdf file name

      NAMELIST/density/input_file,input_Tvar,input_Svar,input_Dvar,outfile
      
      INTEGER,PARAMETER :: inam=1
      INTEGER :: ncid, varid, istat
      INTEGER,DIMENSION(4) :: dimoutids
      INTEGER,PARAMETER :: LX=1442,LY=1021,LZ=46,LT=1
        ! These dimensions need to be changed also in eosbn2.f90 for each new
        ! configuration

      DOUBLE PRECISION,DIMENSION(LX,LY,LZ,LT)   :: PT,PS
      DOUBLE PRECISION,DIMENSION(LZ)            :: PD
      DOUBLE PRECISION,DIMENSION(LX,LY,LZ,LT)   :: PRHO
      DOUBLE PRECISION,DIMENSION(LX,LY,LZ)      :: MASK
      DOUBLE PRECISION,DIMENSION(LX,LY,LZ,LT)   :: Zrho_smow,PO
      INTEGER :: JI,JJ,JK,JT  
      
!
! =============================================================================
!
! Read the namelist
! ==================

      OPEN(UNIT=inam,FILE='namelist_rhop',FORM='FORMATTED',ACCESS=&
           &'SEQUENTIAL',STATUS='OLD',IOSTAT=istat)
      IF (istat /=0 ) THEN
        PRINT*,'IOSTAT = ',istat
        STOP"Problem opening namelist_rhop"
      END IF

      REWIND(inam)
      READ(UNIT=inam,NML=density)


! Read input thermodynamic fields
! =================================
!
      istat = NF_OPEN (input_file,0,ncid )
      call check_stat(istat,'Opening input file')

      istat = NF_INQ_VARID(ncid,input_Tvar,varid)
      call check_stat(istat,'Finding input temperature variable')

      istat = NF_GET_VAR_DOUBLE(ncid,varid,PT)
      call check_stat(istat,'Reading input temperature')

      istat = NF_INQ_VARID(ncid,input_Svar,varid)
      call check_stat(istat,'Finding input salinity variable')

      istat = NF_GET_VAR_DOUBLE(ncid,varid,PS)
      call check_stat(istat,'Reading input temperature')

      istat = NF_INQ_VARID(ncid,input_Dvar,varid)
      call check_stat(istat,'Finding input depth variable')

      istat = NF_GET_VAR_DOUBLE(ncid,varid,PD)
      call check_stat(istat,'Reading input depth')

      istat = NF_CLOSE ( ncid )
      call check_stat(istat,'Closing input file')
!
! compute potential density
!   
      MASK = 1.
      WHERE(PT(:,:,:,1) > 1e19 .or. PT(:,:,:,1)< -9e30)
        MASK = 0.
      END WHERE
     
      DO JT = 1,LT
        call eos_insitu_pot( PT(:,:,:,JT), PS(:,:,:,JT), PD, MASK, &
        &PRHO(:,:,:,JT) )
        WHERE (MASK < 0.5)
          PRHO(:,:,:,JT) = 1e20
        END WHERE
      END DO
!
! Write outputs
! ==============
!
      istat = NF_CREATE(outfile,NF_NOCLOBBER,ncid)
      call check_stat(istat,'Opening output file')

      istat = NF_DEF_DIM(ncid,'x',LX,dimoutids(1))
      call check_stat(istat,'Defining x dimension in output file')

      istat = NF_DEF_DIM(ncid,'y',LY,dimoutids(2))
      call check_stat(istat,'Defining y dimension in output file')

      istat = NF_DEF_DIM(ncid,'z',LZ,dimoutids(3))
      call check_stat(istat,'Defining z dimension in output file')

      istat = NF_DEF_DIM(ncid,'t',LT,dimoutids(4))
      call check_stat(istat,'Defining t dimension in output file')

      istat = NF_DEF_VAR(ncid,'rhop',NF_DOUBLE,4,dimoutids,varid)
      call check_stat(istat,'Defining rhop variable in output file')

      istat = NF_PUT_ATT_DOUBLE(ncid,varid,'missing_value',NF_DOUBLE,1,DBLE(1e20))
      call check_stat(istat,'Defining missing value for rhop')

      istat = NF_ENDDEF(ncid)
      call check_stat(istat,'Closing output file definition')

      istat = NF_PUT_VAR_DOUBLE(ncid,varid,PRHO)
      call check_stat(istat,'Writing rhop in output file')

      istat = NF_CLOSE ( ncid )
      call check_stat(istat,'Closing output file')

end program comp_rhop