program currents_start_128_112

! Purpose
! -------
! Interpolates climatological surface currents from the new, high-resolution
! Mk3L ocean model grid on to the grid used by the sea ice 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
! -----
! ./currents_start_128_112 <infile> <outfile>
!
! where:
!
! infile	file containing climatological surface currents from the OGCM
! outfile	output file
!
! History
! -------
! 2008 Jan 27	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, iom1, iop1, j, jo, jom1, jop1, latuvdid, &
             latuvvid, lonuvdid, lonuvvid, mondid, month, monvid, ocuid, &
             ocustartid, ocvid, ocvstartid, status
  real :: sumdy
  real(kind=4), dimension(imt-2) :: lonuv
  real(kind=4), dimension(jmt-2) :: latuv
  real, dimension(jmt) :: cs, dyt, phi
  real, dimension(nx+2, ny+2, 12) :: ocu, ocu_start, ocv, ocv_start
  real(kind=4), dimension(imt-2, jmt-2, 12) :: uco_in, vco_in
  real, dimension(imt, jmt, 12) :: uco, vco

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

  ! ------------------------------------------------------------------------
  ! The following section calculates the values of CS, the cosine of the
  ! latitude of each velocity 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 velocity gridbox
  phi(1) = swldeg / radian
  sumdy = phi(1)
  do j = 1, jmt
    if (j /= jmt) then
      sumdy = sumdy + dyt(j+1) / radius
      phi(j+1) = sumdy
    end if
    cs(j) = cos(phi(j))
  end do

  ! -----------------------------------------------------
  ! Get the climatological surface currents from the OGCM
  ! -----------------------------------------------------

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

  ! Get the climatological surface currents from the OGCM
  status = nf_open(infile, nf_nowrite, fileid)
  status = nf_inq_varid(fileid, "longitude", lonuvvid)
  status = nf_inq_varid(fileid, "latitude", latuvvid)
  status = nf_inq_varid(fileid, "ocu", ocuid)
  status = nf_inq_varid(fileid, "ocv", ocvid)
  status = nf_get_var_real(fileid, lonuvvid, lonuv)
  status = nf_get_var_real(fileid, latuvvid, latuv)
  status = nf_get_var_real(fileid, ocuid, uco_in)
  status = nf_get_var_real(fileid, ocvid, vco_in)
  status = nf_close(fileid)

  ! Set the currents over land to zero, as in the coupled model
  do month = 1, 12
    do j = 1, jmt-2
      do i = 1, imt-2
        if (uco_in(i, j, month) < -9000.0) uco_in(i, j, month) = 0.0
        if (vco_in(i, j, month) < -9000.0) vco_in(i, j, month) = 0.0
      end do
    end do
  end do

  ! Transfer the data to the arrays UCO and VCO and pad in the x-direction,
  ! as in the coupled model
  uco(2:imt-1, 2:jmt-1, :) = uco_in
  vco(2:imt-1, 2:jmt-1, :) = vco_in
  uco(1, :, :) = uco(imt-1, :, :)
  vco(1, :, :) = vco(imt-1, :, :)
  uco(imt, :, :) = uco(2, :, :)
  vco(imt, :, :) = vco(2, :, :)

  ! Convert the currents to cm/s
  uco = uco * 100.0
  vco = vco * 100.0

  ! -------------------------------------------------------------------------
  ! The following section of code obtains the surface velocities for each ice
  ! model gridbox by calculating the area-weighted average of the surface
  ! velocities for the corresponding region of the OGCM grid.
  !
  ! These calculations are performed in exactly the same way as in ocicurr.f.
  ! -------------------------------------------------------------------------

  do j = 2, ny+1
    jo = 2 * j - 1
    jom1 = jo - 1
    jop1 = jo + 1

    do i = 1, nx+1
      io = 2 * i - 1
      iom1 = io - 1
      iop1 = io + 1
      if (iom1 .eq. 0) iom1 = imt - 2
          
      if (j .eq. ny+1) then

        ocu(i, j, :) = (2.0 * cs(jom1) * uco(io, jom1, :) + &
                        cs(jom1) * uco(iop1, jom1, :) + &
                        cs(jom1) * uco(iom1, jom1, :)) / &
                       (4.0 * cs(jom1) + 8.0 * cs(jo) + &
                        4.0 * cs(jop1))
        ocv(i, j, :) = (2.0 * cs(jom1) * vco(io, jom1, :) + &
                        cs(jom1) * vco(iop1, jom1, :) + &
                        cs(jom1) * vco(iom1, jom1, :)) / &
                       (4.0 * cs(jom1) + 8.0 * cs(jo) + &
                        4.0 * cs(jop1))

      else

        ocu(i, j, :) = (4.0 * cs(jo) * uco(io, jo, :) + &
                        2.0 * cs(jop1) * uco(io, jop1, :) + &
                        2.0 * cs(jo) * uco(iop1, jo, :) + &
                        2.0 * cs(jom1) * uco(io, jom1, :) + &
                        2.0 * cs(jo) * uco(iom1, jo, :) + &
                        cs(jop1) * uco(iop1, jop1, :) + &
                        cs(jom1) * uco(iop1, jom1, :) + &
                        cs(jom1) * uco(iom1, jom1, :) + &
                        cs(jop1) * uco(iom1, jop1, :)) / &
                       (4.0 * cs(jom1) + 8.0 * cs(jo) + &
                        4.0 * cs(jop1))
        ocv(i, j, :) = (4.0 * cs(jo) * vco(io, jo, :) + &
                        2.0 * cs(jop1) * vco(io, jop1, :) + &
                        2.0 * cs(jo) * vco(iop1, jo, :) + &
                        2.0 * cs(jom1) * vco(io, jom1, :) + &
                        2.0 * cs(jo) * vco(iom1, jo, :) + &
                        cs(jop1) * vco(iop1, jop1, :) + &
                        cs(jom1) * vco(iop1, jom1, :) + &
                        cs(jom1) * vco(iom1, jom1, :) + &
                        cs(jop1) * vco(iom1, jop1, :)) / &
                       (4.0 * cs(jom1) + 8.0 * cs(jo) + &
                        4.0 * cs(jop1))
      end if

    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
      ocu_start(:, :, 1) = 0.5 * (ocu(:, :, 12) + ocu(:, :, 1))
      ocv_start(:, :, 1) = 0.5 * (ocv(:, :, 12) + ocv(:, :, 1))
    else
      ocu_start(:, :, month) = 0.5 * (ocu(:, :, month-1) + ocu(:, :, month))
      ocv_start(:, :, month) = 0.5 * (ocv(:, :, month-1) + ocv(:, :, month))
    end if
  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, lonuvdid)
  status = nf_def_dim(fileid, "latitude", ny, latuvdid)
  status = nf_def_dim(fileid, "month", 12, mondid)

  ! (3) Define variables
  status = nf_def_var(fileid, "longitude", nf_float, 1, lonuvdid, lonuvvid)
  status = nf_def_var(fileid, "latitude", nf_float, 1, latuvdid, latuvvid)
  status = nf_def_var(fileid, "month", nf_int, 1, mondid, monvid)
  status = nf_def_var(fileid, "ocu_clim", nf_double, 3, &
                      (/ lonuvdid, latuvdid, mondid /), ocuid)
  status = nf_def_var(fileid, "ocv_clim", nf_double, 3, &
                      (/ lonuvdid, latuvdid, mondid /), ocvid)
  status = nf_def_var(fileid, "ocu", nf_double, 3, &
                      (/ lonuvdid, latuvdid, mondid /), ocustartid)
  status = nf_def_var(fileid, "ocv", nf_double, 3, &
                      (/ lonuvdid, latuvdid, mondid /), ocvstartid)

  ! (4) Put attributes
  status = nf_put_att_text(fileid, lonuvvid, "long_name", 26, &
                           "Longitude of UV gridpoints")
  status = nf_put_att_text(fileid, latuvvid, "long_name", 25, &
                           "Latitude of UV gridpoints")
  status = nf_put_att_text(fileid, ocuid, "long_name", 26, &
                           "Zonal sea surface velocity")
  status = nf_put_att_text(fileid, ocvid, "long_name", 31, &
                           "Meridional sea surface velocity")
  status = nf_put_att_text(fileid, ocustartid, "long_name", 26, &
                           "Zonal sea surface velocity")
  status = nf_put_att_text(fileid, ocvstartid, "long_name", 31, &
                           "Meridional sea surface velocity")
  status = nf_put_att_text(fileid, lonuvvid, "units", 12, "degrees_east")
  status = nf_put_att_text(fileid, latuvvid, "units", 13, "degrees_north")
  status = nf_put_att_text(fileid, monvid, "units", 6, "months")
  status = nf_put_att_text(fileid, ocuid, "units", 4, "cm/s")
  status = nf_put_att_text(fileid, ocvid, "units", 4, "cm/s")
  status = nf_put_att_text(fileid, ocustartid, "units", 4, "cm/s")
  status = nf_put_att_text(fileid, ocvstartid, "units", 4, "cm/s")
  status = nf_put_att_text(fileid, lonuvvid, "point_spacing", 4, "even")
  status = nf_put_att_text(fileid, latuvvid, "point_spacing", 6, "uneven")
  status = nf_put_att_text(fileid, monvid, "point_spacing", 4, "even")
  status = nf_put_att_text(fileid, lonuvvid, "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, lonuvvid, lonuv(2:imt-2:2))
  status = nf_put_var_real(fileid, latuvvid, latuv(2:jmt-2:2))
  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, ocuid, ocu(2:nx+1, 2:ny+1, :))
  status = nf_put_var_double(fileid, ocvid, ocv(2:nx+1, 2:ny+1, :))
  status = nf_put_var_double(fileid, ocustartid, ocu_start(2:nx+1, 2:ny+1, :))
  status = nf_put_var_double(fileid, ocvstartid, ocv_start(2:nx+1, 2:ny+1, :))

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

end program currents_start_128_112