      subroutine mom
c
c=======================================================================
c
c                        GFDL Modular Ocean Model
c
c     A primitive equation ocean model developed by researchers at  
c         the Geophysical Fluid Dynamics Laboratory /NOAA.
c                                
c     For details, refer to the manual.
c
c     author:  r.c.pacanowski   e-mail  rcp@gfdl.gov
c    (mpp items)  v. balaji     e-mail=> vb@gfdl.gov
c=======================================================================
c
      use arrays_2d
      use datamod
      use memory_window
      use asm_x_mod

      use mom3_mpp_mod
c
      external blkdta
c
      include 'mpif.h'
!====================== include file "param.h" =========================
!
!====================== include file "size.h" ==========================
!                                                                       
!-----------------------------------------------------------------------
!     USER INPUT: 
!-----------------------------------------------------------------------
!
!     imt    = number of grid cells in the longitudinal direction
!              (calculated cells are from 2 through imt-1. end cells
!               are boundaries)          
!
!     jmt    = number of grid cells (latitude rows) in the latitudinal
!              direction (calculated cells are from 2 through jmt-1.
!              end cells are boundaries)   
!                                
!     km     = number of grid cells in the vertical direction          
!              (calculated cells are from 1 through km. no boundary
!               cells)
!
!     nt     = number of tracers (temperature, salinity, ...)
!
!     mnisle = maximum number of islands (unconnected land masses)
!
!     maxipp = maximum number of all island perimeter points
!-----------------------------------------------------------------------
!                                                                       

























! window size must be appropriate for options with highest order numerics




! specify lomem or himem on 3D diagnostics. default is himem

!these are used to determine hardware/OS/compiler (V. Balaji vb@gfdl.gov)



!compilers that support Cray pointers

!values of kind: double and long are 8-byte, float and int are 4-byte
!these might be different on non-SGICRAY, I believe


!parallel machine types

      parameter (imt= 362, jmt=  202, km= 40)
      parameter (nt=2)
      parameter (mnisle=50, maxipp=10000)
      parameter (kbot = km)
!
!
!-----------------------------------------------------------------------
!     END USER INPUT
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!     "jmw"  = number of latitude rows within the MW (memory window).
!     "jbuf" = number of northermost and southernmost buffer rows.
!              The MW is symmetric because it has the same number of
!              northermost and southernmost buffer rows "jbuf".
!              total number of buffer rows = 2*jbuf
!-----------------------------------------------------------------------
!
      parameter (jbuf=2)
      parameter (jmw=1 + 2*jbuf)
!
!-----------------------------------------------------------------------
!      num_rows_to_copy = number of northernmost rows to copy into
!                         the southernmost rows of the MW when moving
!                         the MW northward
!-----------------------------------------------------------------------
!
      parameter (num_rows_to_copy = 2*jbuf)
