! source file: /Users/oschlies/UVIC/master/testcase/updates/mom_rest.F
      subroutine mom_rest_in (fname, ids, ide, jds, jde)
!=======================================================================
!     input routine for atmospheric restarts

!     data may be sized differently in x and y from the global fields.
!     fields may be written with or without a time dimension. data
!     should be defined with the routine defvar and written with putvar.
!     if no time dimension, then data is only written once per file.
!     make sure the it, iu, ib, and ic arrays and are defining the
!     correct dimensions. ln may also need to be recalculated.

!   inputs:
!     fname              = file name
!     ids, ide ...       = start and end index for data domain

!     based on code by: M. Eby
!=======================================================================

!      implicit none

      include "param.h"
      include "coord.h"
      include "emode.h"
      include "grdvar.h"
      include "iounit.h"
      include "levind.h"
      include "mw.h"
      include "task_on.h"
      include "tmngr.h"
      include "switch.h"

      character(*) :: fname
      character(3) :: a3
      character(120) :: var1, var2

      integer i, iou, j, ln, n, ntrec, ids, ide, jds, jde, igs, ige, ig
      integer jgs, jge, jg, kgs, kge, kg, lgs, lge, lg, ils, ile, jls
      integer jle, kls, kle, lls, lle, ib(10), ic(10)
      integer nyear, nmonth, nday, nhour, nmin, nsec

      logical inqvardef, exists, warning

      real tmp
      real, allocatable :: tmpij(:,:), tmpik(:,:), tmpil(:,:)

!-----------------------------------------------------------------------
!     open file
!-----------------------------------------------------------------------
      call openfile (fname, iou)
      ntrec = 1

!-----------------------------------------------------------------------
!     local domain size (minimum of data domain and global read domain)
!-----------------------------------------------------------------------
      ils = max(ids,1)
      ile = min(ide,imt)
      jls = max(jds,1)
      jle = min(jde,jmt)
      kls = 1
      kle = km

      allocate ( tmpij(ils:ile,jls:jle) )
      allocate ( tmpik(ils:ile,kls:kle) )

!-----------------------------------------------------------------------
!     read 1d data (t)
!-----------------------------------------------------------------------
      tmp = itt
      call getvars ('itt', iou, 1, tmp, c1, c0)
      itt = tmp
      tmp = irstdy
      call getvars ('irstdy', iou, 1, tmp, c1, c0)
      irstdy = tmp
      tmp = msrsdy
      call getvars ('msrsdy', iou, 1, tmp, c1, c0)
      msrsdy = tmp
      tmp = year0
      call getvars ('year', iou, 1, tmp, c1, c0)
      nyear = tmp
      tmp = month0
      call getvars ('month', iou, 1, tmp, c1, c0)
      nmonth = tmp
      tmp = day0
      call getvars ('day', iou, 1, tmp, c1, c0)
      nday = tmp
      tmp = hour0
      call getvars ('hour', iou, 1, tmp, c1, c0)
      nhour = tmp
      tmp = min0
      call getvars ('minute', iou, 1, tmp, c1, c0)
      nmin = tmp
      tmp = sec0
      call getvars ('second', iou, 1, tmp, c1, c0)
      nsec = tmp
      call mkstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
      if (init_time_in) then
        itt = 0
        irstdy = 0
        msrsdy = 0
        relyr = 0.
        call mkstmp (stamp, year0, month0, day0, hour0, min0, sec0)
      endif

!-----------------------------------------------------------------------
!     read 2d data (x,y)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ile-ils+1
      ib(2) = 1
      ic(2) = jle-jls+1
      ln = ic(1)*ic(2)

!     check if kmt is consistent
      tmpij(:,:) = kmt(ils:ile,jls:jle)
      call getvara ('kmt', iou, ln, ib, ic, tmpij, c1, c0)
      warning = .false.
      do j=jls,jle
        do i=ils,ile
          if (kmt(i,j) .eq. 0 .and. int(tmpij(i,j)) .ne. 0) then
            print*, ' '
            print*, '==> Error: restart kmt is different. ocean point'
            print*, '           has been changed to land. this usually'
            print*, '           causes problems in the ocean solver.'
            print*, ' '
            stop '=>mom_rest_in'
          endif
          if (kmt(i,j) .ne. int(tmpij(i,j))) warning = .true.
        enddo
      enddo
      if (warning) then
        print*, ' '
        print*, '==> Warning: restart kmt is different. this often '
        print*, '             causes problems in the ocean solver.'
      endif

!-----------------------------------------------------------------------
!     read 3d data (x,y,t)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ile-ils+1
      ib(2) = 1
      ic(2) = jle-jls+1
      ib(3) = 1
      ic(3) = 1
      ln = ic(1)*ic(2)*ic(3)
      taum1disk = mod(itt+1,2) + 1
      taudisk   = mod(itt  ,2) + 1
      taup1disk = taum1disk
!     update pointers to tau-1, tau, & tau+1 data in the MW based on itt
      if (wide_open_mw) then
!       rotate time levels instead of moving data
        taum1 = mod(itt+0,3) - 1
        tau   = mod(itt+1,3) - 1
        taup1 = mod(itt+2,3) - 1
      endif

!     first do psi at "tau" then at "tau+1"
      tmpij(ils:ile,jls:jle) = tmpij(ils:ile,jls:jle)
      call getvara ('psi1', iou, ln, ib, ic, tmpij, c1, c0)
      psi(ils:ile,jls:jle,1) = tmpij(ils:ile,jls:jle)
      tmpij(ils:ile,jls:jle) = psi(ils:ile,jls:jle,2)
      call getvara ('psi2', iou, ln, ib, ic, tmpij, c1, c0)
      psi(ils:ile,jls:jle,2) = tmpij(ils:ile,jls:jle)
