! source file: /nfs/tape_cache/smomw258/latin_hypercube_nobio/template/updates/UVic_ESCM.F

      program UVic_ESCM

!=======================================================================

!     UNIVERSITY OF VICTORIA EARTH SYSTEM CLIMATE MODEL (UVic ESCM)

!     The UVic ESCM is a climate model developed by researchers in the
!     Climate Modelling Group, within the School of Earth and Ocean
!     Sciences, located at the University of Victoria, Victoria,
!     British Columbia, Canada.

!     Many people have contributed to the development of this code.
!     This is a collective effort and although individual contributions
!     are appreciated, individual authorship is not indicated.

!     Some of the people who have contributed to the development of the
!     model are:
!       D. Archer, C. Avis, A. Berger, R. Betts, A. Biastoch, C. Bitz,
!       J. Brauch, C. Brennan, B. Bryan, J. Burton, S. Carto,
!       M. Cottet-Puinel, M. Cox, P. Cox, G. Danabasoglu, K. Dixon,
!       J. Dukowicz, M. Eby, R. Essery, T. Ewen, A. Fanning, J. Fyke,
!       R. Gerdes, A. Gnanadesikan, C. Goldberg, D. Goldberg,
!       D. Gregory, J. Gregory, S. Griffies, R. Hanson, R. Hetherington,
!       H. Hickey, M. Holland, G. Holloway, T. Huck, T. Hughes,
!       E. Hunke, W. Hurlin, W. Ingram, R. Key, E. Kluzek, C. Koeberle,
!       K. Kvale, C. Lawson, J. Lewis, M. Loutre, R. Malone,
!       D. Matthews, K. Meissner, A. Montenegro, A. Mouchet, T. Murdock,
!       P. Myers, A. Oschlies, R. Pacanowski, M. Pahlow, P. Poussart
!       S. Rahmstorf, R. Redler,, D. Robitaille, A. Rosati, M. Roth
!       C. Sabine, O. Saenko, A. Schmittner, B. Semtner, W. Sijp
!       T. Silva, H. Simmons, A. Skvortsov, R. Smith, P. Spence
!       D. Stone, R. Tonkonojenkov, C. Tricot, S. Valcke, A. Weaver
!       E. Wiebe, J. Willebrandt, M. Yoshimori, K. Zickfeld

!     Please direct problems or questions to the code contact
!     person at: http://climate.uvic.ca/model

!     Requirements:

!     Standard fortran 90 is used

!     Disclaimer:

!     The UVic Earth System Climate Model is a climate modelling
!     research tool developed at the University of Victoria. Others may
!     use it freely but we assume no responsibility for problems or
!     incorrect use. It is left to the user to ensure that a particular
!     configuration is working correctly.
!=======================================================================

      implicit none

      integer i, j, numots, numats, numseg, n, loop

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "csbc.h"
      include "iounit.h"
      include "emode.h"
      include "levind.h"
      include "scalar.h"
      include "switch.h"
      include "tmngr.h"
      include "cembm.h"

      include "atm.h"

      include "mw.h"

      print*, '== UNIVERSITY OF VICTORIA EARTH SYSTEM CLIMATE MODEL =='
      print*, '                                            '

!-----------------------------------------------------------------------
!     initialize i/o units
!-----------------------------------------------------------------------

      call ioinit

!-----------------------------------------------------------------------
!     setup file renaming
!-----------------------------------------------------------------------

      call file_names

!-----------------------------------------------------------------------
!     Initialize S.B.C. indices
!-----------------------------------------------------------------------

      call sbc_init

!-----------------------------------------------------------------------
!     Initialize tracers
!-----------------------------------------------------------------------

      call tracer_init

!-----------------------------------------------------------------------
!     read namelist variables
!-----------------------------------------------------------------------

      call read_namelist

!-----------------------------------------------------------------------
!     read grid
!-----------------------------------------------------------------------

      call grids

!-----------------------------------------------------------------------
!     read topography
!-----------------------------------------------------------------------

      call topog (kmt, kmu, map, xt, yt, zt, xu, yu, zw, imt, jmt, km

     &,           sg_bathy

     &            )

      call isleperim (kmt, map, iperm, jperm, iofs, nippts, nisle, imt
     &,               jmt, km, mnisle, maxipp, xu, yu, zw)