!
!     nvar   = number of prognostic variables
!     lseg   = maximum number of longitudinal stream function segments
!     nlatpr = maximum number of latitudes for matrix printouts
!              on diagnostic time steps
!     nhreg  = number of regions in the horizontal used for averaging
!              tracers.
!     nvreg  = number of regions in the vertical used for term balance
!              calculations. note "nvreg" isn`t used for tracer averages
!     numreg = total number of regions ( = product of nhreg & nvreg)
!              used for term balance calculations
!
      parameter (lseg=5, nlatpr=10)
      parameter (nhreg=4, nvreg=6, numreg=nhreg*nvreg)
      parameter (nvar=nt+2)
!
      parameter (imtm1=imt-1, imtm2=imt-2)
      parameter (kmm1=km-1, kmp1=km+1)                                
      parameter (jmtm1=jmt-1, jmtm2=jmt-2)
      parameter (imtkm=imt*km, nwds=imt*jmt, nkflds=2)
      parameter (nslab=imt*nvar*km, ntmin2=nt+1/nt)
!
!     parameter constants
!
!====================== include file "pconst.h" ========================
!
!
!     rules for parameter constants
!
!     use prefix of "c" for whole real numbers (eg: c57 for 57.0)
!     use "m" after prefix to designate negative values (minus sign)
!       (eg: cm7 for -7.0)
!     use prefix of "p" for non repeating fractions (eg: p5 for 0.5)
!     use prefix of "r" for reciprocals (eg: r3 for 1/3.0)
!     combine use of prefix above and "e" for scientific notation, with
!       (eg: c5e4 for 5.0e4, c1em10 for 1.0e-10)
!
      parameter (c0=0.0, c1=1.0, c2=2.0, c3=3.0, c4=4.0, c5=5.0, c7=7.0)
      parameter (c8=8.0)
      parameter (c14=14.0, c16=16.0, c360=360.0)
      parameter (p125=0.125, p25=0.25, p5=0.5, p75=0.75)
      parameter (epsln=1.0e-20)
!
      parameter (c24=24.0, c60=60.0, c1440=1440.0)
      parameter (r24=c1/c24, r60=c1/c60, r1440=c1/c1440)
      parameter (secday=c1/(c60*c1440))

!     minimum allowable thickness (cm) of the top cell when using the
!     nonlinear free surface.        
      parameter (thickmin=10.0)      
!
!     standard I/O units
!
!====================== include file "stdunits.h" ======================
!
!     stdin  = unit number for standard input.
!     stdout = unit number for standard output.
!     stderr = unit number for standard error.
!
      integer stdin, stdout, stderr
      parameter (stdin = 5, stderr = 0)
      common /stdunit/ stdout
c====================== include file "emode.h" =========================
c
c
c     mxscan  = max number of allowable scans for poisson solvers
c     mscan   = actual number of scans taken by poisson solvers
c     tolrsf  = tolerance for stream function calculation. 
c               the solution is halted when it is within "tolrsf"
c               of the "true" solution assuming geometric convergence.
c     tolrsp  = tolerance for surface pressure calculation
c               the solution is halted when it is within "tolrsp"
c               of the "true" solution assuming geometric convergence.
c     tolrfs  = tolerance for implicit free surface calculation
c               the solution is halted when it is within "tolrfs"
c               of the "true" solution assuming geometric convergence.
c     esterr  = estimated maximum error in elliptic solver assuming
c               geometric convergence
c    
c
c     nisle = number of land masses
c     nippts= number of land mass perimeter points
c     iperm = "i" coordinate for the land mass perimeter point
c     jperm = "j" coordinate for the land mass perimeter point
c     iofs  = offset for indexing into the land mass perimeter points
c     imask = controls whether calculations get done on perimeters
c     set mask for land mass perimeters on which to perform calculations
c     imask(-n) = .false.  [no equations ever on dry land mass n]
c     imask(0)  = .true.   [equations at all mid ocean points]
c     imask(n)  = .true./.false [controls whether there will be
c                                equations on the ocean perimeter of
c                                land mass n]
c     note: land mass 1 is the northwest-most land mass
c     for the numbering of the other landmasses, see generated map(i,j)
      logical imask
      common /emodei/ nippts(mnisle), iofs(mnisle), iperm(maxipp)
      common /emodei/ jperm(maxipp), nisle, imain
      common /emodel/ imask (-mnisle:mnisle)
c
      character*16 variable
      logical converged
      common /emodel/ converged
      common /emodec/ variable
c
c====================== include file "iounit.h" ========================
c
c     i/o units and related variables
c
c     taum1disk = disk pointer for tau-1 latitude rows
c     taudisk   = disk pointer for tau   latitude rows
c     taup1disk = disk pointer for tau+1 latitude rows
c     kflds     = disk unit used for two dimensional fields
c     latdisk   = disk units for latitude rows (alternately pointed to 
c                by taum1disk, taudisk, and taup1disk)
c              
c     iodoc  = unit for documentation
c     iostab = unit for stability testing
c     iotim  = unit for time means
c     iotim1 = scratch disk (SSD) unit for accumulating time means
c     ionew1 = unit for reading sponge layer data
c     ionew2 = mirror unit of sponge layer data on SSD
c     
c     for the following, a control # < 0 implies that unformatted data
c     will be written to a unit selected by the i/o manager "iomngr.F"
c     and given a hardwired name (grep getunit *.F to see names) 
c     and formatted data (to stdout) will be written. if a # > 0 and
c      # <> stdout, only unformatted data will be written.
c
c     iotavg = control # for tracer averages
c     iotmb  = control # for writing tracer meridional budget.
c     iotrmb = control # for term balances for tracer and momentum
c     ioglen = control # for writing global energetics integrals
c     iovmsf = control # for writing meridional stream function
c     iogyre = control # for writing gyre transport.
c     ioprxz = control # for writing x-z sections from latitudes
c     ioext  = control # for writing external mode (stream function)
c     iodsp  = control # for writing diagnostic surface pressure
c     iotsi  = control # for writing time step integrals
c     ioxbt  = control # for writing time averaged xbt data
c     iozmbc = control # for writing zonal mean surf boundary conditions
c
      integer taum1disk, taudisk, taup1disk
      common /iounit/ taum1disk, taudisk, taup1disk
      common /iounit/ latdisk(3), kflds
      common /iounit/ iodoc, iostab, iotavg, iotmb, iotrmb
      common /iounit/ iotim, iotim1
      common /iounit/ ioglen, iovmsf, iogyre, ioprxz, ioext, iodsp
      common /iounit/ iotsi, iozmbc, ionew1, ionew2, ioxbt
c
c     iotext  = 80 character text string for describing the details
c                of the next unformatted data record.
c     expnam  = 60 character text string for the experiment name
c
      character iotext*80, expnam*60
      common /iounitc/ iotext, expnam
c
c     when writing unformatted data records in MOM, each data record is
c     preceded by a header record which was written as:
c     write(unit) stamp, iotext, expnam
c     where stamp is a 32 character specification of the model date &
c     time corresponding to the time step when the data was written and
c     iotext is a 80 character description of what is in the
c     data record and how it is to be read. expnam is a 60 character
c     experiment name which shows which experiment wrote the data.
c     this makes it easy to decipher any unformatted output from the
c     model by using a program similar to the following: 
c
c      program decifr
c
c-----------------------------------------------------------------------
c      decifer an unformatted file from MOM by showing the header
c      records. the file needs to copied to file "fort.21"
c-----------------------------------------------------------------------
c
c      character*32 stamp
c      character*80 iotext
c      character*60 expnam
c
c      iounit = 21
c      rewind iounit
c      do n=1,100000
c
c        read the header record
c
c        read (iounit, end=110) stamp, iotext, expnam
c        write (*,'(1x,a32,1x,a80)') stamp, iotext
c
c        skip the data record
c
c        read (iounit)
c      enddo
c110   continue
c      write (*,*) " => end of file on fort.",iounit
c      stop
c      end
c
c     note: all unformatted diagnostic MOM data is handled this way.
c     to insure that data is read properly, verify that arrays are
c     dimensioned correctly by comparing the listed variables against
c     those in the *.h files. (grep -i -n "variable" *.h) Also, most
c     data from MOM is written IEEE 32bit so it is read directly by
c     most workstations. However, when trying to read these IEEE files
c     on the CRAY, they must be assigned IEEE before being read.
c     Some diagnostic data is averaged over time before being written.
c     In these cases, the time "stamp" refers to the last time step
c     at the end of the averaging period. An averaging interval is
c     also written as part of the data. Averaging periods = zero 
c     indicate instantaneous data.
c
c======================= include file "mw.h" ===========================
c
c     Refer to the manual for a description of the memory window (MW)
c
c     author:  r.c.pacanowski   e-mail  rcp@gfdl.gov
c=======================================================================
c
c     taum1 = tau-1 time level for variables in MW
c     tau   = tau   time level for variables in MW
c     taup1 = tau+1 time level for variables in MW
c
      integer taum1, tau, taup1
      common /mwi/ taum1, tau, taup1
c
c-----------------------------------------------------------------------
c     MW arrays for prognostic equations:
c-----------------------------------------------------------------------
c
c     u(i,k,j,n,tau) = total velocity where:
c      i   = index for longitude
c      k   = index for depth
c      j   = index for latitude row within MW
c      n   = component (1 = zonal, 2 = meridional)
c      tau = time level (tau-1, tau, tau+1)
c             
c      (only internal modes are on disk and at tau+1 in the MW)
c
c     t(i,k,j,n,tau) = tracer where:
c      i   = index for longitude
c      k   = index for depth
c      j   = index for latitude row within MW
c      n   = component (1 = temperature, 2 = salinity)
c            if nt > 2 then other tracers are allowed.
c      tau = time level (tau-1, tau, tau+1)
c
c     note: temperature is potential temperature in degrees Celsius and 
c           salinity is in "model units", the deviation from 0.035 grams
c           of salt/cm**3 of water, or, assuming a water density of
c           1 gram/cm**3, the deviation from 0.035 g of salt/g of water.
c           one can convert model units to the common units of parts
c           per thousand (ppt) by adding 0.035 grams/cm**3 to the model
c           units and then multiplying by 1000. 
c     
      common /mw/ u(imt,km,jmw,2,-1:1), t(imt,km,jmw,nt,-1:1)
c
c-----------------------------------------------------------------------
c     MW arrays for diagnostic equations and workspace:
c-----------------------------------------------------------------------
c
c     diagnostic advective velocities are in units of cm/sec
c
c     adv_vet  = advective velocity on the eastern face of a T cell
c     adv_vnt  = advective velocity on the northern face of a T cell
c     adv_veu  = advective velocity on the eastern face of a U cell
c     adv_vnu  = advective velocity on the northern face of a U cell
c     adv_vbt  = advective velocity on the bottom face of a T cell
c     adv_vbu  = advective velocity on the bottom face of a U cell
c
c     rho      = density at center of a T cell in units of gm/cm**3
c                note: there is an arbitrary constant which is only a
c                function of depth in "rho". It is related to 
c                subtracting a reference level density for purposes of
c                accuracy.
c
c     pressure = hydrostatic pressure at the grid point within a T cell
c     grad_p   = hydrostatic pressure gradient for U cell. There are
c                two components: (1,2) is for (dp/dx, dp/dy)
c
      common /mw/ adv_vet(imt,km,2:jmw), adv_vnt(imt,km,1:jmw)
      common /mw/ adv_veu(imt,km,2:jmw-1)
      common /mw/ adv_vnu(imt,km,1:jmw-1)
      common /mw/ adv_vbt(imt,0:km,2:jmw)
      common /mw/ adv_vbu(imt,0:km,2:jmw-1)
      common /mw/ rho(imt,km,1:jmw)
      common /mw/ pressure(imt,km,1:jmw)
      common /mw/ grad_p(imt,km,1:jmw-1,2)
c
c     tmask = tracer cell land/sea mask   = (0.0, 1.0) on (land, sea)
c     umask = velocity cell land/sea mask = (0.0, 1.0) on (land, sea)
c
      common /mw/ tmask(imt,km,1:jmw), umask(imt,km,1:jmw)
c
      common /mw/ dht(imt,km,jmw+1), dhu(imt,km,jmw)
      common /mw/ dhwt(imt,0:km,jmw)
      common /mw/ pc_sink(imt,km,2:jmw-1)
      common /mw/ ztp(imt,km,jmw)
c
c
c     adv_f4n = 4th order advective flux
c               
      common /mw/ adv_f4n(imt,km,1:jmw-1,nt)
c
c
c
c     these workspace arrays are recalculated for each component of the
c     equations so do not have to be moved as the MW moves northward. 
c
c     adv_fe   = advective flux across the eastern face of a cell
c     adv_fn   = advective flux across the northern face of a cell
c                (removed in most cases and put directly into the 
c                 statement functions for speed optimization.)
c     adv_fb   = advective flux across the bottom face of a cell
c
c     diff_fe  = diffusive flux across the eastern face of a cell
c     diff_fn  = diffusive flux across the northern face of a cell
c     diff_fb  = diffusive flux across the bottom face of a cell
c     source   = source term
c
      common /mw/ adv_fe(imt,km,2:jmw-1)
      common /mw/ adv_fn(imt,km,1:jmw-1)
      common /mw/ adv_fb(imt,0:km,2:jmw-1)
c
      common /mw/ diff_fe(imt,km,2:jmw-1)
      common /mw/ diff_fn(imt,km,1:jmw-1)
      common /mw/ diff_fb(imt,0:km,2:jmw-1)
      common /mw/ diff_fbiso(imt,0:km,2:jmw-1)
c
      common /mw/ source(imt,km,2:jmw-1)
c
      common /mw/ zzi(imt,km,2:jmw-1)
c
c     these grid factors are for optimizations
c     
      common /mw/ cstdxtr(imt,jmw), cstdxur(imt,jmw)
      common /mw/ cstdxt2r(imt,jmw),ah_cstdxur(imt,jmw)
      common /mw/ csudxtr(imt,jmw), csudxur(imt,jmw)
      common /mw/ csudxu2r(imt,jmw),am_csudxtr(imt,jmw)
c
c     advmet = coeff for metric advection.
c
      common /advec/ advmet(jmt,2)
c
c
c
c     smf = surface momentum flux
c          1 => zonal wind stress (dynes/cm**2)
c          2 => meridional wind stress (dynes/cm**2)
c     bmf = bottom momentum flux
c          1 => zonal bottom drag (dynes/cm**2)
c          2 => meridional bottom drag (dynes/cm**2)
c     stf = surface tracer flux
c          1 => surface heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec)
c          2 => surface salt flux (grams of salt/cm**2/sec)
c     stf_turb 
!         = the turbulent part of surface tracer flux
c     btf = bottom tracer flux (for consistency but normally zero!)
c          1 => bottom heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec)
c          2 => bottom salt flux (grams of salt/cm**2/sec)
c
      common /mw/ smf(imt,1:jmw,2),  bmf(imt,1:jmw,2)
      common /mw/ stf(imt,1:jmw,nt), btf(imt,1:jmw,nt)
      common /mw/ stf_turb(imt,1:jmw,nt)
c
c
c====================== include file "scalar.h" ========================
c
c     various scalar quantities:
c
c     dtts   = time step for density & tracers (in seconds)
c     dtuv   = time step for baroclinic velocity (in seconds)
c     dtsf   = time step for barotropic velocity (in seconds)
c     c2dtts = 2*dtts
c     c2dtuv = 2*dtuv
c     c2dtsf = 2*dtsf
c     smooth = smoothing coefficient for robert time filter
c     acor   = (>0, 0) = (implicit, explicit) treatment of coriolis
c               term for internal and external modes.
c     rho0   = mean density for Boussinesq approximation
c     rho0r  = 1/rho0 
c     omega  = earth`s rotation rate (radians/sec)
c     radius = earth`s radius (cm)
c     grav   = earth`s gravitational acceleration (cm/sec**2)
c     cdbot  = bottom drag coefficient
c     ncon   = number of  passes through convective code in tracer
c     gcor   = time centering for coriolis term
c
c     taux0  = constant zonal windstress (dynes/cm**2) for idealized
c              equatorial studies
c     tauy0  = constant meridional windstress (dynes/cm**2) for
c              idealized equatorial studies
c
      common /scalar/ dtts, dtuv, dtsf, c2dtts, c2dtuv, c2dtsf, acor
      common /scalar/ smooth
      common /scalar/ rho0, rho0r, omega, radius, grav, cdbot, gcor
      common /scalar/ taux0, tauy0
      common /scalri/ ncon