!     guess fields
      tmpij(ils:ile,jls:jle) = ptd(ils:ile,jls:jle)
      call getvara ('ptd1', iou, ln, ib, ic, tmpij, c1, c0)
      ptd(ils:ile,jls:jle) = tmpij(ils:ile,jls:jle)
      call oput (kflds, nwds, nkflds-1, ptd)
      tmpij(ils:ile,jls:jle) = ptd(ils:ile,jls:jle)
      call getvara ('ptd2', iou, ln, ib, ic, tmpij, c1, c0)
      ptd(ils:ile,jls:jle) = tmpij(ils:ile,jls:jle)
      call oput (kflds, nwds, nkflds, ptd)

!     construct depth arrays associated with "u" cells

      call depth_u (kmt, imt, jmt, zw, km, kmu, h, hr)

!-----------------------------------------------------------------------
!     read 4d data (x,y,z,t)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ile-ils+1
      ib(2) = 1
      ic(2) = 1
      ib(3) = 1
      ic(3) = kle-kls+1
      ib(4) = 1
      ic(4) = 1
      ln = ic(1)*ic(2)*ic(3)*ic(4)

!     read the "tau" latitude rows
      do jrow=1,jmt
        if (wide_open_mw) then
          j = jrow
        else
          j = jmw
          call getrow (latdisk(taudisk), nslab, jrow, u(1,1,j,1,tau)
     &,                                               t(1,1,j,1,tau))
        endif
        ib(2) = jrow
        do n=1,nt
          if (n .lt. 1000) write(a3, '(i3)') n
          if (n .lt. 100) write(a3, '(i2)') n
          if (n .lt. 10) write(a3, '(i1)') n
          var1 = 'tracer1_'//trim(a3)
          if (trim(mapt(n)) .eq. 'temp') then
            if (inqvardef('temp1', iou)) var1 = 'temp1'
          elseif (trim(mapt(n)) .eq. 'salt') then
            if (inqvardef('salt1', iou)) var1 = 'salt1'
          elseif (trim(mapt(n)) .eq. 'dic') then
            if (inqvardef('dic1', iou)) var1 = 'dic1'
          elseif (trim(mapt(n)) .eq. 'alk') then
            if (inqvardef('alk1', iou)) var1 = 'alk1'
          elseif (trim(mapt(n)) .eq. 'o2') then
            if (inqvardef('o21', iou)) var1 = 'o21'
          elseif (trim(mapt(n)) .eq. 'n2o') then
            if (inqvardef('n2o1', iou)) var1 = 'n2o1'
          elseif (trim(mapt(n)) .eq. 'abioto2') then
            if (inqvardef('abioto21', iou)) var1 = 'abioto21'
          elseif (trim(mapt(n)) .eq. 'prefo2') then
            if (inqvardef('prefo21', iou)) var1 = 'prefo21'
          elseif (trim(mapt(n)) .eq. 'prefpo4') then
            if (inqvardef('prefpo41', iou)) var1 = 'prefpo41'
          elseif (trim(mapt(n)) .eq. 'prefno3') then
            if (inqvardef('prefno31', iou)) var1 = 'prefno31'
          elseif (trim(mapt(n)) .eq. 'po4') then
            if (inqvardef('po41', iou)) var1 = 'po41'
          elseif (trim(mapt(n)) .eq. 'phyt') then
            if (inqvardef('phyt1', iou)) var1 = 'phyt1'
          elseif (trim(mapt(n)) .eq. 'zoop') then
            if (inqvardef('zoop1', iou)) var1 = 'zoop1'
          elseif (trim(mapt(n)) .eq. 'detr') then
            if (inqvardef('detr1', iou)) var1 = 'detr1'
          elseif (trim(mapt(n)) .eq. 'no3') then
            if (inqvardef('no31', iou)) var1 = 'no31'
          elseif (trim(mapt(n)) .eq. 'diaz') then
            if (inqvardef('diaz1', iou)) var1 = 'diaz1'
          elseif (trim(mapt(n)) .eq. 'c14') then
            if (inqvardef('c141', iou)) var1 = 'c141'
          elseif (trim(mapt(n)) .eq. 'cfc11') then
            if (inqvardef('cfc111', iou)) var1 = 'cfc111'
          elseif (trim(mapt(n)) .eq. 'cfc12') then
            if (inqvardef('cfc121', iou)) var1 = 'cfc121'
          endif
          exists = inqvardef(trim(var1), iou)
          if (exists) then
            tmpik(ils:ile,kls:kle) = t(ils:ile,kls:kle,j,n,tau)
            call getvara(trim(var1), iou, ln, ib, ic, tmpik, c1, c0)
            t(ils:ile,kls:kle,j,n,tau) = tmpik(ils:ile,kls:kle)
          else
            if (jrow .eq. 1) print*, '==> Warning: ', trim(var1)
     &,       ' not found'
          endif
        enddo
        exists = inqvardef('u1', iou)
        if (exists) then
          tmpik(ils:ile,kls:kle) = u(ils:ile,kls:kle,j,1,tau)
          call getvara('u1', iou, ln, ib, ic, tmpik, c1, c0)
          u(ils:ile,kls:kle,j,1,tau) = tmpik(ils:ile,kls:kle)
        else
          if (jrow .eq. 1) print*, '==> Warning: u1 not found'
        endif
        exists = inqvardef('v1', iou)
        if (exists) then
          tmpik(ils:ile,kls:kle) = u(ils:ile,kls:kle,j,2,tau)
          call getvara('v1', iou, ln, ib, ic, tmpik, c1, c0)
          u(ils:ile,kls:kle,j,2,tau) = tmpik(ils:ile,kls:kle)
        else
          if (jrow .eq. 1) print*, '==> Warning: v1 not found'
        endif
!       initialize every latitude
        if (wide_open_mw) then
!         do nothing since "tau" data is in place in the MW
        else
          call putrow (latdisk(taudisk), nslab, jrow, u(1,1,j,1,tau)
     &,                                               t(1,1,j,1,tau))
        endif
      enddo
      if (wide_open_mw) then
