program sss_start_128_112

! Purpose
! -------
! Interpolates climatological SSSs from the new, high-resolution Mk3L ocean
! model grid on to the grid used by the atmosphere model. This interpolation is
! performed in exactly the same fashion as within the coupled model.
!
! The values for consecutive months are then averaged in order to derive best
! estimates of the climatological values for the START of each month.
!
! Usage
! -----
! ./sss_start_128_112 <infile> <outfile>
!
! where:
!
! infile	file containing climatological SSSs from the OGCM
! outfile	output file
!
! History
! -------
! 2008 Feb 19	Steven Phipps	Original version

  implicit none

  include 'netcdf.inc'

  ! Define parameters
  integer, parameter :: imt = 130, &
                        jmt = 114, &
                        nx = 64, &
                        ny = 56
  real, parameter :: pi = 3.14159265, &
                     radian = 180.0/pi, &
                     radius = 6.37122e8, &
                     swldeg = -89.1331111924962

  ! Declare local variables
  character(len=80) :: infile, outfile
  integer :: fileid, i, io, j, jo, latdid, latvid, londid, lonvid, mondid, &
             month, monvid, sssid, sssstartid, status
  real :: sumdy
  real(kind=4), dimension(nx) :: lon
  real(kind=4), dimension(ny) :: lat
  real(kind=4), dimension(imt-2) :: lonts
  real(kind=4), dimension(jmt-2) :: latts
  real, dimension(jmt) :: cst, dyt, phi, phit
  real, dimension(nx, ny, 12) :: asal, sss_start
  real(kind=4), dimension(imt-2, jmt-2, 12) :: sss
  real, dimension(imt, jmt, 12) :: osal

  ! ---------------------
  ! BEGIN EXECUTABLE CODE
  ! ---------------------

  ! ----------------------------------------------------------------------
  ! The following section calculates the values of CST, the cosine of the
  ! latitude of each tracer gridbox, in exactly the same way as in ocean.f
  ! ----------------------------------------------------------------------

  ! Set the meridional dimension of each tracer gridbox in degrees
  dyt(  2) =  1.5717758068763
  dyt(  3) =  1.5717758068763
  dyt(  4) =  1.5873125395736
  dyt(  5) =  1.5873125395736
  dyt(  6) =  1.5904318409297
  dyt(  7) =  1.5904318409297
  dyt(  8) =  1.5915276112852
  dyt(  9) =  1.5915276112852
  dyt( 10) =  1.5920305756336
  dyt( 11) =  1.5920305756336
  dyt( 12) =  1.5923010877064
  dyt( 13) =  1.5923010877064
  dyt( 14) =  1.5924626726494
  dyt( 15) =  1.5924626726494
  dyt( 16) =  1.5925666534402
  dyt( 17) =  1.5925666534402
  dyt( 18) =  1.5926373800807
  dyt( 19) =  1.5926373800807
  dyt( 20) =  1.5926875876849
  dyt( 21) =  1.5926875876849
  dyt( 22) =  1.5927244533195
  dyt( 23) =  1.5927244533195
  dyt( 24) =  1.5927522680130
  dyt( 25) =  1.5927522680130
  dyt( 26) =  1.5927737232154
  dyt( 27) =  1.5927737232154
  dyt( 28) =  1.5927905752292
  dyt( 29) =  1.5927905752292
  dyt( 30) =  1.5928040087625
  dyt( 31) =  1.5928040087625
  dyt( 32) =  1.5928148456915
  dyt( 33) =  1.5928148456915
  dyt( 34) =  1.5928236699541
  dyt( 35) =  1.5928236699541
  dyt( 36) =  1.5928309049419
  dyt( 37) =  1.5928309049419
  dyt( 38) =  1.5928368629406
  dyt( 39) =  1.5928368629406
  dyt( 40) =  1.5928417775583
  dyt( 41) =  1.5928417775583
  dyt( 42) =  1.5928458254894
  dyt( 43) =  1.5928458254894
  dyt( 44) =  1.5928491414079
  dyt( 45) =  1.5928491414079
  dyt( 46) =  1.5928518283288
  dyt( 47) =  1.5928518283288
  dyt( 48) =  1.5928539649014
  dyt( 49) =  1.5928539649014
  dyt( 50) =  1.5928556105851
  dyt( 51) =  1.5928556105851
  dyt( 52) =  1.5928568093168
  dyt( 53) =  1.5928568093168
  dyt( 54) =  1.5928575920821
  dyt( 55) =  1.5928575920821
  dyt( 56) =  1.5928579786506
  dyt( 57) =  1.5928579786506
  dyt( 58) =  1.5928579786506
  dyt( 59) =  1.5928579786506
  dyt( 60) =  1.5928575920821
  dyt( 61) =  1.5928575920821
  dyt( 62) =  1.5928568093168
  dyt( 63) =  1.5928568093168
  dyt( 64) =  1.5928556105851
  dyt( 65) =  1.5928556105851
  dyt( 66) =  1.5928539649014
  dyt( 67) =  1.5928539649014
  dyt( 68) =  1.5928518283288
  dyt( 69) =  1.5928518283288
  dyt( 70) =  1.5928491414079
  dyt( 71) =  1.5928491414079
  dyt( 72) =  1.5928458254894
  dyt( 73) =  1.5928458254894
  dyt( 74) =  1.5928417775583
  dyt( 75) =  1.5928417775583
  dyt( 76) =  1.5928368629406
  dyt( 77) =  1.5928368629406
  dyt( 78) =  1.5928309049419
  dyt( 79) =  1.5928309049419
  dyt( 80) =  1.5928236699541
  dyt( 81) =  1.5928236699541
  dyt( 82) =  1.5928148456915
  dyt( 83) =  1.5928148456915
  dyt( 84) =  1.5928040087625
  dyt( 85) =  1.5928040087625
  dyt( 86) =  1.5927905752292
  dyt( 87) =  1.5927905752292
  dyt( 88) =  1.5927737232154
  dyt( 89) =  1.5927737232154
  dyt( 90) =  1.5927522680130
  dyt( 91) =  1.5927522680130
  dyt( 92) =  1.5927244533195
  dyt( 93) =  1.5927244533195
  dyt( 94) =  1.5926875876849
  dyt( 95) =  1.5926875876849
  dyt( 96) =  1.5926373800807
  dyt( 97) =  1.5926373800807
  dyt( 98) =  1.5925666534402
  dyt( 99) =  1.5925666534402
  dyt(100) =  1.5924626726494
  dyt(101) =  1.5924626726494
  dyt(102) =  1.5923010877064
  dyt(103) =  1.5923010877064
  dyt(104) =  1.5920305756336
  dyt(105) =  1.5920305756336
  dyt(106) =  1.5915276112852
  dyt(107) =  1.5915276112852
  dyt(108) =  1.5904318409297
  dyt(109) =  1.5904318409297
  dyt(110) =  1.5873125395736
  dyt(111) =  1.5873125395736
  dyt(112) =  1.5717758068763
  dyt(113) =  1.5717758068763

  ! Convert DYT to centimetres
  dyt = dyt * radius / radian

  ! Calculate the cosine of the latitude of each tracer gridbox
  phi(1) = swldeg / radian
  phit(1) = phi(1) - 0.5 * dyt(1) / radius
  sumdy = phi(1)
  do j = 1, jmt
    if (j /= jmt) then
      sumdy = sumdy + dyt(j+1) / radius
      phi(j+1) = sumdy
    end if
    if (j /= 1) phit(j) = 0.5 * (phi(j-1) + phi(j))
    cst(j) = cos(phit(j))
  end do

  ! -----------------------------------------
  ! Get the climatological SSSs from the OGCM
  ! -----------------------------------------

  ! Get the name of the input file
  call getarg(1, infile)

  ! Get the climatological SSSs from the OGCM
  status = nf_open(infile, nf_nowrite, fileid)
  status = nf_inq_varid(fileid, "longitude", lonvid)
  status = nf_inq_varid(fileid, "latitude", latvid)
  status = nf_inq_varid(fileid, "sss", sssid)
  status = nf_get_var_real(fileid, lonvid, lonts)
  status = nf_get_var_real(fileid, latvid, latts)
  status = nf_get_var_real(fileid, sssid, sss)
  status = nf_close(fileid)

  ! Set the salinities over land to 35 psu, as in the coupled model
  do month = 1, 12
    do j = 1, jmt-2
      do i = 1, imt-2
        if (sss(i, j, month) < -9000.0) sss(i, j, month) = 35.0
      end do
    end do
  end do

  ! Transfer the data to the array OSAL
  osal(2:imt-1, 2:jmt-1, :) = sss

  ! ------------------------------------------------------------------------
  ! The following section of code obtains the SSS for each atmosphere model
  ! gridbox by calculating the area-weighted average of the SSS for each of
  ! the four corresponding OGCM gridboxes.
  !
  ! These calculations are performed in exactly the same way as in ocinit.f.
  ! ------------------------------------------------------------------------

  do month = 1, 12
    do j = 1, ny
      jo = 2 * j
      do i = 1, nx
        io = 2 * i
        asal(i, j, month) = (cst(jo) * osal(io, jo, month) + &
                             cst(jo) * osal(io+1, jo, month) + &
                             cst(jo+1) * osal(io, jo+1, month) + &
                             cst(jo+1) * osal(io+1, jo+1, month)) / &
                            (2.0 * (cst(jo) + cst(jo+1)))
      end do
    end do
  end do

  ! --------------------------------------------------------------------
  ! Average consecutive monthly values in order to obtain best estimates
  ! of the climatological values for the start of each month
  ! --------------------------------------------------------------------

  do month = 1, 12
    if (month == 1) then
      sss_start(:, :, 1) = 0.5 * (asal(:, :, 12) + asal(:, :, 1))
    else
      sss_start(:, :, month) = 0.5 * (asal(:, :, month-1) + asal(:, :, month))
    end if
  end do

  ! ----------------------------------------------
  ! Derive the axis data for the interpolated SSSs
  ! ----------------------------------------------

  do i = 1, nx
    io = 2 * i
    lon(i) = 0.5 * (lonts(io-1) + lonts(io))
  end do
  do j = 1, ny
    jo = 2 * j
    lat(j) = 0.5 * (latts(jo-1) + latts(jo))
  end do

  ! --------------------------------
  ! Save the data to the output file
  ! --------------------------------
  
  ! Get the name of the output file
  call getarg(2, outfile)

  ! (1) Create netCDF file
  status = nf_create(outfile, nf_noclobber, fileid)

  ! (2) Define dimensions
  status = nf_def_dim(fileid, "longitude", nx, londid)
  status = nf_def_dim(fileid, "latitude", ny, latdid)
  status = nf_def_dim(fileid, "month", 12, mondid)

  ! (3) Define variables
  status = nf_def_var(fileid, "longitude", nf_float, 1, londid, lonvid)
  status = nf_def_var(fileid, "latitude", nf_float, 1, latdid, latvid)
  status = nf_def_var(fileid, "month", nf_int, 1, mondid, monvid)
  status = nf_def_var(fileid, "sss_clim", nf_double, 3, &
                      (/ londid, latdid, mondid /), sssid)
  status = nf_def_var(fileid, "sss", nf_double, 3, &
                      (/ londid, latdid, mondid /), sssstartid)

  ! (4) Put attributes
  status = nf_put_att_text(fileid, sssid, "long_name", 20, &
                           "Sea surface salinity")
  status = nf_put_att_text(fileid, sssstartid, "long_name", 20, &
                           "Sea surface salinity")
  status = nf_put_att_text(fileid, lonvid, "units", 12, "degrees_east")
  status = nf_put_att_text(fileid, latvid, "units", 13, "degrees_north")
  status = nf_put_att_text(fileid, monvid, "units", 6, "months")
  status = nf_put_att_text(fileid, sssid, "units", 3, "psu")
  status = nf_put_att_text(fileid, sssstartid, "units", 3, "psu")
  status = nf_put_att_text(fileid, lonvid, "point_spacing", 4, "even")
  status = nf_put_att_text(fileid, latvid, "point_spacing", 6, "uneven")
  status = nf_put_att_text(fileid, monvid, "point_spacing", 4, "even")
  status = nf_put_att_text(fileid, lonvid, "modulo", 1, " ")
  status = nf_put_att_text(fileid, monvid, "modulo", 1, " ")

  ! (5) Exit define mode
  status = nf_enddef(fileid)

  ! (6) Write data
  status = nf_put_var_real(fileid, lonvid, lon)
  status = nf_put_var_real(fileid, latvid, lat)
  status = nf_put_var_int(fileid, monvid, (/ 1, 2, 3, 4, 5, 6, &
                                             7, 8, 9, 10, 11, 12 /))
  status = nf_put_var_double(fileid, sssid, asal)
  status = nf_put_var_double(fileid, sssstartid, sss_start)

  ! (7) Close output file
  status = nf_close(fileid)

end program sss_start_128_112