c
c     bring in non dimensional constants
c
c====================== include file "ndcon.h" ========================
c
c     various non dimensional quantities:
c
c     radian = degrees per radian
c     pi     = something good to eat
c
      common /ndcon/ radian, pi
c
c
!====================== include file "switch.h" ========================
!
!     all time dependent decisions are made by time manager "tmngr.F"
!     and communicated elsewhere to the model via logical switches.
!
!     inputs: (defaulted in "blkdta.F", optionally reset via namelist)
!
!
!     runlen  = integration period (see rununits). note "runlen" should
!               be an integral number of density time steps. if not,
!               then "runlen" is automatically adjusted to insure this.
!               fractional days are supported but not fractional months
!               or years.
!     rununits= units of "runlen". may be "days", "months", or "years".
!               tmngr will convert "runlen" which is in "rununits"
!               to "rundays" in units of days.
!
!     segtim  = the integration time "runlen" is broken into a number of
!               segments each of length "segtim" days. updated surface
!               boundary conditions are applied to MOM every "segtim" 
!               days. this is useful when coupling to atmospheric models
!               in which case both models exchange surface boundary
!               conditions every "segtim" days where "segtim"
!               is 1/(coupling frequency). without an atmospheric model,
!               when getting surface boundary conditions from data, 
!               "segtim" is set to the time step (in days) by mom.F. in
!               either case, "runlen" (in days) should be an integral
!               number of "segtim".     
!
!     nmix    = number of time steps between mixing timesteps. used
!               to damp timestep splitting due to centered leapfrog.
!
!     init    = (true,false)  indicates that this run is a
!               (start from initial conditions, restart)
!
!     restrt  = (true,false) = (do,don`t) write a restart at the end
!               of the run
!
!     eb      = (true,false) configures for the use of a 
!               (euler backward,forward) type mixing timestep
!
!-----------------------------------------------------------------------
!     inputs to tmngr.F: diagnostic intervals
!-----------------------------------------------------------------------
!
!     note: switches are used to control the interval between doing
!           diagnostics. units for all switches are in days. 
!           setting a switch < 0.0 disables its diagnostic 
!           setting a switch = 0.0 causes the diagnostic to be 
!           done every time step, and setting a switch > 0.0 causes the
!           diagnostic to be done repeatedly at the specified interval
!           from a specified starting date. 
!
!     cmixint = number of days between writing estimated mixing coeffs
!               on faces of T cells and U cells
!
!     crossint = number of days between writing diapycnal and isopycnal
!               components of flow
!
!     fctint = number of days between writing difference between 
!              FCT and leapfrog advection
!
!     pressint = number of days between writing pressure field 
!
!     exconvint = number of days between writing temperature rate of
!                 change due to explicit convection
!
!     glenint =  number of days between global energetics integrals.
!
!     trmbint =  number of days between momentum and tracer term
!                balances (global and regional).
!
!     itrmb   = (true,false) = (do,don`t) write regional mask info for
!               the term balance diagnostic. Typically set true
!               at the beginning of a run; otherwise false since it is
!               not necessary to keep writing a time independent field
!               particularly when it may be a significant part of the
!               time dependent part of the diagnostic.
!
!     gyreint =  number of days between calculation of tracer northward
!                transport.
!     igyre   = (true,false) = (do,don`t) write regional mask info for
!               the gyre diagnostic. Typically set true
!               at the beginning of a run; otherwise false since it is
!               not necessary to keep writing a time independent field
!               particularly when it may be a significant part of the
!               time dependent part of the diagnostic.
!
!     vmsfint =  number of days between calculation of vertical and
!                meridional stream function.
!
!     tyzint  =  number of days between calculation of zonally averaged
!                tracer components.
!
!     rhoint  = number of days between writing instantaneous terms 
!               in the locally referenced potential density equation.
!               Relevant only when 
!               "save_density_terms"
!               is enabled. see "iounit.h" for more details.
!               all data is at time level "tau" so it is one time step
!               before the date shown in "stamp".
!
!     diffint   = number of days between computing nonconstant 
!               diffusivity.
!
!     diffsnapint = number of days between saving nonconstant 
!               diffusivity information to netcdf file.
!
!     prxzint =  number of days between printouts of x-z data.
!
!     extint  =  number of days between printouts of external mode.
!
!     dspint  =  number of days between surface pressure calculation.
!                Note: only when "diagnostic_surface_height" is enabled.
!     dspper  = averaging period for "diagnostic_surface_height"
!
!     tavgint = number of days between regional tracer averages (under
!               horizontal regions).
!
!     itavg   = (true,false) = (do,don`t) write regional mask info for
!               the tracer average diagnostic. Typically set true
!               at the beginning of a run; otherwise false since it is
!               not necessary to keep writing a time independent field
!               particularly when it may be a significant part of the
!               time dependent part of the diagnostic.
!
!     tmbint  = number of days over which tracer equation in averaged
!               in depth and longitude to determine the meridional 
!               balance among storage, divergence, dissipation and 
!               forcing.
!     tmbper  = averaging period for "meridional_tracer_balance"
!
!     itmb    = (true,false) = (do,don`t) write "msktmb" for tracer
!               the meridional balance diagnostic. Typically set true
!               at the beginning of a run; otherwise false since it is
!               not necessary to keep writing a time independent field
!               particularly when it may be a significant part of the
!               time dependent part of the diagnostic.
!
!
!     tsiint  = number of days between printing of time step integrals.
!
!     stabint = number of days between sampling for various stability
!               criteria.
!
!     snapint = number of days between writing instantaneous samples
!               of data from the MOM grid. note: only when "snapshots"
!               is enabled. see "iounit.h" for more details.
!               all data is at time level "tau" so it is one time step
!               before the date shown in "stamp".
!
!     timavgint= interval (days) for writing time mean data
!               (only when "time_averages" is enabled). 
!               if "timavgint" is not an integral number of
!               density time steps,"timavgint" is automatically adjusted
!               to insure this. if the number of days to integrate is
!               not an integral number of "timavgint" then the last
!               averaging period will be less than "timavgint" days.this
!               may lead to one more averaging period than expected. 
!               see "iounit.h" for more details.
!     timavgper= averaging period for "time_averages"
!
!     xbtint  = averaging period (days) for writing XBT data (only when
!               "xbts" is enabled). if "xbtint" is not an integral
!               number of density time steps, "xbtint" is automatically
!               adjusted to insure this. if the number of days to 
!               integrate is not an integral number of "xbtint" then the
!               last averaging period will be less than "xbtint" days.
!               this may lead to one more averaging period than  
!               expected. see "iounit.h" for more details.
!     xbtper  = averaging period for "xbts"
!
!     zmbcint = number of days between calculation of zonal mean
!               surface boundary conditions (and related  quantities)
!
!     trajint = number of days between writing particle trajectories.
!               particle trajectories are always written on the last
!               time step of the run to allow for restarting.
!
!-----------------------------------------------------------------------
!     outputs from tmngr.F: logical switches
!-----------------------------------------------------------------------
!
!     rundays = integration time in days (from "runlen")
!
!     the following are logical counterparts to the above switches and
!     are set every time step. logical switches control all decisions
!     about when to do things in MOM.
!
!     cmixts  = (false,true) = (don`t, do) do write estimated mixing
!               coefficients on this time step.
!               based on "cmixint".
!
!     crossts  = (false,true) = (don`t, do) write diapycnal and 
!               isopycnal components of flow on this time step.
!               based on "crossint".
!
!     fctts    = (false,true) = (don`t, do) write difference between 
!               FCT and leapfrog advection on this time step.
!               based on "fctint".
!
!     pressts  = (false,true) = (don`t, do) write pressure on this time
!               step. based on "pressint".
!               
!
!     exconvts  = (false,true) = (don`t, do) do write temperature change
!               due to explicit convection on this time step.
!               based on "exconvint".
!
!     glents  = (false,true) = (don`t, do) do calculation of global
!               energy integrals on this time step. based on "glenint".
!
!     trmbts  = (false,true) = (don`t, do) do calculation of momentum &
!               tracer term balance on this timestep. based on "trmbint"
!
!     gyrets  = (false,true) = (don`t, do) do calculation of tracer
!               northward transport on this timestep. based on "gyreint"
!
!     vmsfts  = (false,true) = (don`t, do) do calculation of vertical
!               and meridional stream function on this time step.
!               based on "vmsfint"
!
!     tyzts   = (false,true) = (don`t, do) do calculation of zonally
!               averaged tracer components on this time step.
!               based on "tyzint"
!
!     prxzts  = (false,true) = (don`t, do) do printouts of x-z data
!               on this time step. based on "prxzint"
!
!     extts  = (false,true) = (don`t, do) do printout of external mode
!               on this time step. based on "extint"
!
!     dspts  = (false,true) = (don`t, do) do calculation of diagnosti! 
!              surface pressure on this time step. based on "dspint"
!               
!     rhots  = (false,true) =(don`t, do) do calculation of diagnostic
!              for locally referenced potential density terms.
!              Based on "rhoint".
!
!     diffts    = (false,true) =(don`t, do) do calculation of new 
!              nonconstant diffusivity.
!              Based on "diffint".
!
!     diffsnapts = (false,true) =(don`t, do) save fields related 
!              to nonconstant diffusivity calculation.
!              Based on "diffsnapint".
!
!     stabts  = (false,true) = (don`t, do) test for stability on this
!               time step. based on "stabint"
!
!     tavgts  = (false,true) = (don`t do) do tracer averages on this
!               time step. based on "tavgint"
!
!     tmbts   = (false,true) = (don`t, do) write out tracer meridional .
!               balance on this time step. based on "tmbint"
!
!     tsits   = (false,true) = (don`t, do) print time step integrals
!               on this time step. based on "tsiint"
!
!     zmbcts  = (false,true) = (don`t, do) print zonal mean boundary
!               conditions on this time step.  based on "zmbcint"
!
!     trajts  = (false,true) = (don`t, do) print particle trajectories
!               on this time step. based on "trajint"
!
!     snapts  = (false,true) = (don`t, do) save a snapshot of the data
!               on this time step. based on "snapint"
!
!     timats  = (false,true) = (don`t, do) write time mean data
!               on this time step. based on "timavgint"
!                                          
!     xbtts   = (false,true) = (don`t, do) write averaged XBT data on
!               this time step based on "xbtint"
!
!
!     leapfrog= (false,true) on a (mixing, normal leapfrog) time step
!                based on "nmix"
!                                          
!     euler1  = true on the 1st pass of an euler backward time step
!               otherwise false. (applies when "eb" = true)
!     euler2  = true on the 2nd pass of an euler backward time step
!               otherwise false. (applies when "eb" = true)
!     forward = true on a forward time step. otherwise false
!		(applies when "eb" = false)		  
!                                          
!
!     the following logical switches are based on the model time step
!     and not on a specified interval in days. 
!    
!     first   = (true,false) =  when it`s (the first, not the first)
!                               time step of a run
!     eots    = end of a time step. always true except for first
!               pass of an euler backward time step
!     eorun   = last time step of a run. always false except during the
!               last time step of the run.
!
!     eoday   = true when within 1/2 time step of the end of a day
!               else ... false
!     midmon  = true when within 1/2 time step of the middle of a month
!               else ... false
!     eomon   = true when within 1/2 time step of the end of a month
!               else ... false
!     eoyear  = true when within 1/2 time step of the end of a year
!               else ... false
!     osegs   = true on the 1st time step of an ocean segment in mom.F
!               otherwise false.
!     osege  =  true on the last time step of an ocean segment in mom.F
!               otherwise false.
!
!     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
!
!
      logical eb, leapfrog, euler1, euler2, forward, eots
      logical init, first, restrt
      logical itavg, itmb, itrmb, igyre
      character *8 rununits
!
      common /switcc/ rununits
      common /switcr/ runlen, rundays
      common /switci/ nmix
      common /switcl/ eb, leapfrog, euler1, euler2, forward, eots
      common /switcl/ init, first, restrt
      common /switcl/ itavg, itmb, itrmb, igyre
!
      logical          eoday, eomon, midmon, eoyear, eorun
      common /switcl/  eoday, eomon, midmon, eoyear, eorun
!
c====================== include file "casimsw.h" ========================
c
c     holds switching information for assimilation
c     defaults and namelist settings are handled in setassim.F
c     it is "included" in switch.h
c
c     rassim  = true when within 1/2 time step of a specified interval
c               else ... false
c     rassint = interval in days for triggering a new assimilation
c
      logical rassim
      common /switcr/ rassint
      common /switcl/ rassim
c
!
!-----------------------------------------------------------------------
!
!     S W I T C H E S    B A S E D    O N    A N    I N T E R V A L
!
!     each interval switch needs two variables in common. The 
!     following naming convention is used.
!
!         1) an interval (real) for diagnosti! output (e.g,.  glenint)
!         2) a switch (logical) for the interval (e.g.,  glents )
!
!-----------------------------------------------------------------------
!
      logical           tavgts,    glents,    trmbts,   prxzts
      common /switcr/   tavgint,   glenint,   trmbint,  prxzint
      common /switcl/   tavgts,    glents,    trmbts,   prxzts
!
      logical          extts,   trajts,   exconvts,   cmixts
      common /switcr/  extint,  trajint,  exconvint,  cmixint
      common /switcl/  extts,   trajts,   exconvts,   cmixts
!
      logical          vmsfts,   gyrets,    tsits,   snapts,   fctts
      common /switcr/  vmsfint,  gyreint,   tsiint,  snapint,  fctint
      common /switcl/  vmsfts,   gyrets,    tsits,   snapts,   fctts
!
      logical          tyzts
      common /switcr/  tyzint
      common /switcl/  tyzts
!
      logical          stabts,   zmbcts,   crossts,   pressts
      common /switcr/  stabint,  zmbcint,  crossint,  pressint
      common /switcl/  stabts,   zmbcts,   crossts,   pressts
!
      logical          osegs,  osege
      common /switcr/          segtim
      common /switcl/  osegs,  osege
!
      logical          rhots
      common /switcr/  rhoint
      common /switcl/  rhots
!
      logical          diffts,  diffsnapts
      common /switcr/  diffint, diffsnapint
      common /switcl/  diffts,  diffsnapts
!   
!-----------------------------------------------------------------------
!
!     S W I T C H E S    B A S E D    O N    A N    I N T E R V A L
!
!              A N D   A V E R A G I N G   P E R I O D
!
!     each averaging period switch needs four variables in common. The 
!     following naming convention is used.
!
!         1) an interval (real) for diagnostic output    (e.g. xbtint  )
!         2) a switch (logical) for the interval         (e.g. xbtts   )
!         3) an averaging period (real)                  (e.g. xbtper  )
!         4) a switch (logical) for accumulating         (e.g. xbtperts)
!
!
!     The user must specify the interval [e.g., xbtint] for diagnostic
!     output in units of days and the averaging period [e.g., xbtper] 
!     in units of days. The averaging period may be less than or equal
!     to the interval. For example, if the interval is 30.0 days and the
!     averaging period is 5.0 days, results will be averaged over all
!     time steps within days 26, 27, 28, 29, and 30. An averaging period
!     of 0.0 days averages over the last time step of the interval (as
!     does xbtper = dt), and an averaging period less than zero turns
!     the switches off for all time steps.  
!
!     The logical switch for writing output at the specified interval 
!     [e.g., xbtts] is set to true on the last time step of the 
!     averaging period. The logical switch for accumulating results
!     [e.g., xbtperts] is true for all time steps within the averaging
!     period, otherwise it is false. 
!
!-----------------------------------------------------------------------
!
      logical         xbtts,     dspts,    tmbts,    timavgts
      logical         xbtperts,  dspperts, tmbperts, timavgperts
      common /switcr/ xbtint,    dspint,   tmbint,   timavgint
      common /switcr/ xbtper,    dspper,   tmbper,   timavgper
      common /switcl/ xbtts,     dspts,    tmbts,    timavgts
      common /switcl/ xbtperts,  dspperts, tmbperts, timavgperts
!====================== include file "taskrows.h" ======================
!
!     max_tasks      = maximum number of processors
!     num_processors = requested number of processors
!     pn             = this processor number (from 1 to num_processors)
!
!     jstask3(n)  = southernmost buffer jrow for processor "n"
!     jetask3(n)  = northernmost buffer jrow for processor "n"
!     jscomp3(n)  = southernmost computed jrow for processor "n"
!     jecomp3(n)  = northernmost computed jrow for processor "n"
!     num_loads3(n)= number of times processor "n`s" MW is loaded
!     window_size3(n)= the size of the memory window for processor "n"
!
!     jstask  = southernmost buffer jrow for this processor (pn)
!     jetask  = northernmost buffer jrow for this processor (pn)
!     jscomp  = southernmost computed jrow for this processor (pn)
!     jecomp  = northernmost computed jrow for this processor (pn)
!     num_loads= number of times processor "pn`s" MW is loaded
!     window_size= the size of the memory window for this processor (pn)
!                  It is set as follows:
!                  if defined max_window then
!                    jecomp3-jscomp3+1 + 2*jbuf
!                    for second, fourth, and sixth order windows
!                    on processor "pn". The window_size must be the
!                    same on all processors. 
!                  else
!                    3 (for second order windows)
!                    5 (for fourth order windows)
!                    7 (for sixth order windows)
!                    regardless of the number of processors
!                  endif
!
      parameter (max_tasks = 2048)
      integer pn, window_size3, window_size
      common /taskrows/ pn
      common /taskrows/ num_processors
      common /taskrows/ jstask3(max_tasks), jetask3(max_tasks)
      common /taskrows/ jscomp3(max_tasks), jecomp3(max_tasks)
      common /taskrows/ num_loads3(max_tasks)
      common /taskrows/ jstask, jetask, num_loads, jscomp, jecomp
      common /taskrows/ window_size3(max_tasks), window_size