!       Initialize 1st and last latitude row for tau-1 to prevent
!       use of uninitialized values on boundary row.
        do j=1,jmt,jmt-1
          do k=1,km
            do i=1,imt
              u(i,k,j,1,taum1) = u(i,k,j,1,tau)
              u(i,k,j,2,taum1) = u(i,k,j,2,tau)
              do n=1,nvar-2
                t(i,k,j,n,taum1) = t(i,k,j,n,tau)
              enddo
            enddo
          enddo
        enddo
      endif

!     read the "tau+1" latitude rows
      do jrow=1,jmt
        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
        ib(2) = jrow
        do n=1,nt
          if (n .lt. 1000) write(a3, '(i3)') n
          if (n .lt. 100) write(a3, '(i2)') n
          if (n .lt. 10) write(a3, '(i1)') n
          var2 = 'tracer2_'//trim(a3)
          if (trim(mapt(n)) .eq. 'temp') then
            if (inqvardef('temp2', iou)) var2 = 'temp2'
          elseif (trim(mapt(n)) .eq. 'salt') then
            if (inqvardef('salt2', iou)) var2 = 'salt2'
          elseif (trim(mapt(n)) .eq. 'dic') then
            if (inqvardef('dic2', iou)) var2 = 'dic2'
          elseif (trim(mapt(n)) .eq. 'alk') then
            if (inqvardef('alk2', iou)) var2 = 'alk2'
          elseif (trim(mapt(n)) .eq. 'o2') then
            if (inqvardef('o22', iou)) var2 = 'o22'
          elseif (trim(mapt(n)) .eq. 'n2o') then
            if (inqvardef('n2o2', iou)) var2 = 'n2o2'
          elseif (trim(mapt(n)) .eq. 'abioto2') then
            if (inqvardef('abioto22', iou)) var2 = 'abioto22'
          elseif (trim(mapt(n)) .eq. 'prefo2') then
            if (inqvardef('prefo22', iou)) var2 = 'prefo22'
          elseif (trim(mapt(n)) .eq. 'prefpo4') then
            if (inqvardef('prefpo42', iou)) var2 = 'prefpo42'
          elseif (trim(mapt(n)) .eq. 'prefno3') then
            if (inqvardef('prefno32', iou)) var2 = 'prefno32'
          elseif (trim(mapt(n)) .eq. 'po4') then
            if (inqvardef('po42', iou)) var2 = 'po42'
          elseif (trim(mapt(n)) .eq. 'phyt') then
            if (inqvardef('phyt2', iou)) var2 = 'phyt2'
          elseif (trim(mapt(n)) .eq. 'zoop') then
            if (inqvardef('zoop2', iou)) var2 = 'zoop2'
          elseif (trim(mapt(n)) .eq. 'detr') then
            if (inqvardef('detr2', iou)) var2 = 'detr2'
          elseif (trim(mapt(n)) .eq. 'no3') then
            if (inqvardef('no32', iou)) var2 = 'no32'
          elseif (trim(mapt(n)) .eq. 'diaz') then
            if (inqvardef('diaz2', iou)) var2 = 'diaz2'
          elseif (trim(mapt(n)) .eq. 'c14') then
            if (inqvardef('c142', iou)) var2 = 'c142'
          elseif (trim(mapt(n)) .eq. 'cfc11') then
            if (inqvardef('cfc112', iou)) var2 = 'cfc112'
          elseif (trim(mapt(n)) .eq. 'cfc12') then
            if (inqvardef('cfc122', iou)) var2 = 'cfc122'
          endif
          exists = inqvardef(trim(var2), iou)
          if (exists) then
            tmpik(ils:ile,kls:kle) = t(ils:ile,kls:kle,j,n,taup1)
            call getvara(trim(var2), iou, ln, ib, ic, tmpik, c1, c0)
            t(ils:ile,kls:kle,j,n,taup1) = tmpik(ils:ile,kls:kle)
          else
            if (jrow .eq. 1) print*, '==> Warning: ',trim(var2)
     &,       ' not found'
          endif
        enddo
        exists = inqvardef('u2', iou)
        if (exists) then
          tmpik(ils:ile,kls:kle) = u(ils:ile,kls:kle,j,1,taup1)
          call getvara('u2', iou, ln, ib, ic, tmpik, c1, c0)
          u(ils:ile,kls:kle,j,1,taup1) = tmpik(ils:ile,kls:kle)
        else
          if (jrow .eq. 1) print*, '==> Warning: u2 not found'
        endif
        exists = inqvardef('v2', iou)
        if (exists) then
          tmpik(ils:ile,kls:kle) = u(ils:ile,kls:kle,j,2,taup1)
          call getvara('v2', iou, ln, ib, ic, tmpik, c1, c0)
          u(ils:ile,kls:kle,j,2,taup1) = tmpik(ils:ile,kls:kle)
        else
          if (jrow .eq. 1) print*, '==> Warning: v2 not found'
        endif
!       initialize every latitude
        if (wide_open_mw) then
!         do nothing since "tau+1" data is in place in the MW
        else
          call putrow (latdisk(taup1disk), nslab, jrow
     &,                u(1,1,j,1,taup1), t(1,1,j,1,taup1))
        endif
      enddo

!-----------------------------------------------------------------------
!     close the file
!-----------------------------------------------------------------------
      print*, '=> Ocn restart read from ',trim(fname),' on ', stamp

      deallocate (tmpij)
      deallocate (tmpik)

      call closefile (iou)

      return
      end

      subroutine mom_rest_def (fname)
!=======================================================================
!     definition routine for atmospheric restarts

!   inputs:
!     fname = file name

!     based on code by: M. Eby
!=======================================================================