!-----------------------------------------------------------------------
!     common setup
!-----------------------------------------------------------------------

      call setcom (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!     ocean setup
!-----------------------------------------------------------------------

      call setmom (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!     atmosphere setup
!-----------------------------------------------------------------------
      call setembm (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!     land setup (requires embm)
!-----------------------------------------------------------------------
      call setmtlm (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!     initialize data reading routine
!-----------------------------------------------------------------------

      call setdata (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!     compute the number of ocean time steps "numots" for this run and
!     the number of ocean time steps per ocean segment "ntspos".
!     compute the number of atmos time steps "numats" for this run and
!     the number of atmos time steps per atmos segment "ntspas".
!     divide the integration time "days" into "numseg" segments.
!     each will be length "segtim" days. Surface boundary conditions
!     are supplied every "segtim" days.
!-----------------------------------------------------------------------

      numots = nint(rundays/(dtocn*secday))
      ntspos = nint(segtim/(dtocn*secday))
      numats = nint(rundays/(dtatm*secday))
      ntspas = nint(segtim/(dtatm*secday))
      numseg = numots/ntspos

      if (segtim .gt. 1.) then
        ntspls = nint(c1/(dtlnd*secday))
      else
        ntspls = nint(segtim/(dtlnd*secday))
      endif

!-----------------------------------------------------------------------
!     check for consistency in the S.B.C. setup
!-----------------------------------------------------------------------

      call chkcpl

!-----------------------------------------------------------------------
!     get global sums at the start of the run
!-----------------------------------------------------------------------

      call globalsum (1)

!-----------------------------------------------------------------------
!     S T A R T    S E G M E N T    L O O P
!-----------------------------------------------------------------------

      tsicpl = .false.
      tavgcpl = .false.
      do n=1,numseg
        cplts = .false.

!-----------------------------------------------------------------------
!       get the atmospheric S.B.C.
!-----------------------------------------------------------------------

        call gasbc (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!       call the atmospheric model once for each time step until one
!       segment of "segtim" days is complete. hold atmos S.B.C. fixed
!       during each segment and predict average S.B.C. for ocean
!-----------------------------------------------------------------------

        do loop=1,ntspas

          call embm (1, imt, 1, jmt)

        enddo

!-----------------------------------------------------------------------
!       get land S.B.C.s
!-----------------------------------------------------------------------

        call glsbc (1, imt, 1, jmt)

!----------------------------------------------------------------------
!       call the land-surface and vegetation  model once for each time
!       step until one segment of "segtim" days is complete.
!-----------------------------------------------------------------------

        do loop=1,ntspls
          call mtlm (1, imt, 1, jmt)
        enddo

!-----------------------------------------------------------------------
!       get ocean S.B.C.s
!-----------------------------------------------------------------------

        call gosbc (1, imt, 1, jmt)

!-----------------------------------------------------------------------
!       call the ocean model once for each time step until one
!       segment of "segtim" days is complete. hold ocean S.B.C. fixed
!       during each segment and predict average S.B.C. for atmos
!-----------------------------------------------------------------------

        do loop=1,ntspos
          if (loop .eq. ntspos) cplts = .true.
          call mom

          call embmout (1, imt, 1, jmt)

          call mtlmout (1, imt, 1, jmt)

          if (tsits .and. iotsi .ne. stdout .and. iotsi .gt. 0) then
            write (*,'(1x, a3, i7, 1x, a32)') 'ts=',itt, stamp
          endif
        enddo

!-----------------------------------------------------------------------
!       write change in global sums for heat and fresh water
!-----------------------------------------------------------------------

        if (tsits) call globalsum (2)

!-----------------------------------------------------------------------
!       close any open netcdf files to flush buffers
!-----------------------------------------------------------------------

        call closeall

      enddo

!-----------------------------------------------------------------------
!     E N D    S E G M E N T    L O O P
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!     get global sums at the end of the run
!-----------------------------------------------------------------------

      call globalsum (3)

      print*, ' ==>  UVIC_ESCM integration is complete.'

      call closeall

      call release_all

      stop
      end

      subroutine chkcpl

      implicit none

      integer jrow, j, i, k

      logical errorc

      real critv, tmp, t1, r1, r2, r3, r4, r5, r6

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"
      include "cembm.h"
      include "calendar.h"
      include "switch.h"
      include "tmngr.h"

!-----------------------------------------------------------------------
!     adjust output intervals for forcing acceleration
!-----------------------------------------------------------------------

      if (accel .le. 1) then
        accel = 1.
      else
        print*, ' '
        print*, '==> Warning: any forcing that changes on time scales'
     &        ,             ' longer than 1 year has been accelerated'
        print*, '             acceleration factor: ', accel
        print*, '             starts at relyr: ', accel_yr0
        print*, '             current year: ', year0 + relyr
        tmp = year0 + accel_yr0 + (relyr - accel_yr0)*accel
        print*, '             current accelerated year: ', tmp
        print*, ' '
        print*, '==> Warning: adjusting output intervals for'
     &,  ' long term forcing acceleration'
        print*, ' '
        if (tsiint .ge. yrlen*accel) tsiint = tsiint/accel
        if (tavgint .ge. yrlen*accel) tavgint = tavgint/accel
        if (tmbint .ge. yrlen*accel) tmbint = tmbint/accel
        if (stabint .ge. yrlen*accel) stabint = stabint/accel
        if (zmbcint .ge. yrlen*accel) zmbcint = zmbcint/accel
        if (glenint .ge. yrlen*accel) glenint = glenint/accel
        if (trmbint .ge. yrlen*accel) trmbint = trmbint/accel
        if (vmsfint .ge. yrlen*accel) vmsfint = vmsfint/accel
        if (gyreint .ge. yrlen*accel) gyreint = gyreint/accel
        if (extint .ge. yrlen*accel) extint = extint/accel
        if (prxzint .ge. yrlen*accel) prxzint = prxzint/accel
        if (exconvint .ge. yrlen*accel) exconvint = exconvint/accel
        if (dspint .ge. yrlen*accel) dspint = dspint/accel
        if (timavgint .ge. yrlen*accel) timavgint = timavgint/accel
        if (cmixint .ge. yrlen*accel) cmixint = cmixint/accel
        if (xbtint .ge. yrlen*accel) xbtint = xbtint/accel
        if (crossint .ge. yrlen*accel) crossint = crossint/accel
        if (densityint .ge. yrlen*accel) densityint = densityint/accel
        if (fctint .ge. yrlen*accel) fctint = fctint/accel
        if (tyzint .ge. yrlen*accel) tyzint = tyzint/accel
        if (restint .ge. yrlen*accel) restint = restint/accel
        if (tbtint .ge. yrlen*accel) tbtint = tbtint/accel
      endif

      if (timavgint .lt. timavgper) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: averaging period "timavgper" exceeds interval  '
     &,   '             "timavgint". Setting timavgper = timavgint     '
        timavgper = timavgint
      endif
      write (stdout,'(/,1x,a,f10.2,a,/,1x,a,f10.2,a)')
     &  '==> Time averages written every ', timavgint, ' days, '
     &, '    with an averaging period of ', timavgper, ' days. '

      if (tsiint .lt. tsiper) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: averaging period "tsiper" exceeds interval  '
     &,   '             "tsiint". Setting tsiper = tsiint     '
        tsiper = tsiint
      endif
      write (stdout,'(/,1x,a,f10.2,a,/,1x,a,f10.2,a)')
     &  '==> Time step integrals written every ', tsiint, ' days, '
     &, '    with an averaging period of       ', tsiper, ' days. '

      if (mod(namix,2) .ne. mod(nint(segtim*daylen/dtatm),2)) then
        write(*,*) '==> Error: time steps between mixing and coupling'
        write(*,*) '           in the atmosphere must both be even to'
        write(*,*) '           use even_fluxes.'
        stop '=>UVic_ESCM'
      endif

      if (mod(nmix,2) .ne. mod(nint(segtim*daylen/dtocn),2)) then
        write(*,*) '==> Error: time steps between mixing and coupling'
        write(*,*) '           in the ocean must both be even to use'
        write(*,*) '           even_fluxes.'
        stop '=>UVic_ESCM'
      endif

      if (mod(namix,2) .ne. mod(nats,2)) then
        write(*,*) '==> Warning: restart was not saved with even flux'
        write(*,*) '             averaging. Starting with a mixing time'
        write(*,*) '             step in atmosphere.'
        nats = namix
      endif

      if (mod(nmix,2) .ne. mod(nots,2)) then
        write(*,*) '==> Warning: restart was not saved with even flux'
        write(*,*) '             averaging. Starting with a mixing time'
        write(*,*) '             step in ocean.'
        nots = nmix
      endif

!-----------------------------------------------------------------------
!     do consistency checks before allowing model to continue
!-----------------------------------------------------------------------

      errorc = .false.
      write (stdout,*) ' '
      write (stdout,*) '    (checking setup)'

      if (dtatm .eq. c0) then
          write (stdout,9000)
     & '==> Error: the atmospheric time step must be set in "setatm"'
          errorc = .true.
          dtatm = 1.e-6
      endif
!      critv = 1.e-6
      critv = 1.e-4
      if (segtim .ne. c0) then
        r1 = rundays/segtim
      else
        r1 = 0.5
      endif
      if (segtim .eq. c0) then
          write (stdout,9000)
     & '==> Error: coupling period "segtim" must be specified'
          errorc = .true.
      endif
      if (abs(r1-nint(r1)) .gt. critv) then
          write (stdout,9000)
     & '==> Error: there must be an integral number of segments '
     &,'"segtim"  within "rundays" (the length of the run)'
          errorc = .true.
      endif

      r2 = segtim/(dtocn*secday)
      if (abs(r2-nint(r2)) .gt. critv) then
        write (stdout,9000)
     & '==> Error: there must be an integral number of density time '
     &,'steps "dtocn"  within "segtim" (the segment time)'
        errorc = .true.
      endif

      r3 = segtim/(dtatm*secday)
      if (abs(r3-nint(r3)) .gt. critv) then
        write (stdout,9000)
     & '==> Error: there must be an integral number of atmos time '
     &,'steps "dtatm"  within "segtim" (the segment time)'
        errorc = .true.
      endif

      r5 = segtim/(dtlnd*secday)
      if (abs(r5-nint(r5)) .gt. critv) then
        write (stdout,9000)
     & '==> Error: there must be an integral number of land time '
     &,'steps "dtlnd"  within "segtim" (the segment time)'
        errorc = .true.
      endif

          tmp = 0

     &        + 1

      if (tmp .gt. 1) then
        write (stdout,9000)
     & '==> Error: use only one of O_co2ccn_data, O_co2emit_data, '
     &,'or O_co2ccn_user'
        errorc = .true.
      endif

      write (stdout,*) ' '
      write (stdout,*) '    (End of checks) '
      write (stdout,*) ' '
      if (errorc) stop '=> ERRORS found in chkcpl'

9000  format (/,(1x,a))
      return
      end

      subroutine set (index, num, name, text, inc)

!-----------------------------------------------------------------------
!     increment counter, set index and text
!-----------------------------------------------------------------------

      character(*) :: name, text

      name = text
      index = num
      inc = index + 1

      return
      end

      subroutine getst (jrow, ocnout, ntabc)

!-----------------------------------------------------------------------
!     read surface tracers from disk row "jrow"
!-----------------------------------------------------------------------

      implicit none

      integer i, jrow, ntabc

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "iounit.h"
      include "mw.h"
      include "tmngr.h"

      real ocnout(imt,jmt)

      call getrow (latdisk(taup1disk), nslab, jrow
     &,          u(1,1,jmw,1,taup1), t(1,1,jmw,1,taup1))
      do i=1,imt
        ocnout(i,jrow) = t(i,1,jmw,ntabc,taup1)
      enddo

      return
      end

      subroutine sbc_init

!-----------------------------------------------------------------------
!     Initialize S.B.C. indices
!-----------------------------------------------------------------------

      implicit none

      integer m

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"

      sbc(:,:,:) = 0.0
      mapsbc(:) = " "
      itaux = 0
      itauy = 0
      iws = 0
      iaca = 0
      isca = 0
      ihflx = 0
      isflx = 0
      isst = 0
      isss = 0
      iro = 0
      iwa = 0
      iwxq = 0
      iwyq = 0
      iwxt = 0
      iwyt = 0
      iwxc = 0
      iwyc = 0
      ipsw = 0
      isu = 0
      isv = 0
      igu = 0
      igv = 0
      issdic = 0
      idicflx = 0
      issalk = 0
      ialkflx = 0
      isso2 = 0
      io2flx = 0
      isspo4 = 0
      ipo4flx = 0
      issphyt = 0
      iphytflx = 0
      isszoop = 0
      izoopflx = 0
      issdetr = 0
      idetrflx = 0

      issmp = 0
      impflx = 0
      issmpa = 0
      impaflx = 0
      issmpp = 0
      imppflx = 0

      issdetrz = 0
      idetrzflx = 0

      issno3 = 0
      ino3flx = 0
      issdiaz = 0
      idiazflx = 0

      isscocc = 0
      icoccflx = 0

      icaco3flx = 0
      isscaco3 = 0

      issdetr_B = 0
      idetrflx_B = 0

      issc14 = 0
      ic14flx = 0
      isscfc11 = 0
      icfc11flx = 0
      isscfc12 = 0
      icfc12flx = 0
      iat = 0
      irh = 0
      ipr = 0
      ips = 0
      iaws = 0
      iswr = 0
      ilwr = 0
      isens = 0
      ievap = 0
      idtr = 0
      isr = 0
      inpp = 0
      iburn = 0
      ibtemp = 0
      ibsalt = 0
      ibdic = 0
      ibdicfx = 0
      ibalk = 0
      ibalkfx = 0
      ibo2 = 0
      ircal = 0
      irorg = 0

      m = 1
      call set (itaux, m, mapsbc(m), 'taux', m)
      call set (itauy, m, mapsbc(m), 'tauy', m)
      call set (iws, m, mapsbc(m), 'ws', m)
      call set (iaca, m, mapsbc(m), 'a_calb', m)
      call set (isca, m, mapsbc(m), 's_calb', m)
      call set (ihflx, m, mapsbc(m), 'hflx', m)
      call set (isflx, m, mapsbc(m), 'sflx', m)
      call set (isst, m, mapsbc(m), 'sst', m)
      call set (isss, m, mapsbc(m), 'sss', m)
      call set (iro, m, mapsbc(m), 'ro', m)

      call set (iwxq, m, mapsbc(m), 'wx_q', m)
      call set (iwyq, m, mapsbc(m), 'wy_q', m)
      call set (iwxt, m, mapsbc(m), 'wx_t', m)
      call set (iwyt, m, mapsbc(m), 'wy_t', m)

      call set (isu, m, mapsbc(m), 'su', m)
      call set (isv, m, mapsbc(m), 'sv', m)
      call set (igu, m, mapsbc(m), 'gu', m)
      call set (igv, m, mapsbc(m), 'gv', m)

      call set (issdic, m, mapsbc(m), 'ssdic', m)
      call set (idicflx, m, mapsbc(m), 'dicflx', m)

      call set (issc14, m, mapsbc(m), 'ssc14', m)
      call set (ic14flx, m, mapsbc(m), 'c14flx', m)

      call set (issmp, m, mapsbc(m), 'ssmp', m)
      call set (impflx, m, mapsbc(m), 'mpflx', m)
      call set (issmpa, m, mapsbc(m), 'ssmpa', m)
      call set (impaflx, m, mapsbc(m), 'mpaflx', m)
      call set (issmpp, m, mapsbc(m), 'ssmpp', m)
      call set (imppflx, m, mapsbc(m), 'mppflx', m)

      call set (issalk, m, mapsbc(m), 'ssalk', m)
      call set (ialkflx, m, mapsbc(m), 'alkflx', m)

      call set (isso2, m, mapsbc(m), 'sso2', m)
      call set (io2flx, m, mapsbc(m), 'o2flx', m)

      call set (isspo4, m, mapsbc(m), 'sspo4', m)
      call set (ipo4flx, m, mapsbc(m), 'po4flx', m)

      call set (issno3, m, mapsbc(m), 'ssno3', m)
      call set (ino3flx, m, mapsbc(m), 'no3flx', m)

      call set (iat, m, mapsbc(m), 'at', m)
      call set (irh, m, mapsbc(m), 'rh', m)
      call set (ipr, m, mapsbc(m), 'pr', m)
      call set (ips, m, mapsbc(m), 'ps', m)
      call set (iaws, m, mapsbc(m), 'aws', m)
      call set (iswr, m, mapsbc(m), 'swr', m)
      call set (ilwr, m, mapsbc(m), 'lwr', m)
      call set (isens, m, mapsbc(m), 'sens', m)
      call set (ievap, m, mapsbc(m), 'evap', m)
      call set (idtr, m, mapsbc(m), 'dtr', m)

      call set (isr, m, mapsbc(m), 'sr', m)
      call set (inpp, m, mapsbc(m), 'npp', m)
      call set (iburn, m, mapsbc(m), 'burn', m)

      if ( m-1 .gt. numsbc) then
        print*, '==> Error: increase numsbc in csbc.h to ', m-1
        stop '=>UVic_ESCM'
      endif

      return
      end

      subroutine tracer_init

!-----------------------------------------------------------------------
!     Initialize ocean tracer names
!-----------------------------------------------------------------------

      implicit none

      integer m

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "atm.h"

      include "mw.h"

!-----------------------------------------------------------------------
!     Initialize ocean tracer names
!-----------------------------------------------------------------------

      mapt(:) = " "
      itemp = 0
      isalt = 0
      iidealage = 0
      idic = 0
      ic14 = 0
      icfc11 = 0
      icfc12 = 0
      ialk = 0
      io2 = 0
      ipo4 = 0
      ino3 = 0
      iphyt = 0
      izoop = 0
      idetr = 0

      imp = 0
      impa = 0
      impp = 0

      idetrz = 0

      idiaz = 0

      icocc = 0

      icaco3 = 0

      idetr_B = 0

      m = 1
      call set (itemp, m, mapt(m), 'temp', m)
      call set (isalt, m, mapt(m), 'salt', m)

      call set (idic, m, mapt(m), 'dic', m)

      call set (ic14, m, mapt(m), 'c14', m)

      call set (imp, m, mapt(m), 'mp', m)
      call set (impa, m, mapt(m), 'mpa', m)
      call set (impp, m, mapt(m), 'mpp', m)

      call set (ialk, m, mapt(m), 'alk', m)

      call set (io2, m, mapt(m), 'o2', m)

      call set (ipo4, m, mapt(m), 'po4', m)
      call set (iphyt, m, mapt(m), 'phyt', m)
      call set (izoop, m, mapt(m), 'zoop', m)
      call set (idetr, m, mapt(m), 'detr', m)

      call set (idetrz, m, mapt(m), 'detrz', m)

      call set (ino3, m, mapt(m), 'no3', m)
      call set (idiaz, m, mapt(m), 'diaz', m)

      call set (icocc, m, mapt(m), 'cocc', m)

      call set (icaco3, m, mapt(m), 'caco3', m)

      call set (idetr_B, m, mapt(m), 'detr_B', m)

      if ( m-1 .gt. nt) then
        print*, '==> Error: increase nt for tracers in size.h'
        stop '=>UVic_ESCM'
      endif
!-----------------------------------------------------------------------
!     Initialize ocean biology tracer names
!     Added by KK
!-----------------------------------------------------------------------

      mapb(:) = " "
      ibion = 0

      ibiop = 0
      ibioz = 0
      ibiod = 0

      ibiomp = 0.
      ibiompa = 0.
      ibiompp = 0.

      ibiodz = 0.

      ibiono3 = 0
      ibiodiaz = 0

      ibioc = 0

      ibiocaco3 = 0

      ibiod_B = 0

      m = 1

      call set (ibion, m, mapb(m), 'po4', m)

      call set (ibioc, m, mapb(m), 'cocc', m)

      call set (ibiop, m, mapb(m), 'phyt', m)
      call set (ibioz, m, mapb(m), 'zoop', m)
      call set (ibiod, m, mapb(m), 'detr', m)

      call set (ibiomp, m, mapb(m), 'mp', m)
      call set (ibiompa, m, mapb(m), 'mpa', m)
      call set (ibiompp, m, mapb(m), 'mpp', m)

      call set (ibiodz, m, mapb(m), 'detrz', m)

      call set (ibiono3, m, mapb(m), 'no3', m)
      call set (ibiodiaz, m, mapb(m), 'diaz', m)

      call set (ibiocaco3, m, mapb(m), 'caco3', m)

      call set (ibiod_B, m, mapb(m), 'detr_B', m)

      if ( m-1 .gt. ntnpzd) then
        print*, '==> Error: increase ntnpzd for tracers in size.h'
        stop '=>UVic_ESCM'
      endif

!-----------------------------------------------------------------------
!     Initialize ocean tracer source names, must have equivalent tracer
!-----------------------------------------------------------------------

      mapst(:) = " "
      itrc(:) = 0

      m = 1

      call set (isdic, m, mapst(m), 'sdic', m)
      itrc(idic) = m-1

      call set (isc14, m, mapst(m), 'sc14', m)
      itrc(ic14) = m-1

      call set (ismp, m, mapst(m), 'smp', m)
      itrc(imp) = m-1
      call set (ismpa, m, mapst(m), 'smpa', m)
      itrc(impa) = m-1
      call set (ismpp, m, mapst(m), 'smpp', m)
      itrc(impp) = m-1

      call set (isalk, m, mapst(m), 'salk', m)
      itrc(ialk) = m-1

      call set (iso2, m, mapst(m), 'so2', m)
      itrc(io2) = m-1

      call set (ispo4, m, mapst(m), 'spo4', m)
      itrc(ipo4) = m-1
      call set (isphyt, m, mapst(m), 'sphyt', m)
      itrc(iphyt) = m-1
      call set (iszoop, m, mapst(m), 'szoop', m)
      itrc(izoop) = m-1
      call set (isdetr, m, mapst(m), 'sdetr', m)
      itrc(idetr) = m-1

      call set (isdetrz, m, mapst(m), 'sdetrz', m)
      itrc(idetrz) = m-1

      call set (isno3, m, mapst(m), 'sno3', m)
      itrc(ino3) = m-1
      call set (isdiaz, m, mapst(m), 'sdiaz', m)
      itrc(idiaz) = m-1

      call set (iscocc, m, mapst(m), 'scocc', m)
      itrc(icocc) = m-1

      call set (iscaco3, m, mapst(m), 'scaco3', m)
      itrc(icaco3) = m-1

      call set (isdetr_B, m, mapst(m), 'sdetr_B', m)
      itrc(idetr_B) = m-1

      if ( m-1 .gt. nt) then
        print*, '==> Error: increase nsrc for tracer sources in size.h'
        stop '=>UVic_ESCM'
      endif

!-----------------------------------------------------------------------
!     Initialize atmosphere tracer names
!-----------------------------------------------------------------------

      mapat(:) = " "
      isat = 0
      ishum = 0
      ico2 = 0

      m = 1
      call set (isat, m, mapat(m), 'sat', m)
      call set (ishum, m, mapat(m), 'shum', m)

      if ( m-1 .gt. nat) then
        print*, '==> Error: increase nat in size.h'
        stop '=>UVic_ESCM'
      endif

      return
      end

      subroutine read_namelist

!-----------------------------------------------------------------------
!     read all model namelist variables
!-----------------------------------------------------------------------

      implicit none

      character (120) :: fname, new_file_name, logfile

      integer iotraj, i, j, k, ip, kr, jq, n, ioun, num_processors

      logical exists, annlevobc, annlev, initpt

      real dtatms, snapint, snapls, snaple, snapde, ahbkg, runstep
      real slmx, s_dm, afkph, dfkph, sfkph, zfkph, ahs, ahb, trajint
      real crops_yr

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "accel.h"
      include "calendar.h"
      include "cembm.h"
      include "coord.h"
      include "cnep.h"
      include "csbc.h"
      include "cprnts.h"
      include "diag.h"
      include "emode.h"
      include "fwa.h"
      include "hmixc.h"
      include "insolation.h"
      include "iounit.h"

      include "isopyc.h"

      include "mtlm.h"
      include "npzd.h"
      include "scalar.h"
      include "sed.h"
      include "stab.h"
      include "switch.h"
      include "tmngr.h"
      include "veg.h"
      include "vmixc.h"

!-----------------------------------------------------------------------
!     read all model namelist variables
!-----------------------------------------------------------------------

      namelist /contrl/ init, runlen, rununits, restrt, initpt
     &,                 num_processors, runstep, initbgc
      namelist /tsteps/ dtts, dtuv, dtsf, dtatm, dtatms, namix, segtim
     &,                 daylen
      namelist /riglid/ mxscan, sor, tolrsf, tolrsp, tolrfs
      namelist /mixing/ am, ah, ahbkg, ambi, ahbi, kappa_m, kappa_h
     &,                 cdbot, spnep, senep, aidif, ncon, nmix, eb
     &,                 acor, dampts, dampdz, annlev, annlevobc
      namelist /diagn/  tsiint, tsiper, tavgint, itavg, tmbint, tmbper
     &,                 itmb, stabint, zmbcint, glenint, trmbint, itrmb
     &,                 vmsfint, gyreint,igyre, extint, prxzint, trajint
     &,                 exconvint, dspint, dspper, snapint, snapls
     &,                 snaple, snapde, timavgint, timavgper, cmixint
     &,                 prlat, prslon, prelon, prsdpt, predpt
     &,                 slatxy, elatxy, slonxy, elonxy
     &,                 cflons, cflone, cflats, cflate, cfldps, cfldpe
     &,                 maxcfl, xbtint, xbtper, crossint, densityint
     &,                 fctint, tyzint, restint, tbtint, tbtper
     &,                 accel, accel_yr0
      namelist /io/     expnam, logfile, runstamp, iotavg, iotmb, iotrmb
     &,                 ioglen, iovmsf, iogyre, ioprxz, ioext, iodsp
     &,                 iotsi, iozmbc, iotraj, ioxbt, isot1, ieot1
     &,                 isot2, ieot2, jsot, jeot, ksot, keot, mrot
      namelist /ictime/ year0, month0, day0, hour0, min0, sec0, ryear
     &,                 rmonth, rday, rhour, rmin, rsec, refrun, refinit
     &,                 refuser, eqyear, eqmon, monlen, init_time
     &,                 init_time_in, init_time_out, omega
      namelist /npzd/   alpha, kw, kc, abio, bbio, cbio, k1n, nup
     &,                 gamma1, gbio0, nuz, nud0, wdd, alpha_D

     &,                 wmp0, pinp, pbc, pagg, prise, kp, ping

     &,                 wdz0, rpoo

     &,                 wd0, par, dtnpzd, redctn, ki, redptn, capr
     &,                 dcaco3, nupt0, redotn, jdiar, kzoo, geZ
     &,                 zprefP, zprefZ, zprefDet, zprefDiaz, kfe, kfe_D
     &,                 zprefMP, zprefDiat, zprefC, ro2n, o2t

     &,                 abioc, k1n_C, nuct0,nuc, alpha_C
     &,                 kfe_C

     &,                 bapr

     &,                 wc0, dissk0, wdc
     &,                 kc_c

      namelist /fwa/    isfwa1, iefwa1, isfwa2, iefwa2, jsfwa, jefwa
     &,                 mrfwa, fwaflxi, fwayri, fwayrf, fwarate
     &,                 compensate
      namelist /blmix/  Ahv, Ahh
      namelist /hlmix/  hl_depth, hl_back, hl_max
      namelist /isopyc/ slmx, ahisop, athkdf, del_dm, s_dm
      namelist /ppmix/  wndmix, fricmx, diff_cbt_back, visc_cbu_back
     &,                 visc_cbu_limit, diff_cbt_limit
      namelist /smagnl/ diff_c_back
      namelist /sed/    dtsedyr, nsedacc, weath, kmin
      namelist /embm/   rlapse, rf1, rf2, scatter, rhmax, vcsref, vcsfac
     &,                 vcsyri, aggfor_os, adiff, cropf, pastf
      namelist /carbon/ co2ccn, co2emit, co2for, co2_yr, c14_yr, dc14ccn
     &,                 c14prod

      namelist /paleo/  orbit_yr, pyear, eccen, obliq, mvelp, sealev
     &,                 sealev_yr, volcano_yr, sulph_yr, aggfor_yr
     &,                 cfcs_yr
      namelist /ice/    niats, nivts, dampice, tsno, hsno_max, ice_yr
     &,                 landice_yr, ice_calb, sno_calb
      namelist /veg/    veg_alb, veg_rl, veg_rs, veg_smd, agric_yr
     &,                 crops_yr, iagric, icrops, idesert, iice
      namelist /solar/  solarconst, solar_yr
      namelist /mtlm/   TIMESTEP, INT_VEG, VEG_EQUIL, DAY_TRIF, DAY_PHEN
     &,                 BF

!     physical constants

      rho0 = 1.035
      rho0r = c1/rho0
      grav = 980.6
      radius = 6370.0e5
      pi = atan(1.0) * 4.0

!     set defaults for namelist contrl
      init = .true.
      runlen = 365.0
      rununits = 'days'
      restrt = .true.
      initpt = .false.
      num_processors = 1
      runstep = -1.0
      initbgc = .true.

!     set defaults for namelist tsteps
      dtts = 43200.0
      dtuv = 600.0
      dtsf = 600.0
      dtatm = 43200.0
      dtatms = 43200.0
      namix = 16
      segtim = 1.0
      daylen = 86400.0

!     set defaults for namelist riglid
      mxscan = 300
      sor = 1.60
      tolrsf = 5.0e8
      tolrsp = 1.0e4
      tolrfs = 1.0e4

!     set defaults for namelist mixing
      am = 2.0e9
      ah = 2.0e7
      ahbkg = 0.
      ambi = 1.0e23
      ahbi = 5.0e22
      kappa_m = 10.0
      kappa_h = 1.0
      cdbot = 1.3e-3
      spnep = 3.0e5
      senep = 12.0e5

      aidif = 1.0

      ncon = 1
      nmix = 16
      eb = .false.
      acor = 0.0
      do n=1,nt
        dampts(n) = 50.0
        dampdz(n) = 26.575e2
      enddo
      annlev = .false.
      annlevobc = .false.

!     set defaults for namelist diagn
      tsiint = 1.0
      tsiper = 1.0
      tavgint = -36500.0
      itavg = .true.
      tmbint = -36500.0
      tmbper = 365.0
      itmb = .true.
      stabint = -36500.0
      zmbcint = -36500.0
      glenint = -36500.0
      trmbint = -36500.0
      itrmb = .true.
      vmsfint = -36500.0
      gyreint = -36500.0
      igyre = .true.
      extint = -36500.0
      prxzint = -36500.0
      trajint = 0.0
      exconvint = -36500.0
      dspint = -36500.0
      dspper = -365.0
      snapint = 0.0
      snapls = 0.0
      snaple = 0.0
      snapde = 0.0
      timavgint = 36500.0
      timavgper = 365.0
      cmixint = -36500.0
      do n=1, nlatpr
        prlat(n) = 100.0
        prslon(n) = 0.0
        prelon(n) = 0.0
        prsdpt(n) = 0.0
        predpt(n) = 6000.0e2
        if (n. le. 4) then
          prslon(n) = 180.0
          prelon(n) = 250.0
        endif
      enddo
      prlat(1) = -60.0
      prlat(2) = 0.0
      prlat(3) = 27.0
      prlat(4) = 55.0
      slatxy = -90.0
      elatxy = 90.0
      slonxy = 3.0
      elonxy = 357.0
      cflons = 0.0
      cflone = 360.0
      cflats = -90.0
      cflate = 90.0
      cfldps = 0.0
      cfldpe = 6000.0e2
      maxcfl = 3
      xbtint = -36500.0
      xbtper = -365.0
      crossint = 365000.0
      densityint = -36500.0
      fctint = -36500.0
      tyzint = -36500.0
      restint = 36500.0
      tbtint = -36500.0
      tbtper = -365.0
      accel = 1.
      accel_yr0 = 0.

!     set defaults for namelist io
      expnam = ' '
      logfile = 'machine.log'
      runstamp = ' '
      restrt = .false.
      iotavg = -1
      iotmb = -1
      iotrmb = -1
      ioglen = -1
      iovmsf = -1
      iogyre = -1
      ioprxz = -1
      ioext = -1
      iodsp = -1
      iotsi = -1
      iozmbc = -1
      iotraj = -1
      ioxbt = -1
      isot1 = 2
      ieot1 = imtm1
      isot2 = 2
      ieot2 = 1
      jsot = 2
      jeot = jmtm1
      ksot = 1
      keot = km
      mrot = 0

!     set defaults for namelist ictime
      year0 = 1
      month0 = 1
      day0 = 1
      hour0 = 0
      min0 = 0
      sec0 = 0
      ryear = 1
      rmonth = 1
      rday = 1
      rhour = 0
      rmin = 0
      rsec = 0
      refrun = .false.
      refinit = .true.
      refuser = .false.
      eqyear = .true.
      eqmon = .false.
      monlen = 30
      init_time = .false.
      init_time_in = .false.
      init_time_out = .false.
      omega = pi/43082.0

!     set defaults for namelist npzd
      alpha = 0.1  ! Initial slope P-I curve [(W/m^2)^(-1)/day]
      alpha_D = 0.38
      kw = 0.04  ! Light attenuation due to water [1/m]
      kc = 0.07  ! Light atten. by phytoplankton [1/(m*mmol/m^3)]
      ki = 5.0  ! Light attenuation through sea ice & snow

      alpha_C = 0.1  ! Initial slope P-I curve [(W/m^2)^(-1)/day]
      abioc = 1.03  ! a; Maximum growth rate parameter coccolithophore [1/day]

      k1n_C = 0.00074  ! Half sat const for N uptake (coccolithophores) [mmol/m^3]
      nuc = 0.03  ! Specific mortality rate (coccolithophores) [1/day]
      nuct0 = 0.015  ! Specific mortality rate (coccolithophores) [1/day]

      bapr = 0.05 !+-0.011 detritus to carbonate ratio [mg POC/mg PIC]

      kc_c = 0.2  ! Light atten. by calcite [1/(m*mmol/m^3)]
      wc0 = 35.0 ! constant calcite sinking speed [m/day]
      wdc = 5.0e-2 ! increase in calcite sinking speed with depth
!      Gehlen et al 2007 use a surface sink speed of 50 m/day inc w/ depth
!      dissk0 = 0.013 !dissolution rate parameter [1/day] from Gehlen et al 2007
      dissk0 = 0.
      kcal = 100. !mol m^-3 half saturation for PIC dissl from PlankTOM10

      abio = 1.0  ! a; Maximum growth rate parameter [1/day]
      bbio = 1.066  ! b
      cbio = 1.0  ! c [1/deg_C]
      k1n = 0.001  ! Half saturation constant for N uptake [mmol/m^3]
      nup = 0.03  ! Specific mortality rate (Phytoplankton) [1/day]
      nupt0 = 0.015  ! Specific mortality rate (Phytoplankton) [1/day]
      gamma1 = 0.70  ! gama1; Assimilation efficiency (zpk)
      gbio0 = 0.4  ! Maximum grazing rate at 0 deg C [1/day]
      nuz = 0.06  ! Quadratic mortality (zpk) [(mmol/m^3)^(-2)day^(-1)]
      nud0 = 0.07  ! remineralization rate [1/day]
      wd0 = 12.  ! Sinking speed of detritus at surface [m/day]

      wmp0 = 0.019 !rising plastic rate m/s
      pinp = 0.2 !portion of total plastic production into ocean
      pbc = 0.2 !portion of bottom amount returned to water colum
      pagg = 0.03 !aggregation rate
      prise = 0.02 ! portion of free plastic that rises
      ping = 1.0 ! replacement of food by plastic
      kp= 10. !half saturation constant for free plastic uptake by mar snow

      wdz0 = 24. !sinking poo rate m/d
      rpoo = 0.5 ! ratio of zoop detr going in to poo

      wdd = 6.0e-2 ! increase in sinking speed with depth
      par = 0.43  ! fraction of photosythetically active radiation
      dtnpzd = dtts  ! time step of biology [s]
      redctn = 6.625  ! C/N Redfield ratio
      redptn = 1./16.  ! P/N Redfield ratio
      capr = 0.04  ! carbonate to carbon production ratio
      dcaco3 = 650000.0  ! remineralisation depth of calcite [cm]
      redotn = 10.  ! O2/N Redfield ratio
      jdiar = 0.4  ! factor reducing the growth rate of diazotrophs
      kzoo = 0.15  ! half sat. constant for Z grazing [mmol N m^3]
      geZ =  0.40  ! Zooplankton growth efficiency
! graz prefs needed to be restructured- too confusing w/so many choices! KK
      zprefMP = 0. !zooplankton preference for microplastic
      zprefP = 0. ! Zooplankton preference for P
      zprefC = 0. ! Zooplankton preference for C
      zprefZ = 0. ! Zooplankton preference for other Z
      zprefDet = 0. ! Zooplankton preference for Detritus
      zprefDiaz = 0. ! Zooplankton preference for diazotrophs
      zprefDiat = 0. ! Zooplankton preference for diatoms
      ro2n = 800. !ratio of o to n for denitrification, default is 800
      o2t = 5. !threshold for denitrification, default is 5 mmol O2

      kmfe = 3 ! number of depth levels used for Fe limitation
      kfe = 0.12    ! Half saturation constant for Fe limitation
      kfe_D = 0.1 ! Half saturation constant for Diaz Fe limitation

      kfe_C = 0.06

!     set defaults for namelist fwa
      isfwa1 = 1
      iefwa1 = imt
      isfwa2 = 0
      iefwa2 = 0
      jsfwa = 1
      jefwa = jmt
      mrfwa = 1
      fwaflxi = 0.
      fwayri = -1.e20
      fwayrf = 1.e20
      fwarate = 0.
      compensate = .false.

!     set defaults for namelist blmix
!     Reference:
!     A Water Mass Model of the World Ocean  K. Bryan, L.J. Lewis
!     JGR, vol 84, No. C5, May 20, 1979
      afkph = 0.8
      dfkph = 1.05
      sfkph = 4.5e-5
      zfkph = 2500.0e2
!     Use Bryan & Lewis values for vertical tracer diffusion
!     Ahv range of 0.3 to 1.3, crossover at 2500m.
!     compute depth dependent vertical diffusion coefficients for
!     tracers using the relationship of Bryan and Lewis
      pi = 4.0 * atan(1.0)

!     compute depth dependent horizontal diffusion coefficients for
!     tracers using the relationship of Bryan and Lewis
      ahs = 5.0e+3
      ahb = 1.0e+3

!     set defaults for namelist hlmix
      hl_depth = 500.0e2
      hl_back  = 1.e4
      hl_max   = 1.e9

!     set defaults for namelist isopyc
      slmx  = 1.0/100.0   ! maximum isopycnal slope
      ahisop = 1.e7       ! isopycnal diffusion coefficient
      athkdf = 1.0e7      ! isopycnal thickness diffusion coefficient
      del_dm = 4.0/1000.0 ! transition for scaling diffusion coefficient
      s_dm = 1.0/1000.0   ! half width scaling for diffusion coefficient

!     set defaults for namelist ppmix
      wndmix    = 10.0
      fricmx    = 50.0
      diff_cbt_back =  0.1
      visc_cbu_back =  1.0

!     in regions of gravitational instability set mixing limits to the
!     maximum consistant with the "cfl" criterion. convective adjustment
!     will also act on the instability.
      visc_cbu_limit = fricmx
      diff_cbt_limit = fricmx

!     set defaults for namelist smagnl
      diff_c_back = c0

!     set defaults for namelist sed
      dtsedyr = 1.   ! time step of sediment model in years
      nsedacc = 1    ! number of steps for accelerating sediments
      weath = 2.e20  ! weathering flux in Pg year-1
! set high to check for constant namelist input
      kmin = 8       ! minimum ocean level for sediments

!     set defaults namelist embm
      rlapse     = 5.e-5
      rf1        = 0.3
      rf2        = 3.e5
      scatter    = 0.23
      rhmax      = 0.85
      vcsref     = 0.
      vcsfac     = 0.
      vcsyri     = 0.
      aggfor_os  = 0.
      adiff      = 0.03
      cropf      = 1.
      pastf      = 0.5

!     set defaults namelist carbon
      co2ccn     = 280.
      co2emit    = 0.
      co2for     = 5.35e03
      co2_yr     = 2.e20
      c14_yr     = 2.e20
      dc14ccn    = 0.
      c14prod    = 0.

!     set defaults namelist paleo
      pyear      = 2.e20
      orbit_yr   = 2.e20
!     default user orbit is for 1950
      eccen      = 1.672393E-02
      obliq      = 23.44627
      mvelp      = 102.0391
      sealev     = 0.
      sealev_yr  = 2.e20
      volcano_yr = 2.e20
      sulph_yr   = 2.e20
      aggfor_yr  = 2.e20
      cfcs_yr    = 2.e20

!     set defaults namelist ice
      niats      = 1
      nivts      = 1
      dampice    = 5.
      tsno       = 0.
      hsno_max   = 1000.

      ice_yr     = 2.e20
      landice_yr = 2.e20
      ice_calb   = 0.25
      sno_calb   = 0.2

!     set defaults namelist veg
      veg_alb(1:7) = (/0.17, 0.17, 0.22, 0.22, 0.22, 0.30, 0.60/)
      veg_rl(1:7) = (/1.5, 1.0, 0.1, 0.3, 0.01, 0.01, 0.005/)
      veg_rs(1:7) = (/0.75, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
      veg_smd(1:7) = (/2.0, 3.0, 0.1, 0.1, 0.5, 0.01, 0.01/)
      agric_yr = 2.e20
      crops_yr = 2.e20
      iagric = 3
      icrops = 3
      idesert = 6
      iice = 7

!     set defaults namelist solar
      solarconst = 1.368e6  ! mw m-2
      solar_yr = 2.e20

!     set defaults namelist mtlm
      TIMESTEP = 3600.
      INT_VEG = .true.
      VEG_EQUIL = .false.
      DAY_TRIF = 30
      DAY_PHEN = 1
      BF = 1.

!-----------------------------------------------------------------------
!     provide for change in above presets using "namelist"
!-----------------------------------------------------------------------

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, contrl, end=101)
101   continue
      runstep = float(int(runstep))
      write (stdout, contrl)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, tsteps, end=102)
102   continue
      write (stdout, tsteps)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, riglid, end=103)
103   continue
      write (stdout, riglid)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, mixing, end=104)
104   continue
      write (stdout, mixing)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, diagn, end=105)
105   continue
      write (stdout, diagn)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, io, end=106)
106   continue
      write (stdout, io)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, ictime, end=107)