!====================== include file "tmngr.h" =========================
!
!                       time manager variables
!
!     time manager inputs:
!
!     refrun  = T ==> logical switches referenced to the start of run
!     refinit = T ==> logical switches referenced to initial condition
!                     time: year0, month0, day0, hour0, min0, sec0
!     refuser = T ==> logical switches referenced to user specified 
!                     time: ryear, rmonth, rday, rhour, rmin, rsec
!
!     dayoyr  = relative day number referenced to the beginning
!               of the current year.  (real)
!     relyr   = number of years (and fractional years) of model
!               integration (for time tau+1 {itt}) relative to 
!               initial condition
!     prelyr  = relyr for previous time step
!
!     stamp   = 32 character date and time for current model timestep
!     pstamp  = 32 character date and time for previous model timestep
!     ic_stamp= 32 character date and time for initial conditions
!     stop_stamp = 32 character date and time for last time step
!
!     itt     = current time step counter (from initial cond.)
!     itt0    = time step at start of current run
!
!               variables used for initialization
!
!     irstdy  = integer number of days at start of run
!     msrsdy  = fractional day in millisec at start of run
!
!     year0   = year of initial conditions
!     month0  = month of initial conditions
!     day0    = day of initial conditions
!     hour0   = hour of initial conditions
!     min0    = minute of initial conditions
!     sec0    = second of initial conditions  
!     
!     ryear   = year of user specified reference time
!     rmonth  = month of user specified reference time
!     rday    = day of user specified reference time
!     rhour   = hour of user specified reference time
!     rmin    = minute of user specified reference time
!     rsec    = second of user specified reference time  
!
!-----------------------------------------------------------------------
!
!
      logical refrun, refinit, refuser, eqmon, eqyear