!      implicit none

      include "param.h"
      include "iounit.h"
      include "mw.h"
      include "tmngr.h"

      character(*) :: fname
      character(3) :: a3
      character(32) :: lstamp

      integer i, iou, j, ln, n, ntrec, igs, ige, ig, jgs, jge, jg, kgs
      integer kge, kg, lgs, lge, lg, it(10), iu(10), id_time, id_xt
      integer id_xu, id_yt, id_yu, id_zt, id_zw, id_zl, id_xt_e, id_xu_e
      integer id_yt_e, id_yu_e, id_zt_e, id_zw_e, id_zl_e, id_part

      real c100, c1e3, c1e20

      c100 = 100.
      c1e3 = 1.e3
      c1e20 = 1.e20
      lstamp = stamp

!-----------------------------------------------------------------------
!     open file
!-----------------------------------------------------------------------
      call openfile (fname, iou)
      ntrec = 1

!-----------------------------------------------------------------------
!     set global write domain size
!-----------------------------------------------------------------------
      igs = 1
      ige = imt
      ig  = ige-igs+1
      jgs = 1
      jge = jmt
      jg  = jge-jgs+1
      kgs = 1
      kge = km
      kg  = kge-kgs+1

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

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

!-----------------------------------------------------------------------
!     define dimensions
!-----------------------------------------------------------------------
      call defdim ('time', iou, 0, id_time)
      call defdim ('xt', iou, ig, id_xt)
      call defdim ('yt', iou, jg, id_yt)
      call defdim ('zt', iou, kg, id_zt)
      call defdim ('xu', iou, ig, id_xu)
      call defdim ('yu', iou, jg, id_yu)
      call defdim ('zw', iou, kg, id_zw)
      call defdim ('xt_edges', iou, ig+1, id_xt_e)
      call defdim ('yt_edges', iou, jg+1, id_yt_e)
      call defdim ('zt_edges', iou, kg+1, id_zt_e)
      call defdim ('xu_edges', iou, ig+1, id_xu_e)
      call defdim ('yu_edges', iou, jg+1, id_yu_e)
      call defdim ('zw_edges', iou, kg+1, id_zw_e)

!-----------------------------------------------------------------------
!     define 1d data (t)
!-----------------------------------------------------------------------
      it(1) = id_time
      call defvar ('time', iou, 1, it, c0, c0, 'T', 'D'
     &, 'time since initial condition', 'time', trim(timunit))
      call defvar ('itt', iou, 1, it, c0, c0, ' ', 'D'
     &, 'itt', ' ',' ')
      call defvar ('irstdy', iou, 1, it, c0, c0, ' ', 'D'
     &, 'irstdy', ' ',' ')
      call defvar ('msrsdy', iou, 1, it, c0, c0, ' ', 'D'
     &, 'msrsdy', ' ',' ')
      call defvar ('year', iou, 1, it, c0, c0, ' ', 'D'
     &, 'year', ' ',' ')
      call defvar ('month', iou, 1, it, c0, c0, ' ', 'D'
     &, 'month', ' ',' ')
      call defvar ('day', iou, 1, it, c0, c0, ' ', 'D'
     &, 'day', ' ',' ')
      call defvar ('hour', iou, 1, it, c0, c0, ' ', 'D'
     &, 'hour', ' ',' ')
      call defvar ('minute', iou, 1, it, c0, c0, ' ', 'D'
     &, 'minute', ' ',' ')
      call defvar ('second', iou, 1, it, c0, c0, ' ', 'D'
     &, 'second', ' ',' ')

!-----------------------------------------------------------------------
!     define 1d data (x, y or z)
!-----------------------------------------------------------------------
      it(1) = id_xt
      call defvar ('xt', iou, 1, it, c0, c0, 'X', 'D'
     &, 'longitude of the t grid', 'grid_longitude', 'degrees_east')
      it(1) = id_yt
      call defvar ('yt', iou, 1, it, c0, c0, 'Y', 'D'
     &, 'latitude of the t grid', 'grid_latitude', 'degrees_north')
      it(1) = id_zt
      call defvar ('zt', iou, 1, it, c0, c0, 'Z', 'D'
     &, 'depth of the t grid', 'depth', 'm')
      it(1) = id_xu
      call defvar ('xu', iou, 1, it, c0, c0, 'X', 'D'
     &, 'longitude of the u grid', 'grid_longitude', 'degrees_east')
      it(1) = id_yu
      call defvar ('yu', iou, 1, it, c0, c0, 'Y', 'D'
     &, 'latitude of the u grid', 'grid_latitude', 'degrees_north')
      it(1) = id_zw
      call defvar ('zw', iou, 1, it, c0, c0, 'Z', 'D'
     &, 'depth of the w grid', 'depth', 'm')
      it(1) = id_xt_e
      call defvar ('xt_edges', iou, 1, it, c0, c0, ' ', 'D'
     &, 'longitude of t grid edges', ' ', 'degrees')
      it(1) = id_yt_e
      call defvar ('yt_edges', iou, 1, it, c0, c0, ' ', 'D'
     &, 'latitude of t grid edges', ' ', 'degrees')
      it(1) = id_zt_e
      call defvar ('zt_edges', iou, 1, it, c0, c0, ' ', 'D'
     &, 'depth of t grid edges', ' ', 'm')
      it(1) = id_xu_e
      call defvar ('xu_edges', iou, 1, it, c0, c0, ' ', 'D'
     &, 'longitude of u grid edges', ' ', 'degrees')
      it(1) = id_yu_e
      call defvar ('yu_edges', iou, 1, it, c0, c0, ' ', 'D'
     &, 'latitude of u grid edges', ' ', 'degrees')
      it(1) = id_zw_e
      call defvar ('zw_edges', iou, 1, it, c0, c0, ' ', 'D'
     &, 'depth of w grid edges', ' ', 'm')

!-----------------------------------------------------------------------
!     define 2d data (x,y)
!-----------------------------------------------------------------------
      it(1) = id_xt
      iu(1) = id_xu
      it(2) = id_yt
      iu(2) = id_yu
      call defvar ('kmt', iou, 2, it, c0, c1e3, ' ', 'I'
     &,  'kmt', ' ' ,' ')