107   continue
      write (stdout, ictime)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, npzd, end=108)
108   continue
      write (stdout, npzd)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, fwa, end=109)
109   continue
      write (stdout, fwa)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, blmix, end=110)
110   continue
      write (stdout, blmix)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read (ioun, hlmix, end=111)
      write (stdout, hlmix)
      call relunit (ioun)
111   continue

      call getunit (ioun, 'control.in', 'f s r')
      read (ioun, isopyc, end=112)
112   continue
      write (stdout, isopyc)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read (ioun, ppmix, end=113)
113   continue
      write (stdout, ppmix)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read (ioun, smagnl, end=114)
114   continue
      write (stdout, smagnl)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, sed, end=115)
115   continue
      write (stdout, sed)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, embm, end=116)
116   continue
      write (stdout, embm)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, carbon, end=117)
117   continue
      write (stdout, carbon)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, paleo, end=118)
118   continue
      write (stdout, paleo)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, ice, end=119)
119   continue
      write (stdout, ice)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, veg, end=120)
120   continue
      write (stdout, veg)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, solar, end=121)
121   continue
      write (stdout, solar)
      call relunit (ioun)

      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, mtlm, end=122)
122   continue
      write (stdout, mtlm)
      call relunit (ioun)

!     limit or convert some variables