!
      integer year0, month0, day0, hour0, min0, sec0
      integer ryear, rmonth, rday, rhour, rmin, rsec
      character *32 pstamp, stamp, ic_stamp, stop_stamp, xstamp
!
      common /tmngrl/ refrun, refinit, refuser, eqmon, eqyear
!
      common /tmngrii/ itt0, itt, irstdy, msrsdy
!
      common /tmngrr/ dayoyr, relyr, prelyr, relday
!
      common /tmngrc/ stamp, pstamp, ic_stamp, stop_stamp, xstamp
!
      common /tmngrii/ year0, month0, day0, hour0, min0, sec0
      common /tmngrii/ ryear, rmonth, rday, rhour, rmin, rsec
!
      character*10 dayname
      character*12 monname
      common /tmngrc/ dayname(7), monname(12)
!
      real yrlen
      common /tmngrr/ yrlen
!
      character(len=80) :: timeunits
      common /tmngrc/ timeunits
      
      save cps, itt1st
c
      logical goto_2nd_pass
c     
c-----------------------------------------------------------------------
c     integrate equations for one time step
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c     update timestep counter, set corresponding model time, and set
c     time dependent logical switches which determine program flow.
c-----------------------------------------------------------------------
c
      itt = itt + 1
      call tmngr (dtts)
      call set_time_switches
