program generate_sssa

! Purpose
! -------
! Generates the auxiliary file sssa.nc for the Mk3L atmosphere model.
!
! This file contains climatological sea surface salinities [SSS].
!
! Note that this program assumes that the climatological data is provided at
! 64-bit precision.
!
! Usage
! -----
! ./generate_sssa <sss_file> <sss_var> <lon_var> <lat_var>
!
! where:
!
! sss_file	file containing climatological SSSs
! sss_var	variable containing climatological SSSs
! lon_var	variable containing longitude data for the AGCM
! lat_var	variable containing latitude data for the AGCM
!
! History
! -------
! 2008 Feb 27	Steven Phipps	Original version

  implicit none

  include 'netcdf.inc'

  ! Define local parameters
  integer, parameter :: nx = 64, &
                        ny = 56

  ! Declare local variables
  character(len=80) :: lon_var, lat_var, sss_file, sss_var, title
  character(len=160) :: history
  integer :: datid, latdid, latvid, londid, lonvid, mondid, monvid, ncid, &
             status
  real(kind=4), dimension(nx) :: lon
  real(kind=4), dimension(ny) :: lat
  real, dimension(nx, ny, 12) :: sss

  ! Get command-line arguments
  call getarg(1, sss_file)
  call getarg(2, sss_var)
  call getarg(3, lon_var)
  call getarg(4, lat_var)

  ! Get the climatological SSSs and the grid data
  status = nf_open(trim(sss_file), nf_nowrite, ncid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_inq_varid(ncid, trim(lon_var), lonvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_inq_varid(ncid, trim(lat_var), latvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_inq_varid(ncid, trim(sss_var), datid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_get_var_real(ncid, lonvid, lon)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_get_var_real(ncid, latvid, lat)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_get_var_double(ncid, datid, sss)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_close(ncid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! Convert the SSSs from psu to kg/kg
  sss = 0.001 * sss

  ! Get title and history data for the output file
  write (*, *)
  write (*, *) "Please enter the title and history for the file sssa.nc :"
  write (*, *)
  write (*, '(a)', advance='no') "Title   : "
  read (*, '(a)') title
  write (*, '(a)', advance='no') "History : "
  read (*, '(a)') history

  ! Generate the auxiliary file sssa.nc
  !
  ! (1) Create netCDF file
  status = nf_create("sssa.nc", nf_noclobber, ncid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! (2) Define dimensions
  status = nf_def_dim(ncid, "longitude", nx, londid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_dim(ncid, "latitude", ny, latdid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_dim(ncid, "month", 12, mondid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! (3) Define variables
  status = nf_def_var(ncid, "longitude", nf_float, 1, londid, lonvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(ncid, "latitude", nf_float, 1, latdid, latvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(ncid, "month", nf_int, 1, mondid, monvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(ncid, "sss", nf_double, 3, &
                      (/ londid, latdid, mondid /), datid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! (4) Create attributes
  status = nf_put_att_text(ncid, lonvid, "units", 12, "degrees_east")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, latvid, "units", 13, "degrees_north")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, monvid, "units", 6, "months")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, lonvid, "modulo", 1, " ")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, monvid, "modulo", 1, " ")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, lonvid, "point_spacing", 4, "even")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, latvid, "point_spacing", 6, "uneven")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, monvid, "point_spacing", 4, "even")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, nf_global, "title", len(trim(title)), &
                           trim(title))
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(ncid, nf_global, "history", len(trim(history)), &
                           trim(history))
  if (status /= nf_noerr) stop "***  netCDF error"

  ! (5) Exit define mode
  status = nf_enddef(ncid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! (6) Write data
  status = nf_put_var_real(ncid, lonvid, lon)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_var_real(ncid, latvid, lat)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_var_int(ncid, monvid, (/ 1, 2, 3, 4, 5, 6, &
                                           7, 8, 9, 10, 11, 12 /))
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_var_double(ncid, datid, sss)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! (7) Close netCDF file
  status = nf_close(ncid)
  if (status /= nf_noerr) stop "***  netCDF error"

end program generate_sssa