      subroutine globalsum (ind)

#if defined O_global_sums || defined O_co2emit_diag
!=======================================================================
!     calculate the global sums of heat and fresh water

!     input:
!       ind = index for calculation (1=initial, 2=current, 3=final)
!=======================================================================

      implicit none

      character(120) :: fname
      character(32) :: nstamp

      integer i, iou, id_time, ind, j, jrow, k, l, n, nc, ntrec
      integer nyear, nmonth, nday, nhour, nmin, nsec, it(10)
      integer layer

      logical defined

      real area, cstdyt, dtoih1, dtoih2, dtoic1, dtoic2
      real dtah1, dtah2, dtaf1, dtaf2, dtac1, dtac2
      real dtsh1, dtsh2, dtsf1, dtsf2, dtsc1, dtsc2
      real dtih1, dtih2, dtif1, dtif2, dtic1, dtic2
      real dtlh1, dtlh2, dtlf1, dtlf2, dtlc1, dtlc2
      real dtoh1, dtof1, dtoh2, dtof2, dtoc1, dtoc2
      real dth1,  dth2,  dtf1,  dtf2,  dtc1,  dtc2
      real tah(3), taf(3), tac(3), tsh(3), tsf(3), tsc(3), tih(3)
      real tif(3), tic(3), tlh(3), tlf(3), tlc(3), toh(3), tof(3)
      real toc(3), th, tf, tc, RFTIME, time, tmp, c1e20, vol, avgper

      save tah, taf, tac, tsh, tsf, tsc, tih, tif, tic, tlh, tlf, tlc
      save toh, tof, toc, dtoih1, dtoic1

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "cembm.h"
# if defined O_embm
      include "atm.h"
# endif
# if defined O_ice && defined O_embm
#  if defined O_ice_cpts
      include "cpts.h"
#  endif
      include "ice.h"
# endif
      include "calendar.h"
      include "coord.h"
      include "grdvar.h"
      include "levind.h"
      include "csbc.h"
      include "iounit.h"
      include "switch.h"
# if defined O_mom
      include "mw.h"
#  if defined O_npzd
      include "npzd.h"
#  endif
# endif
# if defined O_mtlm
      include "mtlm.h"
# endif
      include "tmngr.h"

!-----------------------------------------------------------------------
!     sum heat, water, carbon in atmosphere, snow, ice, land, ocean
!-----------------------------------------------------------------------

      c1e20 = 1.e20
      if (ind .eq. 1) then
        dtoih = 0.0
        dtoih1 = 0.0
        dtoic = 0.
        dtoic1 = 0.0
      endif

!     zero accumulators

      do n=min(3,ind*2-1),3
        tah(n) = 0.
        taf(n) = 0.
        tac(n) = 0.
        tsh(n) = 0.
        tsf(n) = 0.
        tsc(n) = 0.
        tih(n) = 0.
        tif(n) = 0.
        tic(n) = 0.
        tlh(n) = 0.
        tlf(n) = 0.
        tlc(n) = 0.
        toh(n) = 0.
        tof(n) = 0.
        toc(n) = 0.
      enddo

# if defined O_mtlm
#  if defined O_mtlm_segday
      RFTIME = DAY_TRIF*segtim/DAY_YEAR
#  else
      RFTIME = DAY_TRIF/DAY_YEAR
#  endif
# endif

!     calculate totals

      L = 0
      do jrow=2,jmtm1
# if defined O_mom
        if (wide_open_mw) then
          j = jrow
        else
          j = jmw
          call getrow (latdisk(taup1disk), nslab, jrow
     &,                u(1,1,j,1,taup1), t(1,1,j,1,taup1))
        endif
# endif
        cstdyt = cst(jrow)*dyt(jrow)
        do i=2,imtm1
          area = cstdyt*dxt(i)
          if (kmt(i,jrow) .eq. 0) L = L + 1

!         atmosphere
# if defined O_embm
          tah(3) = tah(3) + at(i,jrow,2,isat)*area
          taf(3) = taf(3) + at(i,jrow,2,ishum)*area
#  if defined O_carbon_co2_2d
          tac(3) = tac(3) + at(i,jrow,2,ico2)*area
#  else
          tac(3) = tac(3) + co2ccn*area
#  endif
# endif