c
      write(stdout,'(a,i7,1x,a32)') ' ts=',itt,stamp
      call flush(stdout)
c
c-----------------------------------------------------------------------
c     read system clock to calculate compute_time/grid_point/timestep 
c     (neglecting overhead of "setocn")
c-----------------------------------------------------------------------
c
      if (first) then
	itt1st = itt
        cps    = cpusec()
      endif
c
c-----------------------------------------------------------------------
c     initialize diagnostic variables
c-----------------------------------------------------------------------
c
      call diag_initial (jscomp, jecomp)
c
c-----------------------------------------------------------------------
c     adjust various quantities for leapfrog/mixing timesteps
c
c     leapfrog----------> h(tau+1) = h(tau-1) + 2dt*F(tau)
c                          
c     forward-----------> tau-1 <= tau
c                           h(tau+1) = h(tau-1) + dt*F(tau)
c
c     euler backward:     tau-1 <= tau
c       euler1----------> h(tau` ) = h(tau-1) + dt*F(tau)
c       euler2----------> h(tau+1) = h(tau-1) + dt*F(tau`)
c-----------------------------------------------------------------------
c
      if (leapfrog) then
c
c       normal leapfrog time step
c
        euler1  = .false.
	euler2  = .false.
	forward = .false.
	eots    = .true.