!     ensure pass is between zero and one.
      pass =  min(max((1. - scatter), 0.), 1.)

      dtocn = dtts
      nots = nmix
      isot1 = max(isot1, 1)
      ieot1 = min(ieot1, imt)
      isot2 = max(isot1, 1)
      ieot2 = min(ieot1, imt)
      jsot = max(jsot, 1)
      jeot = min(jeot, jmt)
      ksot = max(ksot, 1)
      keot = min(keot, km)
      do k=1,km
        dtxcel(k) = 1.0
      enddo

      ah = ahbkg
      slmxr = 0.
      if (slmx .gt. 0.) slmxr = c1/slmx
      s_dm = 0.
      if (s_dm .gt. 0.) s_dmr = c1/s_dm

!     sort out forcing years
      if (pyear  .gt. 1.e20) pyear = 1800.

      if (c14_yr .gt. 1.e20) c14_yr = pyear

      if (orbit_yr .gt. 1.e20) orbit_yr = pyear

      if (volcano_yr .gt. 1.e20) volcano_yr = pyear

      if (sulph_yr .gt. 1.e20) sulph_yr = pyear

      if (aggfor_yr .gt. 1.e20) aggfor_yr = pyear

      if (cfcs_yr .gt. 1.e20) cfcs_yr = pyear

      if (sealev_yr .gt. 1.e20) sealev_yr = pyear

      if (agric_yr .gt. 1.e20 .and. crops_yr .lt. 1.e20)
     &  agric_yr = crops_yr

      if (landice_yr .gt. 1.e20 .and. ice_yr .lt. 1.e20)
     &  landice_yr = ice_yr

      ice_yr = landice_yr

      tsiper = min(tsiper, tsiint)
      tbtper = min(tbtper, tbtint)
      if (accel .le. 1) accel = 1.
      runlen = runlen/accel
      runstep = runstep/accel
      if (runstep .gt. 0.0) runlen = min(runstep, runlen)
      if (init_time) then
        init_time_in = .true.
        init_time_out = .true.
      endif

      monlen = 30

      eqyear = .true.
      eqmon = .false.
      calendar = 'noleap'

!-----------------------------------------------------------------------
!     set runstamp
!-----------------------------------------------------------------------
      if (runstamp .eq. ' ') then
!       read runstamp from last line of logfile if it exists
        fname = new_file_name (logfile)
        inquire (file=trim(fname), exist=exists)
        if (exists) then
          call getunit (ioun, fname, 'f s r')
          do while (exists)
            read  (ioun, '(a130)', end=130) runstamp
          enddo
130       continue
          call relunit (ioun)
          print*, " "
          print*, "runstamp: ", trim(runstamp)
        endif
      endif

      return
      end