!         snow and ice
# if defined O_ice && defined O_embm
          tsf(3) = tsf(3) + hsno(i,jrow,2)*area
          tif(3) = tif(3) + hice(i,jrow,2)*area
#  if defined O_ice_cpts
          do nc=1,ncat
            tsf(3) = tsf(3) + hseff(i,jrow,2,nc)*area
            tif(3) = tif(3) + heff(i,jrow,2,nc)*area
          enddo
          do layer=1,ntilay
            tih(3) = tih(3) + E(i,jrow,2,layer)*area
     &             + hice(i,jrow,2)*area
          enddo
#  endif
# endif
# if defined O_mtlm && defined O_embm
          if (kmt(i,jrow) .eq. 0)
     &      tsf(3) = tsf(3) + LYING_SNOW(L)*area*0.1/rhosno
# endif

# if defined O_embm
!         land
          tlf(3) = tlf(3) + soilm(i,jrow,2)*area*1.0e-3
#  if defined O_mtlm
          if (kmt(i,jrow) .eq. 0) then
            tlh(3) = tlh(3) + (TS1(L)*HCAP_SOIL*ROOTDEP
!           add in most recent fluxes not felt by atmosphere yet
     &             + (sbc(i,jrow,ievap)*vlocn + sbc(i,jrow,isens)
     &             + sbc(i,jrow,ilwr))*segtim*SEC_DAY*1.e-3)*area*1.e-4
            tlf(3) = tlf(3) + (M(L) + MNEG(L)
!           add back in most recent fluxes not felt by atmosphere yet
     &             + sbc(i,jrow,ievap)*segtim*SEC_DAY*10.)*area*1.e-4
            tlc(3) = tlc(3) + (CV(L) + CS(L)
#   if defined O_carbon
!           add in most recent fluxes not felt by atmosphere yet
#    if defined O_mtlm_segday
     &             + (sbc(i,jrow,isr) - sbc(i,jrow,inpp)
     &             + sbc(i,jrow,iburn)*(DAY_TRIF-ntlbc))*segtim*SEC_DAY
#    else
     &             + ((sbc(i,jrow,isr) - sbc(i,jrow,inpp))*segtim
     &             + sbc(i,jrow,iburn)*(DAY_TRIF-ntlbc*segtim))*SEC_DAY
#    endif
#   endif
!           add in driving fluxes not felt by triffid yet
     &             - RESP_S_DR(L)*RFTIME
     &               )*area*1.e-4
            do N=1,NPFT
              tlc(3) = tlc(3) + NPP_DR(L,N)*FRAC(L,N)*RFTIME*area*1.0e-4
            enddo
          endif
#  endif
# endif

!         ocean
# if defined O_mom
          do k=1,km
            vol = dzt(k)*area
            toh(3) = toh(3) + t(i,k,j,itemp,taup1)*vol
            tof(3) = tof(3) + t(i,k,j,isalt,taup1)*vol
#  if defined O_carbon
            toc(3) = toc(3) + t(i,k,j,idic,taup1)*vol
#   if defined O_npzd
            toc(3) = toc(3) + t(i,k,j,iphyt,taup1)*vol*redctn
            toc(3) = toc(3) + t(i,k,j,izoop,taup1)*vol*redctn
            toc(3) = toc(3) + t(i,k,j,idetr,taup1)*vol*redctn
#    if defined O_npzd_nitrogen
            toc(3) = toc(3) + t(i,k,j,idiaz,taup1)*vol*redctn
#    endif
#   endif
#  endif
          enddo
#  if defined O_sealev_data && defined O_sealev_salinity
          if (kmt(i,jrow) .ne. 0)
     &      tof(3) = tof(3) + sealev*area*socn
#  endif
# endif
        enddo
      enddo

!     convert units to Joules and kilograms

# if defined O_embm
!     atmosphere
      tah(3) = (taf(3)*rhoatm*shq*vlocn+tah(3)*cpatm*rhoatm*sht)*1.0e-7
      taf(3) = taf(3)*rhoatm*shq*1.0e-3
!     4.138e-7 => 12e-6 g/umol carbon / 29 g/mol air
      tac(3) = tac(3)*4.138e-7*rhoatm*shc*1.0e-3