c
        c2dtts  = c2*dtts
        c2dtuv  = c2*dtuv
        c2dtsf  = c2*dtsf
      else
c
c       mixing time step (forward step or euler backward step)
c
        if (eb) then
          euler1  = .true.
	  euler2  = .false.
	  forward = .false.
	  eots    = .false.
	else
          euler1  = .false.
	  euler2  = .false.
	  forward = .true.
	  eots    = .true.
	endif
c
        c2dtts = dtts
        c2dtuv = dtuv
        c2dtsf = dtsf
c
      endif
c
        if (.not.leapfrog) then
          do jrow=jstask,jetask
            do i=1,imt
              ubarm1(i,jrow,1) = ubar(i,jrow,1)
              ubarm1(i,jrow,2) = ubar(i,jrow,2)
            enddo
          enddo
        endif
c
c
c-----------------------------------------------------------------------
c     set time centering "gcor" for coriolis term
c-----------------------------------------------------------------------
c
      if (acor .eq. c0) then
        gcor = c1
      elseif (acor .ne. c0) then
	gcor = c0
      endif
c
c
c-----------------------------------------------------------------------
c     update pointers to tau-1, tau, & tau+1 data on disk.
c     for latitude rows they point to latdisk
c     for 2D fields they point to records on kflds
c-----------------------------------------------------------------------
c
      taum1disk = mod(itt+1,2) + 1
      taudisk   = mod(itt  ,2) + 1
      taup1disk = taum1disk
