program ts_stats_128_112_21

! Purpose
! -------
! Reads in annual-mean output from the Mk3L OGCM, and produces a single netCDF
! file containing the following temperature and salinity statistics:
!
! STFHT		mean surface heat flux [W/m^2]
! STFSAL        mean surface salinity tendency [s^-1]
!
! ISTFHT	integrated surface heat flux [J]
! ISTFSAL	integrated surface salinity tendency [kg/kg]
!
! TEMPYZ	zonal-mean potential temperature [degC]
! SALYZ		zonal-mean salinity [psu]
!
! TEMPZ		zonal- and meridional-mean potential temperature [degC]
! SALZ		zonal- and meridional-mean salinity [psu]
!
! TEMP		mean potential temperature [degC]
! SAL		mean salinity [psu]
!
! TEMP1		mean potential temperature of the upper ocean [0-800m, degC]
! TEMP2		mean potential temperature of the mid-ocean [800-2350m, degC]
! TEMP3		mean potential temperature of the deep ocean [2350-4600m, degC]
!
! SAL1		mean salinity of the upper ocean [0-800m, psu]
! SAL2		mean salinity of the mid-ocean [800-2350m, psu]
! SAL3		mean salinity of the deep ocean [2350-4600m, psu]
!
! DTE		equivalent change in mean potential temperature [K]
! DSE		equivalent change in mean salinity [psu]
!
! TERR		error in mean potential temperature [K]
! SERR		error in mean salinity [psu]
!
! Note that, in order to give accurate results, this program must be compiled
! using the same compiler options that were used to compile the model.
! Specifically, it must be compiled with a default real precision of 64 bits.
!
! History
! -------
! 2008 Dec 15	Steven Phipps	Original Fortran 90 version, based on the
!				original implementation in IDL
! 2009 Apr 22	Steven Phipps	Modified for five-character experiment names
! 2009 Aug 28	Steven Phipps	Correct command which writes title attribute

  implicit none

  include 'netcdf.inc'

  ! Define parameters
  integer, parameter :: imt = 130, &
                        jmt = 114, &
                        nx = 128, &
                        ny = 112, &
                        nz = 21, &
                        years_per_file = 50
  real(kind=4), parameter :: missing_value4 = -1.0e34_4
  real, parameter :: missing_value8 = -1.0e34, &
                     pi = 3.14159265, &
                     radian = 180.0/pi, &
                     radius = 6.37122e8, &
                     secs_year = 3.1536e7, &
                     swldeg = -89.1331111924962, &
                     rhocw = 4.1e6
  real, dimension(nz), parameter :: dz = (/ 25.0, 25.0, 30.0, 37.0, 43.0, &
                                            50.0, 60.0, 80.0, 120.0, 150.0, &
                                            180.0, 210.0, 240.0, 290.0, &
                                            360.0, 450.0, 450.0, 450.0, &
                                            450.0, 450.0, 450.0 /)

  ! Declare local variables
  character(len=5) :: run_name, year1s, year2s, yr1s, yr2s, zone
  character(len=8) :: date
  character(len=10) :: time
  character(len=23) :: outfile
  character(len=28) :: infile
  character(len=36) :: title
  character(len=64) :: history
  integer :: datid, depdid, depvid, dseid, dteid, i, inid, istfhtid, &
             istfsalid, iyear, j, k, l, latdid, latvid, nvalues, nyears, &
             outid, salid, sal1id, sal2id, sal3id, salyzid, salzid, serrid, &
             stfhtid, stfsalid, status, tempid, temp1id, temp2id, temp3id, &
             tempyzid, tempzid, terrid, timdid, timvid, year, year1, year2, &
             yr1, yr2
  integer, dimension(:), allocatable :: years
  integer, dimension(8) :: date_time
  logical, dimension(nx, ny, nz) :: mask
  real(kind=4), dimension(nz) :: zts
  real(kind=4), dimension(nx, ny, years_per_file) :: stfht_in, stfsal_in
  real(kind=4), dimension(nx, ny, nz, years_per_file) :: sal_in, temp_in
  real :: area, dse_old, dte_old, sal_old, serr_old, stfht_old, stfsal_old, &
          sumdy, temp_old, terr_old, volume, vtotal
  real, dimension(:), allocatable :: istfht, istfsal
  real, dimension(imt) :: dxt
  real, dimension(jmt) :: cst, dyt, phi, phit
  real, dimension(years_per_file) :: dse, dte, sal, sal1, sal2, sal3, serr, &
                                     stfht, stfsal, temp, temp1, temp2, &
                                     temp3, terr
  real, dimension(nx, ny) :: dats
  real, dimension(nz, years_per_file) :: salz, tempz
  real, dimension(nx, ny, nz) :: dvts
  real, dimension(ny, nz, years_per_file) :: salyz, tempyz

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

  ! ------------------------------------------------------------------
  ! The following section calculates the values of CST, DXT and DYT in
  ! exactly the same way as within the model
  ! ------------------------------------------------------------------

  ! 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

  ! ---------------------------------------------------------------------
  ! Calculate DXT, the x-dimension of each gridbox, exactly as in ocean.f
  ! ---------------------------------------------------------------------

  ! Set the x-dimension of each box in degrees, and convert to centimetres
  do i = 2, imt-1
    dxt(i) = 360.0 / float(imt-2)
    dxt(i) = dxt(i) * radius / radian
  end do

  ! Set cyclic conditions on DXT
  dxt(1) = dxt(imt-1)
  dxt(imt) = dxt(2)

  ! --------------------------------------------
  ! Now begin processing the output of the model
  ! --------------------------------------------

  ! Get the command-line arguments
  call getarg(1, run_name)
  call getarg(2, year1s)
  call getarg(3, year2s)
  read (year1s, *) year1
  read (year2s, *) year2

  ! Calculate the number of years of data to read
  nyears = year2 - year1 + 1

  ! -------------------------
  ! Initialise various arrays
  ! -------------------------

  ! Allocate output arrays
  allocate (istfht(nyears), istfsal(nyears), years(nyears))

  ! Calculate the surface area of each gridbox in m^2, remembering that DXT and
  ! DYT are in cm (as within the model)
  do j = 1, ny
    do i = 1, nx
      dats(i, j) = 1.0e-4 * cst(j+1) * dxt(i+1) * dyt(j+1)
    end do
  end do

  ! Calculate the volume of each gridbox in m^3, remembering that DXT and DYT
  ! are in cm (as within the model) and that DZ is in m
  do k = 1, nz
    do j = 1, ny
      do i = 1, nx
        dvts(i, j, k) = 1.0e-4 * cst(j+1) * dxt(i+1) * dyt(j+1) * dz(k)
      end do
    end do
  end do

  ! Generate depth axis data for the output file
  zts(1) = real(0.5*dz(1), 4)
  do k = 2, nz
    zts(k) = zts(k-1) + real(0.5*(dz(k-1) + dz(k)), 4)
  end do

  ! Generate time axis data for the output file
  do i = 1, nyears
    years(i) = year1 + i - 1
  end do

  ! Get current date and time
  call date_and_time(date, time, zone, date_time)

  ! ----------------------
  ! Create the output file
  ! ----------------------

  ! Derive name of output file
  write (year1s, '(i5)') year1
  write (year2s, '(i5)') year2
  if (year1 < 10000) write (year1s, '(a1, i4)'), "0", year1
  if (year1 < 1000) write (year1s, '(a2, i3)'), "00", year1
  if (year1 < 100) write (year1s, '(a3, i2)'), "000", year1
  if (year1 < 10) write (year1s, '(a4, i1)'), "0000", year1
  if (year2 < 10000) write (year2s, '(a1, i4)'), "0", year2
  if (year2 < 1000) write (year2s, '(a2, i3)'), "00", year2
  if (year2 < 100) write (year2s, '(a3, i2)'), "000", year2
  if (year2 < 10) write (year2s, '(a4, i1)'), "0000", year2
  outfile = "ts." // run_name // "." // year1s // "-" // year2s // ".nc"

  ! Create netCDF file
  status = nf_create(trim(outfile), nf_noclobber, outid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! Define dimensions
  status = nf_def_dim(outid, "latitude", ny, latdid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_dim(outid, "depth", nz, depdid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_dim(outid, "year", nyears, timdid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! Define variables
  status = nf_def_var(outid, "latitude", nf_float, 1, latdid, latvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "depth", nf_float, 1, depdid, depvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "year", nf_int, 1, timdid, timvid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "stfht", nf_float, 1, timdid, stfhtid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "stfsal", nf_float, 1, timdid, stfsalid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "istfht", nf_float, 1, timdid, istfhtid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "istfsal", nf_float, 1, timdid, istfsalid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "tempyz", nf_float, 3, &
                      (/ latdid, depdid, timdid /), tempyzid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "salyz", nf_float, 3, &
                      (/ latdid, depdid, timdid /), salyzid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "tempz", nf_float, 2, (/ depdid, timdid /), &
                      tempzid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "salz", nf_float, 2, (/ depdid, timdid /), &
                      salzid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "temp", nf_float, 1, timdid, tempid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "sal", nf_float, 1, timdid, salid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "temp1", nf_float, 1, timdid, temp1id)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "temp2", nf_float, 1, timdid, temp2id)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "temp3", nf_float, 1, timdid, temp3id)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "sal1", nf_float, 1, timdid, sal1id)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "sal2", nf_float, 1, timdid, sal2id)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "sal3", nf_float, 1, timdid, sal3id)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "dte", nf_float, 1, timdid, dteid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "dse", nf_float, 1, timdid, dseid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "terr", nf_float, 1, timdid, terrid)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_def_var(outid, "serr", nf_float, 1, timdid, serrid)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! Put attributes
  title = "Mk3L run " // run_name // ", years " // year1s // " to " // year2s
  write (history, "('Generated by ts_stats_128_112_21 on ', i4.4, '/', i2.2, &
                    '/', i2.2, ' at ', i2.2, ':', i2.2, ':', i2.2, ' ', a)") &
    date_time(1), date_time(2), date_time(3), date_time(5), date_time(6), &
    date_time(7), zone
  status = nf_put_att_text(outid, nf_global, "title", 36, title)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, nf_global, "history", 64, history)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, stfhtid, "long_name", 22, &
                           "Mean surface heat flux")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, stfsalid, "long_name", 30, &
                           "Mean surface salinity tendency")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, istfhtid, "long_name", 28, &
                           "Integrated surface heat flux")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, istfsalid, "long_name", 36, &
                           "Integrated surface salinity tendency")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, tempyzid, "long_name", 32, &
                           "Zonal-mean potential temperature")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, salyzid, "long_name", 19, &
                           "Zonal-mean salinity")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, tempzid, "long_name", 48, &
                           "Zonal- and meridional-mean potential temperature")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, salzid, "long_name", 48, &
                           "Zonal- and meridional-mean salinity")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, tempid, "long_name", 26, &
                           "Mean potential temperature")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, salid, "long_name", 13, &
                           "Mean salinity")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, temp1id, "long_name", 54, &
           "Mean potential temperature of the upper ocean (0-800m)")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, temp2id, "long_name", 55, &
           "Mean potential temperature of the mid-ocean (800-2350m)")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, temp3id, "long_name", 57, &
           "Mean potential temperature of the deep ocean (2350-4600m)")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, sal1id, "long_name", 41, &
           "Mean salinity of the upper ocean (0-800m)")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, sal2id, "long_name", 42, &
           "Mean salinity of the mid-ocean (800-2350m)")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, sal3id, "long_name", 44, &
           "Mean salinity of the deep ocean (2350-4600m)")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, dteid, "long_name", 47, &
           "Equivalent change in mean potential temperature")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, dseid, "long_name", 34, &
           "Equivalent change in mean salinity")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, terrid, "long_name", 35, &
           "Error in mean potential temperature")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, serrid, "long_name", 22, &
           "Error in mean salinity")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, latvid, "units", 13, "degrees_north")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, depvid, "units", 1, "m")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, timvid, "units", 5, "years")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, stfhtid, "units", 5, "W m-2")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, stfsalid, "units", 3, "s-1")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, istfhtid, "units", 1, "J")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, istfsalid, "units", 7, "kg kg-1")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, tempyzid, "units", 4, "degC")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, salyzid, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, tempzid, "units", 4, "degC")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, salzid, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, tempid, "units", 4, "degC")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, salid, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, temp1id, "units", 4, "degC")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, temp2id, "units", 4, "degC")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, temp3id, "units", 4, "degC")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, sal1id, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, sal2id, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, sal3id, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, dteid, "units", 1, "K")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, dseid, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, terrid, "units", 1, "K")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, serrid, "units", 3, "psu")
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_real(outid, tempyzid, "missing_value", nf_float, 1, &
                           missing_value4)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_real(outid, salyzid, "missing_value", nf_float, 1, &
                           missing_value4)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_att_text(outid, depvid, "positive", 4, "down")
  if (status /= nf_noerr) stop "***  netCDF error"

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

  ! Write axis data
  status = nf_put_var_real(outid, latvid, real(radian*phit(2:jmt-1), 4))
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_var_real(outid, depvid, zts)
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_var_int(outid, timvid, years)
  if (status /= nf_noerr) stop "***  netCDF error"

  ! -------------------------------
  ! Process the output of the model
  ! -------------------------------

  ! Loop over years, reading in one input file at a time
  do year = year1, year2, years_per_file
    iyear = year - year1 + 1

    ! Derive name of input file
    yr1 = year
    yr2 = year + years_per_file - 1
    write (yr1s, '(i5)') yr1
    write (yr2s, '(i5)') yr2
    if (yr1 < 10000) write (yr1s, '(a1, i4)'), "0", yr1
    if (yr1 < 1000) write (yr1s, '(a2, i3)'), "00", yr1
    if (yr1 < 100) write (yr1s, '(a3, i2)'), "000", yr1
    if (yr1 < 10) write (yr1s, '(a4, i1)'), "0000", yr1
    if (yr2 < 10000) write (yr2s, '(a1, i4)'), "0", yr2
    if (yr2 < 1000) write (yr2s, '(a2, i3)'), "00", yr2
    if (yr2 < 100) write (yr2s, '(a3, i2)'), "000", yr2
    if (yr2 < 10) write (yr2s, '(a4, i1)'), "0000", yr2
    infile = "com.ann." // run_name // "." // yr1s // "-" // yr2s // ".nc"
    write (*, *) "Reading data from file ", infile, " ..."

    ! Open input file
    status = nf_open(infile, nf_nowrite, inid)
    if (status /= nf_noerr) stop "***  netCDF error"

    ! Read data
    status = nf_inq_varid(inid, "stfht", datid)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_get_var_real(inid, datid, stfht_in)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_inq_varid(inid, "stfsal", datid)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_get_var_real(inid, datid, stfsal_in)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_inq_varid(inid, "temp", datid)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_get_var_real(inid, datid, temp_in)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_inq_varid(inid, "sal", datid)
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_get_var_real(inid, datid, sal_in)
    if (status /= nf_noerr) stop "***  netCDF error"

    ! If this is the first iteration of the loop, derive a three-dimensional
    ! land/sea mask and calculate the surface area and volume of the ocean
    if (year == year1) then

      ! (1) Derive a three-dimensional land/sea mask
      mask = .false.
      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            if (temp_in(i, j, k, 1) > -9000.0_4) mask(i, j, k) = .true.
          end do
        end do
      end do

      ! (2) Calculate the surface area of the ocean
      area = 0.0
      do j = 1, ny
        do i = 1, nx
          if (mask(i, j, 1)) area = area + dats(i, j)
        end do
      end do

      ! (3) Calculate the volume of the ocean
      volume = 0.0
      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            if (mask(i, j, k)) volume = volume + dvts(i, j, k)
          end do
        end do
      end do

      ! (4) Display the results
      write (*, *)
      write (*, *) "Surface area of ocean = ", area, " m^2"
      write (*, *) "Volume of ocean       = ", volume, " m^3"
      write (*, *)

    end if

    ! Calculate the mean surface fluxes
    stfht = 0.0
    stfsal = 0.0
    do j = 1, ny
      do i = 1, nx
        if (mask(i, j, 1)) then
          stfht = stfht +  dats(i, j) * stfht_in(i, j, :)
          stfsal = stfsal + dats(i, j) * stfsal_in(i, j, :)
        end if
      end do
    end do
    stfht = stfht / area
    stfsal = stfsal / area

    ! Integrate the surface fluxes with respect to time
    do l = 1, years_per_file
      if ((iyear == 1) .and. (l == 1)) then
        istfht(1) = secs_year * area * stfht(1)
        istfsal(1) = secs_year * stfsal(1)
      else
        istfht(iyear+l-1) = istfht(iyear+l-2) + secs_year * area * stfht(l)
        istfsal(iyear+l-1) = istfsal(iyear+l-2) + secs_year * stfsal(l)
     end if
    end do

    ! Calculate the zonal-mean temperature and salinity
    do k = 1, nz
      do j = 1, ny
        nvalues = 0
        tempyz(j, k, :) = 0.0
        salyz(j, k, :) = 0.0
        do i = 1, nx
          if (mask(i, j, k)) then
            nvalues = nvalues + 1
            tempyz(j, k, :) = tempyz(j, k, :) + temp_in(i, j, k, :)
            salyz(j, k, :) = salyz(j, k, :) + sal_in(i, j, k, :)
          end if
        end do
        if (nvalues == 0) then
          tempyz(j, k, :) = missing_value8
          salyz(j, k, :) = missing_value8
        else
          tempyz(j, k, :) = tempyz(j, k, :) / real(nvalues, 8)
          salyz(j, k, :) = salyz(j, k, :) / real(nvalues, 8)
        end if
      end do
    end do

    ! Calculate the zonal- and meridional-mean temperature and salinity
    do k = 1, nz
      vtotal = 0.0
      tempz(k, :) = 0.0
      salz(k, :) = 0.0
      do j = 1, ny
        do i = 1, nx
          if (mask(i, j, k)) then
            vtotal = vtotal + dvts(i, j, k)
            tempz(k, :) = tempz(k, :) + dvts(i, j, k) * temp_in(i, j, k, :)
            salz(k, :) = salz(k, :) + dvts(i, j, k) * sal_in(i, j, k, :)
          end if
        end do
      end do
      tempz(k, :) = tempz(k, :) / vtotal
      salz(k, :) = salz(k, :) / vtotal
    end do

    ! Calculate the mean temperature and salinity
    temp = 0.0
    sal = 0.0
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (mask(i, j, k)) then
            temp = temp + dvts(i, j, k) * temp_in(i, j, k, :)
            sal = sal + dvts(i, j, k) * sal_in(i, j, k, :)
          end if
        end do
      end do
    end do
    temp = temp / volume
    sal = sal / volume

    ! Calculate the mean temperature and salinity of the upper ocean
    temp1 = 0.0
    sal1 = 0.0
    vtotal = 0.0
    do k = 1, 11
      do j = 1, ny
        do i = 1, nx
          if (mask(i, j, k)) then
            vtotal = vtotal + dvts(i, j, k)
            temp1 = temp1 + dvts(i, j, k) * temp_in(i, j, k, :)
            sal1 = sal1 + dvts(i, j, k) * sal_in(i, j, k, :)
          end if
        end do
      end do
    end do
    temp1 = temp1 / vtotal
    sal1 = sal1 / vtotal

    ! Calculate the mean temperature and salinity of the mid-ocean
    temp2 = 0.0
    sal2 = 0.0
    vtotal = 0.0
    do k = 12, 16
      do j = 1, ny
        do i = 1, nx
          if (mask(i, j, k)) then
            vtotal = vtotal + dvts(i, j, k)
            temp2 = temp2 + dvts(i, j, k) * temp_in(i, j, k, :)
            sal2 = sal2 + dvts(i, j, k) * sal_in(i, j, k, :)
          end if
        end do
      end do
    end do
    temp2 = temp2 / vtotal
    sal2 = sal2 / vtotal

    ! Calculate the mean temperature and salinity of the deep ocean
    temp3 = 0.0
    sal3 = 0.0
    vtotal = 0.0
    do k = 17, 21
      do j = 1, ny
        do i = 1, nx
          if (mask(i, j, k)) then
            vtotal = vtotal + dvts(i, j, k)
            temp3 = temp3 + dvts(i, j, k) * temp_in(i, j, k, :)
            sal3 = sal3 + dvts(i, j, k) * sal_in(i, j, k, :)
          end if
        end do
      end do
    end do
    temp3 = temp3 / vtotal
    sal3 = sal3 / vtotal

    ! Calculate the errors in the mean temperature and salinity
    !
    ! (1) The equivalent change in mean temperature between one year and the
    !     next is given by [K]:
    !
    !     mean surface heat flux * secs_year * area / (volume * rhocw)
    !
    !     where rhocw = 4.1e6 J/m^3/K is the volumetric heat capacity of
    !     seawater.
    !
    ! (2) The equivalent change in mean salinity between one year and the next
    !     is given by [psu]:
    !
    !     1000.0 * mean surface salinity tendency * secs_year * 25.0 * area /
    !                                                                  volume
    !
    !     where 1000.0 converts from kg/kg to psu, and 25.0 m is the thickness
    !     of the upper layer of the ocean model
    do l = 1, years_per_file
      if ((iyear == 1) .and. (l == 1)) then
        dte(1) = 0.0
        dse(1) = 0.0
        terr(1) = 0.0
        serr(1) = 0.0
        dte_old = 0.0
        dse_old = 0.0
        terr_old = 0.0
        serr_old = 0.0
      else
        dte(l) = dte_old + 0.5 * secs_year * area * &
                           (stfht_old + stfht(l)) / (volume * rhocw)
        dse(l) = dse_old + 1.25e4 * secs_year * area * &
                           (stfsal_old + stfsal(l)) / volume
        terr(l) = terr_old + temp(l) - temp_old - dte(l) + dte_old
        serr(l) = serr_old + sal(l) - sal_old - dse(l) + dse_old
        dte_old = dte(l)
        dse_old = dse(l)
        terr_old = terr(l)
        serr_old = serr(l)
      end if
      stfht_old = stfht(l)
      stfsal_old = stfsal(l)
      temp_old = temp(l)
      sal_old = sal(l)
    end do

    ! Write data to output file
    status = nf_put_vara_real(outid, stfhtid, iyear, years_per_file, &
                              real(stfht, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, stfsalid, iyear, years_per_file, &
                              real(stfsal, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, tempyzid, (/ 1, 1, iyear /), &
                             (/ ny, nz, years_per_file /), real(tempyz, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, salyzid, (/ 1, 1, iyear /), &
                             (/ ny, nz, years_per_file /), real(salyz, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, tempzid, (/ 1, iyear /), &
                             (/ nz, years_per_file /), real(tempz, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, salzid, (/ 1, iyear /), &
                             (/ nz, years_per_file /), real(salz, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, tempid, iyear, years_per_file, &
                              real(temp, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, salid, iyear, years_per_file, &
                              real(sal, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, temp1id, iyear, years_per_file, &
                              real(temp1, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, temp2id, iyear, years_per_file, &
                              real(temp2, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, temp3id, iyear, years_per_file, &
                              real(temp3, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, sal1id, iyear, years_per_file, &
                              real(sal1, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, sal2id, iyear, years_per_file, &
                              real(sal2, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, sal3id, iyear, years_per_file, &
                              real(sal3, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, dteid, iyear, years_per_file, &
                              real(dte, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, dseid, iyear, years_per_file, &
                              real(dse, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, terrid, iyear, years_per_file, &
                              real(terr, 4))
    if (status /= nf_noerr) stop "***  netCDF error"
    status = nf_put_vara_real(outid, serrid, iyear, years_per_file, &
                              real(serr, 4))
    if (status /= nf_noerr) stop "***  netCDF error"

  end do

  ! Write integrated surface fluxes to output file
  status = nf_put_var_real(outid, istfhtid, real(istfht, 4))
  if (status /= nf_noerr) stop "***  netCDF error"
  status = nf_put_var_real(outid, istfsalid, real(istfsal, 4))
  if (status /= nf_noerr) stop "***  netCDF error"

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

  ! De-allocate output arrays
  deallocate (istfht, istfsal, years)

end program ts_stats_128_112_21