# endif
# if defined O_ice && defined O_embm
!     snow and ice
      tsh(3) = -tsf(3)*rhosno*flice*1.0e-7
      tsf(3) = tsf(3)*rhosno*1.0e-3
#  if defined O_ice_cpts
      tih(3) = tih(3)*1.0e-7   ! cpts ice energy of melting is negative
#  else
      tih(3) = -tif(3)*rhoice*flice*1.0e-7
#  endif
      tif(3) = tif(3)*rhoice*1.0e-3

# endif
!     land
      tlh(3) = tlh(3)
      tlf(3) = tlf(3)
      tlc(3) = tlc(3)

# if defined O_mom
!     ocean
      toh(3) = toh(3)/0.2389
      tof(3) = -tof(3)/(socn*1.e3)
!     12e-9 kg/umol
      toc(3) = toc(3)*12.e-9

# endif
!     dtoic includes any sediment burial, weathering or emissions
!     12e-9 kg/umol
      dtoic = dtoic*12e-9

!     calculate differences from start of run

      dtah1 = tah(3) - tah(1)
      dtaf1 = taf(3) - taf(1)
      dtac1 = tac(3) - tac(1)
      dtsh1 = tsh(3) - tsh(1)
      dtsf1 = tsf(3) - tsf(1)
      dtsc1 = tsc(3) - tsc(1)
      dtih1 = tih(3) - tih(1)
      dtif1 = tif(3) - tif(1)
      dtic1 = tic(3) - tic(1)
      dtlh1 = tlh(3) - tlh(1)
      dtlf1 = tlf(3) - tlf(1)
      dtlc1 = tlc(3) - tlc(1)
      dtoh1 = toh(3) - toh(1)
      dtof1 = tof(3) - tof(1)
      dtoc1 = toc(3) - toc(1)
      dtoih1 = dtoih1 + dtoih
      dtoic1 = dtoic1 + dtoic

!     calculate differences from last time step

      dtah2 = tah(3) - tah(2)
      dtaf2 = taf(3) - taf(2)
      dtac2 = tac(3) - tac(2)
      dtsh2 = tsh(3) - tsh(2)
      dtsf2 = tsf(3) - tsf(2)
      dtsc2 = tsc(3) - tsc(2)
      dtih2 = tih(3) - tih(2)
      dtif2 = tif(3) - tif(2)
      dtic2 = tic(3) - tic(2)
      dtlh2 = tlh(3) - tlh(2)
      dtlf2 = tlf(3) - tlf(2)
      dtlc2 = tlc(3) - tlc(2)
      dtoh2 = toh(3) - toh(2)
      dtof2 = tof(3) - tof(2)
      dtoc2 = toc(3) - toc(2)
      dtoih2 = dtoih
      dtoic2 = dtoic

      dtoih = 0.0
      dtoic = 0.0
      do n=ind,2
        tah(n) = tah(3)
        taf(n) = taf(3)
        tac(n) = tac(3)
        tsh(n) = tsh(3)
        tsf(n) = tsf(3)
        tsc(n) = tsc(3)
        tih(n) = tih(3)
        tif(n) = tif(3)
        tic(n) = tic(3)
        tlh(n) = tlh(3)
        tlf(n) = tlf(3)
        tlc(n) = tlc(3)
        toh(n) = toh(3)
        tof(n) = tof(3)
        toc(n) = toc(3)
      enddo

!     write differences or totals

      if (ind .eq. 2) then

        dth1 = dtah1 + dtsh1 + dtih1 + dtlh1 + dtoh1 + dtoih1
        dtf1 = dtaf1 + dtsf1 + dtif1 + dtlf1 + dtof1
        dtc1 = dtac1 + dtsc1 + dtic1 + dtlc1 + dtoc1 + dtoic1
        dth2 = dtah2 + dtsh2 + dtih2 + dtlh2 + dtoh2 + dtoih2
        dtf2 = dtaf2 + dtsf2 + dtif2 + dtlf2 + dtof2
        dtc2 = dtac2 + dtsc2 + dtic2 + dtlc2 + dtoc2 + dtoic2

        if (tsits) then

          time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
          call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
          nyear = time
          call mkstmp (nstamp, nyear, nmonth, nday, nhour, nmin, nsec)

          call def_tsi
          call def_tsi_gsums (fname)