c
c-----------------------------------------------------------------------
c     update pointers (indices) to tau-1, tau, & tau+1 data in the MW
c-----------------------------------------------------------------------
c
c
c       time levels are being held constant in time.
c
c
c=======================================================================
c
c               SOLVE THE BAROCLINIC AND TRACER EQUATIONS
c
c     Latitude rows are solved using a MW (memory window) approach.
c     Refer to the manual for  details. 
c     
c=======================================================================
c
      goto_2nd_pass = .true.
      do while (goto_2nd_pass)
c
c
c-----------------------------------------------------------------------
c     Load the processor`s memory window in preparation for solving
c     tracer and baroclinic equations. (See the manual)
c-----------------------------------------------------------------------
c
      do mw=1,num_loads
c
        ncrows = window_size-2*jbuf
	joff = jscomp - jbuf - 1 + (mw-1)*ncrows
        js = (jstask - joff)*(1/mw) + (1 + num_rows_to_copy)*(1-1/mw)
	je = min(window_size,jetask-joff)
c
	is = 2
	ie = imt - 1
c
c       set rows  within the MW for solving prognostic equations
c       (rows are clipped at jrow=2 and jrow=jmt-1)
c
        js_calc = 1+jbuf
        je_calc = ncrows+jbuf
c
c-----------------------------------------------------------------------
c       load prognostic and related variables into the MW
c-----------------------------------------------------------------------
c
        call loadmw (joff, js, je, is, ie, latdisk(taum1disk)
     &,                latdisk(taudisk))
c
c-----------------------------------------------------------------------
c       add correction from assimilation
c-----------------------------------------------------------------------
c
!       call add_correction (joff, js_calc, je_calc)
        call add_correction (joff, js, je)
c-----------------------------------------------------------------------
c       compute pressure gradients
c-----------------------------------------------------------------------
c
        call pressure_gradients (joff, js, je)
c
c-----------------------------------------------------------------------
c       calculate advection velocities for T cells and U cells     
c-----------------------------------------------------------------------
c
        call adv_vel (joff, js, je, is, ie)
c
c-----------------------------------------------------------------------
c       set horizontal mixing coefficients for momentum and tracers
c-----------------------------------------------------------------------
c
        call hmixc (joff, js, je, is, ie)
c
c
c-----------------------------------------------------------------------
c       calculate isopycnal diffusion tensor components (and
c       1 advective velocities) for use with tracers
c-----------------------------------------------------------------------
c
        call isopyc (joff, js, je, is, ie)
c
c-----------------------------------------------------------------------
c       set vertical boundary conditions for momentum and tracers
c-----------------------------------------------------------------------
c
	call setvbc (joff, js, je, is, ie)
c
c-----------------------------------------------------------------------
c       set vertical mixing coefficients for momentum and tracers
c-----------------------------------------------------------------------
c
        call vmixc (joff, js, je, is, ie)
c
c
c
c
c
c-----------------------------------------------------------------------
c       compute tracers and internal mode velocities
c-----------------------------------------------------------------------
c
        call tracer     (joff, js_calc, je_calc, is, ie)
        call baroclinic (joff, js_calc, je_calc, is, ie)
c
c-----------------------------------------------------------------------
c       extract first guess for assimilation
c-----------------------------------------------------------------------
c
        if (nassim .ge. 1 .and. nassim .le. nassdur) then
!         call extract_guess (joff, js_calc, je_calc)
          call extract_guess (joff, js, je)
        endif
c
c-----------------------------------------------------------------------
c       calculate diagnostics
c-----------------------------------------------------------------------
c
        call diag_mw_inside (joff, js_calc, je_calc, is, ie)
c
c
c-----------------------------------------------------------------------
c       robert time filter for velocity and tracers
c-----------------------------------------------------------------------
c
        call r_time_filter (joff, js_calc, je_calc, is, ie
     &,                     latdisk(taudisk))

c
c-----------------------------------------------------------------------
c       write prognostic variables from the MW to disk "tau+1"
c-----------------------------------------------------------------------
c
          call putmw (joff, js_calc, je_calc, latdisk(taup1disk))
c
      enddo
c
c-----------------------------------------------------------------------
c     update boundary latitudes when using more than one processor
c-----------------------------------------------------------------------
c
!     update tau and tau+1 disk
      call mpp_update_domains( ramdisk3d, mom_ramdisk_domain )
c
c=======================================================================
c
c               SOLVE THE BAROTROPIC EQUATION
c
c=======================================================================
c
      call barotropic (cori, itt, jstask, jetask)
c
c-----------------------------------------------------------------------
c     if this is the 1st pass of an euler backward timestep, set the
c     disk pointers so the proper time levels are read on the 2nd pass 
c     and go back to do the 2nd pass.
c-----------------------------------------------------------------------
c
      if (euler1) then
        eots      = .true.
	euler1    = .false.
	euler2    = .true.
        ntaux     = taum1disk
        taum1disk = taudisk
        taudisk   = taup1disk
        taup1disk = ntaux
      else
        goto_2nd_pass = .false.
      endif
c
      enddo
c
c
c-----------------------------------------------------------------------
c     output all remaining diagnostics
c-----------------------------------------------------------------------
c
      call diag_mw_outside(jscomp, jecomp)
c
c-----------------------------------------------------------------------
c     if it`s the last timestep then clean things up otherwise return
c-----------------------------------------------------------------------
c
      if (eorun) then
c
c       calculate cp_time/grid_point/timestep (neglecting setocn)
c
        ncalls = itt - itt1st + 1
        cpe    = cpusec()
        write (stdout,'(/,1x, a, f9.3, a,1x,a,1pg14.7,a)')
     &    'Timing Summary: ', cpe-cps,  ' cp seconds,'
     &,    '(',(cpe-cps)/(ncalls*imt*jmt*km)
     &,    ' cp seconds per grid point per time step)'
c
c       write a restart
c
        if (restrt) call write_restart (jstask, jetask)
c
        write (stdout,'(1x,a)') 'MOMdone'
        call flush(stdout)
        call oclose (kflds)
	call oclose (latdisk(1))
	call oclose (latdisk(2))
      endif
c
c
      return
      end


      function cpusec ()
c
c     user supplied function to measure cpu time in seconds
c
c     real(kind=4) :: second
      cpusec = second()
      return
      end