!-----------------------------------------------------------------------
!     define 3d data (x,y,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      iu(1) = id_xu
      it(2) = id_yt
      iu(2) = id_yu
      it(3) = id_time
      iu(3) = id_time
      call defvar ('psi1', iou, 3, it,  -c1e20, c1e20, ' ', 'D'
     &, 'psi1', ' ', ' ')
      call defvar ('psi2', iou, 3, it,  -c1e20, c1e20, ' ', 'D'
     &, 'psi2', ' ', ' ')
      call defvar ('ptd1', iou, 3, it,  -c1e20, c1e20, ' ', 'D'
     &, 'ptd1', ' ', ' ')
      call defvar ('ptd2', iou, 3, it,  -c1e20, c1e20, ' ', 'D'
     &, 'ptd2', ' ', ' ')

!-----------------------------------------------------------------------
!     define 4d data (x,y,z,t)
!-----------------------------------------------------------------------
      it(1) = id_xt
      iu(1) = id_xu
      it(2) = id_yt
      iu(2) = id_yu
      it(3) = id_zt
      iu(3) = id_zt
      it(4) = id_time
      iu(4) = id_time
      do n=1,nt
        if (trim(mapt(n)) .eq. 'temp') then
          call defvar ('temp1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'ocean potential temperature at tau', ' ', ' ')
          call defvar ('temp2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'ocean potential temperature at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'salt') then
          call defvar ('salt1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'ocean salinity at tau', '  ', ' ')
          call defvar ('salt2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'ocean salinity at tau+1', '  ', ' ')
        elseif (trim(mapt(n)) .eq. 'dic') then
          call defvar ('dic1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'ocean carbon at tau', ' ', ' ')
          call defvar ('dic2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'ocean carbon at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'alk') then
          call defvar ('alk1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'alkalinity at tau',' ', ' ')
          call defvar ('alk2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'alkalinity at tau+1',' ', ' ')
        elseif (trim(mapt(n)) .eq. 'o2') then
          call defvar ('o21', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'oxygen at tau',' ', ' ')
          call defvar ('o22', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'oxygen at tau+1',' ', ' ')
       elseif (trim(mapt(n)) .eq. 'n2o') then
          call defvar ('n2o1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'N2O at tau',' ', ' ')
          call defvar ('n2o2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'N2O at tau+1',' ', ' ')
       elseif (trim(mapt(n)) .eq. 'abioto2') then
          call defvar ('abioto21', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'abiotO2 at tau',' ', ' ')
          call defvar ('abioto22', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'abiotO2 at tau+1',' ', ' ')
       elseif (trim(mapt(n)) .eq. 'prefo2') then
          call defvar ('prefo21', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'prefO2 at tau',' ', ' ')
          call defvar ('prefo22', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'prefO2 at tau+1',' ', ' ')
       elseif (trim(mapt(n)) .eq. 'prefpo4') then
          call defvar ('prefpo41', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'prefPO4 at tau',' ', ' ')
          call defvar ('prefpo42', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'prefPO4 at tau+1',' ', ' ')
       elseif (trim(mapt(n)) .eq. 'prefno3') then
          call defvar ('prefno31', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'prefNO3 at tau',' ', ' ')
          call defvar ('prefno32', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'prefNO3 at tau+1',' ', ' ')
        elseif (trim(mapt(n)) .eq. 'po4') then
          call defvar ('po41', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'phosphate at tau', ' ', ' ')
          call defvar ('po42', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'phosphate at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'phyt') then
          call defvar ('phyt1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'phytoplankton at tau', ' ', ' ')
          call defvar ('phyt2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'phytoplankton at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'zoop') then
          call defvar ('zoop1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'zooplankton at tau', ' ', ' ')
          call defvar ('zoop2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'zooplankton at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'detr') then
          call defvar ('detr1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'detritus at tau', ' ', ' ')
          call defvar ('detr2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'detritus at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'no3') then
          call defvar ('no31', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'nitrate at tau', ' ', ' ')
          call defvar ('no32', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'nitrate at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'diaz') then
          call defvar ('diaz1', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'diazotrophs at tau', ' ', ' ')
          call defvar ('diaz2', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'diazotrophs at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'c14') then
          call defvar ('c141', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'carbon 14 at tau', ' ', ' ')
          call defvar ('c142', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'carbon 14 at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'cfc11') then
          call defvar ('cfc111', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'CFC11 at tau', ' ', ' ')
          call defvar ('cfc112', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'CFC11 at tau+1', ' ', ' ')
        elseif (trim(mapt(n)) .eq. 'cfc12') then
          call defvar ('cfc121', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'CFC12 at tau', ' ', ' ')
          call defvar ('cfc122', iou, 4, it, -c1e20, c1e20, ' ', 'D'
     &,     'CFC12 at tau+1', ' ', ' ')
        else
          if (n .lt. 1000) write(a3,'(i3)') n
          if (n .lt. 100) write(a3,'(i2)') n
          if (n .lt. 10) write(a3,'(i1)') n
          call defvar ('tracer1_'//trim(a3), iou ,4, it, -c1e6, c1e6
     &,     ' ', 'D', 'tracer '//trim(a3)//' at tau', ' ', ' ')
          call defvar ('tracer2_'//trim(a3), iou ,4, it, -c1e6, c1e6
     &,     ' ', 'D', 'tracer '//trim(a3)//' at tau+1', ' '
     &,     ' ')
        endif
      enddo
      call defvar ('u1', iou , 4, iu, -c1e20, c1e20, ' '
     &, 'D', 'u1', ' ', ' ')
      call defvar ('u2', iou , 4, iu, -c1e20, c1e20, ' '
     &, 'D', 'u2', ' ', ' ')
      call defvar ('v1', iou , 4, iu, -c1e20, c1e20, ' '
     &, 'D', 'v1', ' ', ' ')
      call defvar ('v2', iou , 4, iu, -c1e20, c1e20, ' '
     &, 'D', 'v2', ' ', ' ')

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

!-----------------------------------------------------------------------
!     close the file
!-----------------------------------------------------------------------
      call closefile (iou)

      return
      end

      subroutine mom_rest_out (fname, ids, ide, jds, jde)
!=======================================================================
!     output routine for atmospheric restarts

!     data may be sized differently in x and y from the global fields.
!     fields may be written with or without a time dimension. data
!     should be defined with the routine defvar and written with putvar.
!     if no time dimension, then data is only written once per file.
!     make sure the it, iu, ib, and ic arrays and are defining the
!     correct dimensions. ln may also need to be recalculated.

!   inputs:
!     fname              = file name
!     ids, ide ...       = start and end index for data domain

!     based on code by: M. Eby
!=======================================================================

!      implicit none

      include "param.h"
      include "coord.h"
      include "emode.h"
      include "grdvar.h"
      include "iounit.h"
      include "levind.h"
      include "mw.h"
      include "switch.h"
      include "task_on.h"
      include "tmngr.h"

      character(*) :: fname
      character(3) :: a3
      character(32) :: lstamp
      character(120) :: var1, var2

      integer i, iou, j, ln, n, ntrec, ids, ide, jds, jde, igs, ige, ig
      integer jgs, jge, jg, kgs, kge, kg, lgs, lge, lg, ils, ile, jls
      integer jle, kls, kle, lls, lle, ib(10), ic(10), it(10), iu(10)
      integer id_time, id_xt, id_xu, id_yt, id_yu, id_zt, id_zw, id_zl
      integer id_xt_e, id_xu_e, id_yt_e, id_yu_e, id_zt_e, id_zw_e
      integer id_zl_e, id_part, nyear, nmonth, nday, nhour, nmin, nsec

      real xt_e(imt+1), xu_e(imt+1), yt_e(jmt+1), yu_e(jmt+1)
      real zt_e(km+1), zw_e(km+1), time, tmp, bufsl(imt,km,2)
      real ext(imt,2), c100, c1e3, c1e20
      real, allocatable :: tmpij(:,:), tmpik(:,:), tmpil(:,:)
      real, allocatable :: tmpi(:), tmpj(:), tmpk(:), tmpl(:)
      real, allocatable :: tmpie(:), tmpje(:), tmpke(:), tmple(:)

      c100 = 100.
      c1e3 = 1.e3
      c1e20 = 1.e20
      lstamp = stamp

!-----------------------------------------------------------------------
!     open file
!-----------------------------------------------------------------------
      call openfile (fname, iou)
      ntrec = 1

!-----------------------------------------------------------------------
!     set global write domain size
!-----------------------------------------------------------------------
      igs = 1
      ige = imt
      ig  = ige-igs+1
      jgs = 1
      jge = jmt
      jg  = jge-jgs+1
      kgs = 1
      kge = km
      kg  = kge-kgs+1

!-----------------------------------------------------------------------
!     local domain size (minimum of data domain and global write domain)
!-----------------------------------------------------------------------
      ils = max(ids,igs)
      ile = min(ide,ige)
      jls = max(jds,jgs)
      jle = min(jde,jge)
      kls = max(1,kgs)
      kle = min(km,kge)

      allocate ( tmpij(ils:ile,jls:jle) )
      allocate ( tmpik(ils:ile,kls:kle) )
      allocate ( tmpi(igs:ige) )
      allocate ( tmpj(jgs:jge) )
      allocate ( tmpk(kgs:kge) )
      allocate ( tmpie(igs:ige+1) )
      allocate ( tmpje(jgs:jge+1) )
      allocate ( tmpke(kgs:kge+1) )

!-----------------------------------------------------------------------
!     write 1d data (t)
!-----------------------------------------------------------------------
      if (init_time_out) then
        tmp = 0.
        call putvars ('time', iou, ntrec, 0., c1, c0)
        tmp = 0.
        call putvars ('itt', iou, ntrec, tmp, c1, c0)
        tmp = 0.
        call putvars ('irstdy', iou, ntrec, tmp, c1, c0)
        tmp = 0.
        call putvars ('msrsdy', iou, ntrec, tmp, c1, c0)
        call mkstmp (lstamp, year0, month0, day0, hour0, min0, sec0)
      else
        tmp = relyr
        call putvars ('time', iou, ntrec, 0., c1, c0)
        tmp = itt
        call putvars ('itt', iou, ntrec, tmp, c1, c0)
        tmp = iday(imodeltime)
        call putvars ('irstdy', iou, ntrec, tmp, c1, c0)
        tmp = msday(imodeltime)
        call putvars ('msrsdy', iou, ntrec, tmp, c1, c0)
      endif
      call rdstmp (lstamp, nyear, nmonth, nday, nhour, nmin, nsec)
      tmp = nyear
      call putvars ('year', iou, 1, tmp, c1, c0)
      tmp = nmonth
      call putvars ('month', iou, 1, tmp, c1, c0)
      tmp = nday
      call putvars ('day', iou, 1, tmp, c1, c0)
      tmp = nhour
      call putvars ('hour', iou, 1, tmp, c1, c0)
      tmp = nmin
      call putvars ('minute', iou, 1, tmp, c1, c0)
      tmp = nsec
      call putvars ('second', iou, 1, tmp, c1, c0)

!-----------------------------------------------------------------------
!     write 1d data (x, y or z)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ig
      tmpi(igs:ige) = xt(igs:ige)
      call putvara ('xt', iou, ig, ib, ic, tmpi, c1, c0)
      tmpi(igs:ige) = xu(igs:ige)
      call putvara ('xu', iou, ig, ib, ic, tmpi, c1, c0)

      ic(1) = jg
      tmpj(jgs:jge) = yt(jgs:jge)
      call putvara ('yt', iou, jg, ib, ic, tmpj, c1, c0)
      tmpj(jgs:jge) = yu(jgs:jge)
      call putvara ('yu', iou, jg, ib, ic, tmpj, c1, c0)

      ic(1) = kg
      tmpk(kgs:kge) = zt(kgs:kge)
      call putvara ('zt', iou, kg, ib, ic, tmpk, c100, c0)
      tmpk(kgs:kge) = zw(kgs:kge)
      call putvara ('zw', iou, kg, ib, ic, tmpk, c100, c0)

      ic(1) = ig + 1
      call edge_maker (1, xt_e, xt, dxt, xu, dxu, imt)
      tmpie(igs:ige+1) = xt_e(igs:ige+1)
      call putvara ('xt_edges', iou, ig+1, ib, ic, tmpie, c1, c0)
      call edge_maker (2, xu_e, xt, dxt, xu, dxu, imt)
      tmpie(igs:ige+1) = xu_e(igs:ige+1)
      call putvara ('xu_edges', iou, ig+1, ib, ic, tmpie, c1, c0)

      ic(1) = jg + 1
      call edge_maker (1, yt_e, yt, dyt, yu, dyu, jmt)
      tmpje(jgs:jge+1) = yt_e(jgs:jge+1)
      call putvara ('yt_edges', iou, jg+1, ib, ic, tmpje, c1, c0)
      call edge_maker (2, yu_e, yt, dyt, yu, dyu, jmt)
      tmpje(jgs:jge+1) = yu_e(jgs:jge+1)
      call putvara ('yu_edges', iou, jg+1, ib, ic, tmpje, c1, c0)

      ic(1) = kg + 1
      call edge_maker (1, zt_e, zt, dzt, zw, dzw, km)
      tmpke(kgs:kge+1) = zt_e(kgs:kge+1)
      call putvara ('zt_edges', iou, kg+1, ib, ic, tmpke, c100, c0)
      call edge_maker (2, zw_e, zt, dzt, zw, dzw, km)
      tmpke(kgs:kge+1) = zw_e(kgs:kge+1)
      call putvara ('zw_edges', iou, kg+1, ib, ic, tmpke, c100, c0)

!-----------------------------------------------------------------------
!     write 2d data (x,y)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ile-ils+1
      ib(2) = 1
      ic(2) = jle-jls+1
      ln = ic(1)*ic(2)
      tmpij(ils:ile,jls:jle) = kmt(ils:ile,jls:jle)
      call putvara ('kmt', iou, ln, ib, ic, tmpij, c1, c0)

!-----------------------------------------------------------------------
!     write 3d data (x,y,t)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ile-ils+1
      ib(2) = 1
      ic(2) = jle-jls+1
      ib(3) = ntrec
      ic(3) = 1
      ln = ic(1)*ic(2)*ic(3)
!     first do psi at "tau" then at "tau+1"
      tmpij(ils:ile,jls:jle) = psi(ils:ile,jls:jle,1)
      call putvara ('psi1', iou, ln, ib, ic, tmpij, c1, c0)
      tmpij(ils:ile,jls:jle) = psi(ils:ile,jls:jle,2)
      call putvara ('psi2', iou, ln, ib, ic, tmpij, c1, c0)
!     guess fields
      call oget (kflds, nwds, nkflds-1, ptd)
      tmpij(ils:ile,jls:jle) = ptd(ils:ile,jls:jle)
      call putvara ('ptd1', iou, ln, ib, ic, tmpij, c1, c0)
      call oget (kflds, nwds, nkflds, ptd)
      tmpij(ils:ile,jls:jle) =  ptd(ils:ile,jls:jle)
      call putvara ('ptd2', iou, ln, ib, ic, tmpij, c1, c0)

!-----------------------------------------------------------------------
!     write 4d data (x,y,z,t)
!-----------------------------------------------------------------------
      ib(1) = 1
      ic(1) = ile-ils+1
      ib(2) = 1
      ic(2) = 1
      ib(3) = 1
      ic(3) = kle-kls+1
      ib(4) = 1
      ic(4) = 1
      ln = ic(1)*ic(2)*ic(3)*ic(4)

!     save the "tau" latitude rows
      do jrow=1,jmt
        if (wide_open_mw) then
          j = jrow
!         remove external mode from "tau". since psi has been updated
!         psi(,,2) is at "tau"
          if (jrow .lt. jmt) then
            do i=2,imt-1
              diag1 = psi(i+1,jrow+1,2) - psi(i  ,jrow,2)
              diag0 = psi(i  ,jrow+1,2) - psi(i+1,jrow,2)
              ext(i,1) = -(diag1+diag0)*dyu2r(jrow)*hr(i,jrow)
              ext(i,2) =  (diag1-diag0)*dxu2r(i)*hr(i,jrow)*csur(jrow)
             enddo
            do k=1,km
              do i=2,imt-1
                if (k .le. kmu(i,jrow)) then
                  bufsl(i,k,1) = (u(i,k,j,1,tau) - ext(i,1))
                  bufsl(i,k,2) = (u(i,k,j,2,tau) - ext(i,2))
                else
                  bufsl(i,k,1) = c0
                  bufsl(i,k,2) = c0
                endif
              enddo
            enddo
            call setbcx (bufsl(1,1,1), imt, km)
            call setbcx (bufsl(1,1,2), imt, km)
          else
            do k=1,km
              do i=1,imt
                bufsl(i,k,1) = c0
                bufsl(i,k,2) = c0
              enddo
            enddo
          endif
        else
          j = jmw
          call getrow (latdisk(taudisk), nslab, jrow, u(1,1,j,1,tau)
     &,                                               t(1,1,j,1,tau))
          do k=1,km
            do i=1,imt
              bufsl(i,k,1) = u(i,k,j,1,tau)
              bufsl(i,k,2) = u(i,k,j,2,tau)
            enddo
          enddo
        endif

        ib(2) = jrow
        do n=1,nt
          if (trim(mapt(n)) .eq. 'temp') then
            var1 = 'temp1'
          elseif (trim(mapt(n)) .eq. 'salt') then
            var1 = 'salt1'
          elseif (trim(mapt(n)) .eq. 'dic') then
            var1 = 'dic1'
          elseif (trim(mapt(n)) .eq. 'alk') then
            var1 = 'alk1'
          elseif (trim(mapt(n)) .eq. 'o2') then
            var1 = 'o21'
          elseif (trim(mapt(n)) .eq. 'n2o') then
            var1 = 'n2o1'
          elseif (trim(mapt(n)) .eq. 'abioto2') then
            var1 = 'abioto21'
          elseif (trim(mapt(n)) .eq. 'prefo2') then
            var1 = 'prefo21'
          elseif (trim(mapt(n)) .eq. 'prefpo4') then
            var1 = 'prefpo41'
          elseif (trim(mapt(n)) .eq. 'prefno3') then
            var1 = 'prefno31'
          elseif (trim(mapt(n)) .eq. 'po4') then
            var1 = 'po41'
          elseif (trim(mapt(n)) .eq. 'phyt') then
            var1 = 'phyt1'
          elseif (trim(mapt(n)) .eq. 'zoop') then
            var1 = 'zoop1'
          elseif (trim(mapt(n)) .eq. 'detr') then
            var1 = 'detr1'
          elseif (trim(mapt(n)) .eq. 'no3') then
            var1 = 'no31'
          elseif (trim(mapt(n)) .eq. 'diaz') then
            var1 = 'diaz1'
          elseif (trim(mapt(n)) .eq. 'c14') then
            var1 = 'c141'
          elseif (trim(mapt(n)) .eq. 'cfc11') then
            var1 = 'cfc111'
          elseif (trim(mapt(n)) .eq. 'cfc12') then
            var1 = 'cfc121'
          else
            if (n .lt. 1000) write(a3, '(i3)') n
            if (n .lt. 100) write(a3, '(i2)') n
            if (n .lt. 10) write(a3, '(i1)') n
            var1 = 'tracer1_'//trim(a3)
          endif
          tmpik(ils:ile,kls:kle) = t(ils:ile,kls:kle,j,n,tau)
          call putvara(trim(var1), iou, ln, ib, ic, tmpik, c1, c0)
        enddo
        tmpik(ils:ile,kls:kle) = bufsl(ils:ile,kls:kle,1)
        call putvara('u1', iou, ln, ib, ic, tmpik, c1, c0)
        tmpik(ils:ile,kls:kle) = bufsl(ils:ile,kls:kle,2)
        call putvara('v1', iou, ln, ib, ic, tmpik, c1, c0)
      enddo

!     save the "tau+1" latitude rows
      do jrow=1,jmt
        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
        ib(2) = jrow
        do n=1,nt
          if (trim(mapt(n)) .eq. 'temp') then
            var2 = 'temp2'
          elseif (trim(mapt(n)) .eq. 'salt') then
            var2 = 'salt2'
          elseif (trim(mapt(n)) .eq. 'dic') then
            var2 = 'dic2'
          elseif (trim(mapt(n)) .eq. 'alk') then
            var2 = 'alk2'
          elseif (trim(mapt(n)) .eq. 'o2') then
            var2 = 'o22'
          elseif (trim(mapt(n)) .eq. 'n2o') then
            var2 = 'n2o2'
          elseif (trim(mapt(n)) .eq. 'abioto2') then
            var2 = 'abioto22'
          elseif (trim(mapt(n)) .eq. 'prefo2') then
            var2 = 'prefo22'
          elseif (trim(mapt(n)) .eq. 'prefpo4') then
            var2 = 'prefpo42'
          elseif (trim(mapt(n)) .eq. 'prefno3') then
            var2 = 'prefno32'
          elseif (trim(mapt(n)) .eq. 'po4') then
            var2 = 'po42'
          elseif (trim(mapt(n)) .eq. 'phyt') then
            var2 = 'phyt2'
          elseif (trim(mapt(n)) .eq. 'zoop') then
            var2 = 'zoop2'
          elseif (trim(mapt(n)) .eq. 'detr') then
            var2 = 'detr2'
          elseif (trim(mapt(n)) .eq. 'no3') then
            var2 = 'no32'
          elseif (trim(mapt(n)) .eq. 'diaz') then
            var2 = 'diaz2'
          elseif (trim(mapt(n)) .eq. 'c14') then
            var2 = 'c142'
          elseif (trim(mapt(n)) .eq. 'cfc11') then
            var2 = 'cfc112'
          elseif (trim(mapt(n)) .eq. 'cfc12') then
            var2 = 'cfc122'
          else
            if (n .lt. 1000) write(a3, '(i3)') n
            if (n .lt. 100) write(a3, '(i2)') n
            if (n .lt. 10) write(a3, '(i1)') n
            var2 = 'tracer2_'//trim(a3)
          endif
          tmpik(ils:ile,kls:kle) = t(ils:ile,kls:kle,j,n,taup1)
          call putvara(trim(var2), iou, ln, ib, ic, tmpik, c1, c0)
        enddo
        tmpik(ils:ile,kls:kle) = u(ils:ile,kls:kle,j,1,taup1)
        call putvara('u2', iou, ln, ib, ic, tmpik, c1, c0)
        tmpik(ils:ile,kls:kle) = u(ils:ile,kls:kle,j,2,taup1)
        call putvara('v2', iou, ln, ib, ic, tmpik, c1, c0)
      enddo

!-----------------------------------------------------------------------
!     close the file
!-----------------------------------------------------------------------
      print*, '=> Ocn restart written to ',trim(fname),' on ', lstamp

      deallocate ( tmpij )
      deallocate ( tmpik )
      deallocate ( tmpi )
      deallocate ( tmpj )
      deallocate ( tmpk )
      deallocate ( tmpie )
      deallocate ( tmpje )
      deallocate ( tmpke )

      call closefile (iou)

      return
      end