#  if !defined O_save_time_relyear0
!         make output time relative to year 1
          time = time - 1.
#  endif
          avgper = tsiper*accel
          if (avgper .le. 1e-6) avgper = 0.
#  if defined O_save_time_endper
          tmp = 0.
#  elif defined O_save_time_startper
          tmp = 1.
#  else
          tmp = 0.5
#  endif
#  if defined O_units_time_years
#   if defined O_calendar_360_day
          time = time - tmp*avgper/360.
#   elif defined O_calendar_gregorian
          time = time - tmp*avgper/365.25
#   else
          time = time - tmp*avgper/365.
#   endif
#  else
#   if defined O_calendar_360_day
          time = time*360. - tmp*avgper
#   elif defined O_calendar_gregorian
          time = time*365.25 - tmp*avgper
#   else
          time = time*365. - tmp*avgper
#   endif
#  endif

          call opennext (fname, time, ntrec, iou)
          if (ntrec .le. 0) ntrec = 1

          call putvars ('time', iou, ntrec, time, c1, c0)
          call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
          call putvars ('T_avgper', iou, ntrec, avgper, c1, c0)
# if defined O_co2emit_diag
          call putvars ('F_co2diag', iou, ntrec
#  if defined O_save_carbon_totals
     &,     dtc2, tsiint*1.e12/yrlen, c0)
#  else
     &,     dtc2, tsiint*daylen, c0)
#  endif
          carbemit = carbemit + dtc2*1e-12
          dtoic1 = dtoic1 - dtc1
          dtoic2 = dtoic2 - dtc2
          dtc1 = 0.
          dtc2 = 0.
# endif
# if defined O_global_sums
          call putvars ('G_dtheat', iou, ntrec, dth1, c1, c0)
          call putvars ('G_dtwater', iou, ntrec, dtf1, c1, c0)
          call putvars ('G_dtcarb', iou, ntrec, dtc1, c1, c0)
# endif

# if defined O_global_sums
          if (iotsi .eq. stdout .or. iotsi .lt. 0) then
            write (*,'(a,a)') 'Changes in heat, fresh water and carbon '
     &,       'from start of run and last calculation'
            write (*,'(a,e22.14,a,e22.14,a,e22.14,a
     &,                  e22.14,a,e22.14,a,e22.14,a)')
     &        ' d atm    ',dtah1, ' J ',dtaf1,' kg ',dtac1, ' kg'
     &,                    dtah2, ' J ',dtaf2,' kg ',dtac2, ' kg'
     &,       ' d snow   ',dtsh1, ' J ',dtsf1,' kg ',dtsc1, ' kg'
     &,                    dtsh2, ' J ',dtsf2,' kg ',dtsc2, ' kg'
     &,       ' d ice    ',dtih1, ' J ',dtif1,' kg ',dtic1, ' kg'
     &,                    dtih2, ' J ',dtif2,' kg ',dtic2, ' kg'
     &,       ' d lnd    ',dtlh1, ' J ',dtlf1,' kg ',dtlc1, ' kg'
     &,                    dtlh2, ' J ',dtlf2,' kg ',dtlc2, ' kg'
     &,       ' d ocn    ',dtoh1, ' J ',dtof1,' kg ',dtoc1, ' kg'
     &,                    dtoh2, ' J ',dtof2,' kg ',dtoc2, ' kg'
     &,       ' d out-in ',dtoih1,' J ',0.0,  ' kg ',dtoic1,' kg'
     &,                    dtoih2,' J ',0.0,  ' kg ',dtoic2,' kg'
     &,       ' d total  ',dth1,  ' J ',dtf1, ' kg ',dtc1,  ' kg'
     &,                    dth2,  ' J ',dtf2, ' kg ',dtc2,  ' kg'
          endif
        endif

      else

        th = tah(ind) + tsh(ind) + tih(ind) + tlh(ind) + toh(ind)
        tf = taf(ind) + tsf(ind) + tif(ind) + tlf(ind) + tof(ind)
        tc = tac(ind) + tsc(ind) + tic(ind) + tlc(ind) + toc(ind)
        write (*,'(/,a,a)') 'Total heat (in Joules referenced to 0 C'
     &,   ' and no ice or snow) and fresh water (in kg)'
#  if defined O_embm
        write (*,'(a,a)') 'Total ocean fresh water is the equivalent'
     &,   ' difference from the ocean volume referenced to socn'
#  else
        write (*,'(a,a)') 'Total ocean fresh water is the equivalent'
     &,   ' difference from the ocean volume referenced to 34.9 ppt'
#  endif
        write (*,'(a,e22.14,a,e22.14,a,e22.14,a)')
     &    ' t atm    ', tah(ind),' J ',taf(ind),' kg ',tac(ind),' kg'
     &,   ' t snow   ', tsh(ind),' J ',tsf(ind),' kg ',tsc(ind),' kg'
     &,   ' t ice    ', tih(ind),' J ',tif(ind),' kg ',tic(ind),' kg'
     &,   ' t lnd    ', tlh(ind),' J ',tlf(ind),' kg ',tlc(ind),' kg'
     &,   ' t ocn    ', toh(ind),' J ',tof(ind),' kg ',toc(ind),' kg'
     &,   ' t total  ', th,      ' J ',tf,      ' kg ',tc,      ' kg'

# else
        endif
# endif
      endif

      return
      end

      subroutine gsums_tsi_def (fname, calendar, expnam, runstamp)

!=======================================================================
!     output routine for atmospheric global sums step integrals

!   inputs:
!     fname      = file name
!     calendar   = calendar
!     expnam     = experiment name
!     runstamp   = run stamp
!=======================================================================

      implicit none

      character(*) :: fname, calendar, expnam, runstamp

      integer id(1), id_time, iou

      real c0, c1, c1e3, c1e20

      c0 = 0.
      c1 = 1.
      c1e3 = 1.e3
      c1e20 = 1.e20

!-----------------------------------------------------------------------
!     open file
!-----------------------------------------------------------------------
      call openfile (fname, iou)

!-----------------------------------------------------------------------
!     start definitions
!-----------------------------------------------------------------------
      call redef (iou)

!-----------------------------------------------------------------------
!     write global attributes
!-----------------------------------------------------------------------
      call putatttext (iou, 'global', 'Conventions', 'CF-1.0')
      call putatttext (iou, 'global', 'experiment_name', expnam)
      call putatttext (iou, 'global', 'run_stamp', runstamp)

!-----------------------------------------------------------------------
!     define dimensions
!-----------------------------------------------------------------------
      call defdim ('time', iou, 0, id_time)
      id(1) = id_time

!-----------------------------------------------------------------------
!     define data
!-----------------------------------------------------------------------
      call defvar ('time', iou, 1, id, c0, c0, 'T', 'D'
# if defined O_units_time_years
#  if !defined O_save_time_relyear0
     &, 'time', 'time', 'years since 1-1-1')
#  else
     &, 'time', 'time', 'years since 0-1-1')
#  endif
# else
#  if !defined O_save_time_relyear0
     &, 'time', 'time', 'days since 1-1-1')
#  else
     &, 'time', 'time', 'days since 0-1-1')
#  endif
# endif
      call putatttext (iou, 'time', 'calendar', calendar)
      call defvar ('T_avgper', iou, 1, id, c0, c0, ' ', 'F'
     &, 'averaging period', ' ','day')
# if defined O_co2emit_diag
      call defvar ('F_co2diag', iou, 1, id, -c1e20, c1e20, ' '
#  if defined O_save_carbon_totals
     &, 'F', 'diagnosed carbon emissions', ' ', 'Pg yr-1')
#  else
     &, 'F', 'diagnosed carbon emissions', ' ', 'kg s-1')
#  endif
# endif
# if defined O_global_sums
      call defvar ('G_dtheat', iou, 1, id, -c1e20
     &, c1e20, ' ', 'F', 'heat conservation error', ' ', 'J')
      call defvar ('G_dtwater', iou, 1, id, -c1e20
     &, c1e20, ' ', 'F', 'water conservation error', ' ', 'kg')
      call defvar ('G_dtcarb', iou, 1, id, -c1e20
     &, c1e20, ' ', 'F', 'carbon conservation error', ' ', 'kg')
# endif

!-----------------------------------------------------------------------
!     end definitions
!-----------------------------------------------------------------------
      call enddef (iou)

#endif
      return
      end
