[Dart-dev] DART/branches Revision: 10740

dart at ucar.edu dart at ucar.edu
Fri Nov 11 10:21:37 MST 2016


thoar at ucar.edu
2016-11-11 10:21:36 -0700 (Fri, 11 Nov 2016)
339
Moving shell scripts to the shell_scripts directory.
Removing pieces for dart_to_pop and pop_to_dart.
Removing unused 'broken_mkmf_model_mod_check' script since the mkmf_model_mod_check script exists and works.
Removing pop_model_mod.f90 since this strategy has been replaced for strongly coupled DA. That version was just a stub anyway.




Deleted: DART/branches/rma_trunk/models/POP/dart_to_pop.f90
===================================================================
--- DART/branches/rma_trunk/models/POP/dart_to_pop.f90	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/dart_to_pop.f90	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,136 +0,0 @@
-! DART software - Copyright 2004 - 2013 UCAR. This open source software is
-! provided by UCAR, "as is", without charge, subject to all terms of use at
-! http://www.image.ucar.edu/DAReS/DART/DART_download
-!
-! $Id$
-
-program dart_to_pop
-
-!----------------------------------------------------------------------
-! purpose: interface between DART and the POP model
-!
-! method: Read DART state vector and overwrite values in a POP restart file.
-!         If the DART state vector has an 'advance_to_time' present, a
-!         file called pop_in.DART is created with a time_manager_nml namelist 
-!         appropriate to advance POP to the requested time.
-!
-!         The dart_to_pop_nml namelist setting for advance_time_present 
-!         determines whether or not the input file has an 'advance_to_time'.
-!         Typically, only temporary files like 'assim_model_state_ic' have
-!         an 'advance_to_time'.
-!
-! author: Tim Hoar 25 Jun 09, revised 12 July 2010
-!----------------------------------------------------------------------
-
-use        types_mod, only : r8
-use    utilities_mod, only : initialize_utilities, finalize_utilities, &
-                             find_namelist_in_file, check_namelist_read, &
-                             logfileunit
-use  state_vector_io_mod, only : open_restart_read, aread_state_restart, close_restart
-use time_manager_mod, only : time_type, print_time, print_date, operator(-)
-use        model_mod, only : static_init_model, sv_to_restart_file, &
-                             get_model_size, get_pop_restart_filename
-use     dart_pop_mod, only : write_pop_namelist
-
-implicit none
-
-! version controlled file description for error handling, do not edit
-character(len=256), parameter :: source   = &
-   "$URL$"
-character(len=32 ), parameter :: revision = "$Revision$"
-character(len=128), parameter :: revdate  = "$Date$"
-
-!------------------------------------------------------------------
-! The namelist variables
-!------------------------------------------------------------------
-
-character (len = 128) :: dart_to_pop_input_file   = 'dart_restart'
-logical               :: advance_time_present     = .false.
-
-namelist /dart_to_pop_nml/ dart_to_pop_input_file, &
-                           advance_time_present
-
-!----------------------------------------------------------------------
-
-integer               :: iunit, io, x_size
-type(time_type)       :: model_time, adv_to_time
-real(r8), allocatable :: statevector(:)
-character (len = 128) :: pop_restart_filename = 'no_pop_restart_file'
-
-!----------------------------------------------------------------------
-
-call initialize_utilities(progname='dart_to_pop')
-
-!----------------------------------------------------------------------
-! Call model_mod:static_init_model() which reads the POP namelists
-! to set grid sizes, etc.
-!----------------------------------------------------------------------
-
-call static_init_model()
-
-x_size = get_model_size()
-allocate(statevector(x_size))
-
-! Read the namelist to get the input filename. 
-
-call find_namelist_in_file("input.nml", "dart_to_pop_nml", iunit)
-read(iunit, nml = dart_to_pop_nml, iostat = io)
-call check_namelist_read(iunit, io, "dart_to_pop_nml")
-
-call get_pop_restart_filename( pop_restart_filename )
-
-write(*,*)
-write(*,'(''dart_to_pop:converting DART file '',A, &
-      &'' to POP restart file '',A)') &
-     trim(dart_to_pop_input_file), trim(pop_restart_filename)
-
-!----------------------------------------------------------------------
-! Reads the valid time, the state, and the target time.
-!----------------------------------------------------------------------
-
-iunit = open_restart_read(dart_to_pop_input_file)
-
-if ( advance_time_present ) then
-   call aread_state_restart(model_time, statevector, iunit, adv_to_time)
-else
-   call aread_state_restart(model_time, statevector, iunit)
-endif
-call close_restart(iunit)
-
-!----------------------------------------------------------------------
-! update the current POP state vector
-! Convey the amount of time to integrate the model ...
-! time_manager_nml: stop_option, stop_count increments
-!----------------------------------------------------------------------
-
-call sv_to_restart_file(statevector, pop_restart_filename, model_time)
-
-if ( advance_time_present ) then
-   call write_pop_namelist(model_time, adv_to_time)
-endif
-
-!----------------------------------------------------------------------
-! Log what we think we're doing, and exit.
-!----------------------------------------------------------------------
-
-call print_date( model_time,'dart_to_pop:POP  model date')
-call print_time( model_time,'dart_to_pop:DART model time')
-call print_date( model_time,'dart_to_pop:POP  model date',logfileunit)
-call print_time( model_time,'dart_to_pop:DART model time',logfileunit)
-
-if ( advance_time_present ) then
-call print_time(adv_to_time,'dart_to_pop:advance_to time')
-call print_date(adv_to_time,'dart_to_pop:advance_to date')
-call print_time(adv_to_time,'dart_to_pop:advance_to time',logfileunit)
-call print_date(adv_to_time,'dart_to_pop:advance_to date',logfileunit)
-endif
-
-call finalize_utilities('dart_to_pop')
-
-end program dart_to_pop
-
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$

Deleted: DART/branches/rma_trunk/models/POP/dart_to_pop.html
===================================================================
--- DART/branches/rma_trunk/models/POP/dart_to_pop.html	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/dart_to_pop.html	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,270 +0,0 @@
-<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
-"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
-<HTML>
-<HEAD>
-<TITLE>program dart_to_pop</TITLE>
-<link rel="stylesheet" type="text/css" href="../../doc/html/doc.css" />
-<link href="../../doc/images/dart.ico" rel="shortcut icon" />
-</HEAD>
-<BODY>
-<A NAME="TOP"></A>
-
-<H1>PROGRAM <em class=program>dart_to_pop</em></H1>
-
-<table border=0 summary="" cellpadding=5>
-<tr>
-    <td valign=middle>
-    <img src="../../doc/images/Dartboard7.png" alt="DART project logo" height=70 />
-    </td>
-    <td>
-       <P>Jump to <a href="../../index.html">DART Documentation Main Index</a><br />
-          <small><small>version information for this file: <br />
-          <!-- version tag follows, do not edit -->
-          $Id$</small></small>
-       </P></td>
-</tr>
-</table>
-
-<A HREF="#Namelist">NAMELIST</A> /
-<A HREF="#Modules">MODULES</A> /
-<A HREF="#FilesUsed">FILES</A> /
-<A HREF="#References">REFERENCES</A> /
-<A HREF="#Errors">ERRORS</A> /
-<A HREF="#FuturePlans">PLANS</A> /
-<A HREF="#Legalese">TERMS OF USE</A>
-
-<H2>Overview</H2>
-
-<P>
-   <em class=program>dart_to_pop</em> is the program that updates a POP
-   netCDF-format restart file (usually <em class=file>pop.r.nc</em>) 
-   with the state information contained in a DART output/restart file 
-   (e.g. <em class=file>perfect_ics, filter_ics, ... </em>).
-   Only the CURRENT values in the POP restart file will be updated: 
-   specifically -
-   <em class=code>SALT_CUR</em>,
-   <em class=code>TEMP_CUR</em>,
-   <em class=code>UVEL_CUR</em>,
-   <em class=code>VVEL_CUR</em>, and
-   <em class=code>PSURF_CUR</em>. It is <strong>necessary</strong> to 
-   perform a forward euler timestep since the '<em class=code>_OLD</em>'
-   variables are untouched by DART.  The DART model time is compared to the
-   time in the POP restart file. If they are not identical, the program
-   issues an error message and aborts.
-   <br />
-   <br />
-   From the user perspective, most of the time 
-   <em class=program>dart_to_pop</em> will be used on DART files that
-   have a header containing one time stamp followed by the model state.
-   <br />
-   <br />
-   The <a href=#Namelist>dart_to_pop_nml</a> namelist allows 
-   <em class=program>dart_to_pop</em> to read the 
-   <em class=file>assim_model_state_ic</em> files that have
-   <em class=italic>two</em> timestamps in the header. These files are 
-   temporarily generated when DART is used to advance the model. 
-   One timestamp is the 'advance_to' time, the other is the 'valid_time' 
-   of the model state. In this case, a namelist for POP (called 
-   <em class=file>pop_in.DART</em>) is written that contains the 
-   <em class=code>&amp;time_manager_nml</em> settings appropriate to 
-   advance POP to the time requested by DART. The repository version 
-   of the <em class=program>advance_model.csh</em> script has a section 
-   to ensure the proper DART namelist settings for this case.
-   <br />
-   <br />
-   Conditions required for successful execution of <em class=program>dart_to_pop</em>:
-</P>
-
-<UL>
-   <LI>a valid <em class=file>input.nml</em> namelist file for DART</LI>
-   <LI>a DART file (typically <em class=file>filter_restart.xxxx</em> or
-                              <em class=file>filter_ics.xxxx</em>)</LI>
-   <LI>a valid <em class=file>pop_in</em> namelist file for POP</LI>
-   <LI>the POP geometry files mentioned in the <em class=file>pop_in</em> namelist file</LI>
-   <LI>a POP restart file (typically <em class=file>pop.r.nc</em>). 
-   Keep in mind that the DART/POP interface requires the existence of 
-   a <em class=file>pop.r.nc</em> to harvest the geometry dimensions.</LI>
-</UL>
-
-<P>
-Since this program is called repeatedly for every ensemble member,
-we have found it convenient to link the DART input file
-to the default input filename (<em class=file>dart_restart</em>). The same
-thing goes true for the POP output filename <em class=file>pop.r.nc</em>.
-</P>
-
-<!--==================================================================-->
-<!--=================== DESCRIPTION OF A NAMELIST ====================-->
-<!--==================================================================-->
-
-<A NAME="Namelist"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>NAMELIST</H2>
-<P>
-This namelist is read from the file <em class=file>input.nml</em>.
-Namelists start with an ampersand
-'&amp;' and terminate with a slash '/'.
-Character strings that contain a '/' must be
-enclosed in quotes to prevent them from 
-prematurely terminating the namelist.
-</P>
-
-<div class=namelist>
-<pre>
-&amp;dart_to_pop_nml
-   dart_to_pop_input_file = 'dart_restart',
-   advance_time_present   = .false.  
-/
-</pre>
-</div>
-
-<br />
-<br />
-
-<div>
-<TABLE border=0 cellpadding=10 width=100% summary='namelist description'>
-<THEAD align=left>
-<TR><TH> Item </TH>
-    <TH> Type </TH>
-    <TH> Description </TH> </TR>
-</THEAD>
-
-<TBODY valign=top>
-
-<TR><TD>dart_to_pop_input_file   </TD>
-    <TD>character(len=128)   </TD>
-    <TD>The name of the DART file containing the model state
-to insert into the POP restart file.
-</TD></TR>
-
-<TR><TD>advance_time_present</TD>
-    <TD>logical</TD>
-    <TD>If you are converting a DART initial conditions or restart file this
-should be <em class=code>.false.</em>; these files have a single timestamp
-describing the valid time of the model state.  If <em class=code>.true.</em>,
-TWO timestamps are expected to be the DART file header.  In this case, a
-namelist for POP (called <em class=file>pop_in.DART</em>) is created that
-contains the <em class=code>&amp;time_manager_nml</em> settings appropriate
-to advance POP to the time requested by DART.
-</TD></TR>
-
-</TBODY> 
-</TABLE>
-</div>
-
-<br />
-<br />
-
-
-<!--==================================================================-->
-
-<A NAME="Modules"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>MODULES USED</H2>
-<PRE>
-assim_model_mod
-dart_pop_mod
-location_mod
-model_mod
-null_mpi_utilities_mod
-obs_kind_mod
-random_nr_mod
-random_seq_mod
-time_manager_mod
-types_mod
-utilities_mod
-</PRE>
-
-<!--==================================================================-->
-<!-- Describe the Files Used by this module.                          -->
-<!--==================================================================-->
-
-<A NAME="FilesUsed"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>FILES Read</H2>
-<UL><LI>DART initial conditions/restart file; e.g. <em class=file>filter_ic</em></LI>
-    <LI>DART namelist file; <em class=file>input.nml</em></LI>
-    <LI>POP namelist file; <em class=file>pop_in</em></LI>
-    <LI>POP geometry files specified in <em class=file>pop_in</em></LI>
-    <LI>POP restart file <em class=file>pop.r.nc</em> (to get grid values)</LI>
-</UL>
-
-<H2>FILES Written</H2>
-<UL><LI>POP restart file; <em class=file>pop.r.nc</em></LI>
-    <LI>POP namelist file; <em class=file>pop_in.DART</em></LI>
-</UL>
-
-
-<!--==================================================================-->
-<!-- Cite references, if need be.                                     -->
-<!--==================================================================-->
-
-<A NAME="References"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>REFERENCES</H2>
-<ul>
-<li>Anderson,&nbsp;J.,&nbsp;T.&nbsp;Hoar,&nbsp;K.&nbsp;Raeder,
-    H.&nbsp;Liu,&nbsp;N.&nbsp;Collins,&nbsp;R.&nbsp;Torn,
-    and&nbsp;A.&nbsp;Arellano,&nbsp;2009:<br />
-    The Data Assimilation Research Testbed: A Community Facility. 
-    <span style="font-style: italic;">Bull. Amer. Meteor. Soc.</span>,
-    <span style="font-weight: bold;">90</span>, 1283-1296.<br />
-    <a href="http://ams.allenpress.com/perlserv/?doi=10.1175%2F2009BAMS2618.1&amp;request=get-abstract">DOI: 10.1175/2009BAMS2618.1</a></li>
-</ul>
-<br />
-
-<!--==================================================================-->
-<!-- Describe all the error conditions and codes.                     -->
-<!--==================================================================-->
-
-<A NAME="Errors"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>ERROR CODES and CONDITIONS</H2>
-<P>
-none - all error messages come from modules that have their own documentation.
-</P>
-
-<H2>KNOWN BUGS</H2>
-<P>
-none
-</P>
-
-<!--==================================================================-->
-<!-- Describe Future Plans.                                           -->
-<!--==================================================================-->
-
-<A NAME="FuturePlans"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>FUTURE PLANS</H2>
-<P>
-None.
-</P>
-
-<!--==================================================================-->
-<!-- Legalese & Metadata                                              -->
-<!--==================================================================-->
-
-<A NAME="Legalese"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>Terms of Use</H2>
-
-<P>
-DART software - Copyright 2004 - 2013 UCAR.<br />
-This open source software is provided by UCAR, "as is",<br />
-without charge, subject to all terms of use at<br />
-<a href="http://www.image.ucar.edu/DAReS/DART/DART_download">
-http://www.image.ucar.edu/DAReS/DART/DART_download</a>
-</P>
-
-<TABLE border=0 cellpadding=0 width=100% summary="">
-<TR><TD valign=top>Contact:       </TD><TD> DART core group   </TD></TR>
-<TR><TD valign=top>Revision:      </TD><TD> $Revision$ </TD></TR>
-<TR><TD valign=top>Source:        </TD><TD> $URL$ </TD></TR>
-<TR><TD valign=top>Change Date:   </TD><TD> $Date$ </TD></TR>
-<TR><TD valign=top>Change&nbsp;history:&nbsp;</TD><TD> try "svn&nbsp;log" or "svn&nbsp;diff" </TD></TR>
-</TABLE>
-
-<!--==================================================================-->
-
-</BODY>
-</HTML>

Deleted: DART/branches/rma_trunk/models/POP/dart_to_pop.nml
===================================================================
--- DART/branches/rma_trunk/models/POP/dart_to_pop.nml	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/dart_to_pop.nml	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,4 +0,0 @@
-&dart_to_pop_nml
-   dart_to_pop_input_file = 'dart_restart',
-   advance_time_present   = .false.  /
-

Deleted: DART/branches/rma_trunk/models/POP/pop_model_mod.f90
===================================================================
--- DART/branches/rma_trunk/models/POP/pop_model_mod.f90	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/pop_model_mod.f90	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,3630 +0,0 @@
-! DART software - Copyright 2004 - 2013 UCAR. This open source software is
-! provided by UCAR, "as is", without charge, subject to all terms of use at
-! http://www.image.ucar.edu/DAReS/DART/DART_download
-!
-! $Id$
-
-module pop_model_mod
-
-! This is the interface between the POP ocean model and DART.
-
-! Modules that are absolutely required for use are listed
-use        types_mod, only : r4, r8, SECPERDAY, MISSING_R8, rad2deg, PI
-use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
-                             print_time, print_date,                           &
-                             operator(*),  operator(+), operator(-),           &
-                             operator(>),  operator(<), operator(/),           &
-                             operator(/=), operator(<=)
-use     location_mod, only : location_type, get_dist,                          &
-                             get_close_maxdist_init, get_close_obs_init,       &
-                             set_location,                                     &
-                             VERTISHEIGHT, get_location, vert_is_height,       &
-                             vert_is_level, vert_is_surface,                   &
-                             loc_get_close_obs => get_close_obs, get_close_type
-use    utilities_mod, only : register_module, error_handler,                   &
-                             E_ERR, E_WARN, E_MSG, logfileunit, get_unit,      &
-                             nc_check, do_output, to_upper,                    &
-                             find_namelist_in_file, check_namelist_read,       &
-                             open_file, file_exist, find_textfile_dims,        &
-                             file_to_text, do_output
-use     obs_kind_mod, only : KIND_TEMPERATURE, KIND_SALINITY, KIND_DRY_LAND,   &
-                             KIND_U_CURRENT_COMPONENT,                         &
-                             KIND_V_CURRENT_COMPONENT, KIND_SEA_SURFACE_HEIGHT, &
-                             KIND_SEA_SURFACE_PRESSURE, KIND_POTENTIAL_TEMPERATURE
-use mpi_utilities_mod, only: my_task_id
-use    random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
-use      dart_pop_mod, only: set_model_time_step,                              &
-                             get_horiz_grid_dims, get_vert_grid_dim,           &
-                             read_horiz_grid, read_topography, read_vert_grid, &
-                             get_pop_restart_filename
-
-use typesizes
-use netcdf 
-
-implicit none
-private
-
-! these routines must be public and you cannot change
-! the arguments - they will be called *from* the DART code.
-public :: pop_get_model_size,         &
-          pop_adv_1step,              &
-          pop_get_state_meta_data,    &
-          pop_model_interpolate,      &
-          pop_get_model_time_step,    &
-          pop_static_init_model,      &
-          pop_end_model,              &
-          pop_init_time,              &
-          pop_init_conditions,        &
-          pop_nc_write_model_atts,    &
-          pop_nc_write_model_vars,    &
-          pop_pert_model_state,       &
-          get_close_maxdist_init,     &
-          get_close_obs_init,         &
-          pop_get_close_obs,          &
-          pop_ens_mean_for_model
-
-! generally useful routines for various support purposes.
-! the interfaces here can be changed as appropriate.
-public :: pop_get_gridsize, pop_restart_file_to_sv, pop_sv_to_restart_file, &
-          get_pop_restart_filename, pop_test_interpolation
-
-! version controlled file description for error handling, do not edit
-character(len=256), parameter :: source   = &
-   "$URL$"
-character(len=32 ), parameter :: revision = "$Revision$"
-character(len=128), parameter :: revdate  = "$Date$"
-
-character(len=256) :: msgstring
-logical, save :: module_initialized = .false.
-
-! Storage for a random sequence for perturbing a single initial state
-type(random_seq_type) :: random_seq
-
-! things which can/should be in the model_nml
-logical  :: output_state_vector = .true.
-integer  :: assimilation_period_days = 1
-integer  :: assimilation_period_seconds = 0
-real(r8) :: model_perturbation_amplitude = 0.2
-logical  :: update_dry_cell_walls = .false.
-integer  :: debug = 0   ! turn up for more and more debug messages
-
-! FIXME: currently the update_dry_cell_walls namelist value DOES
-! NOTHING.  it needs additional code to detect the cells which are
-! wet, but within 1 cell of the bottom/sides/etc.  
-
-namelist /pop_model_nml/  &
-   output_state_vector,         &
-   assimilation_period_days,    &  ! for now, this is the timestep
-   assimilation_period_seconds, &
-   model_perturbation_amplitude,&
-   update_dry_cell_walls,       &
-   debug
-
-!------------------------------------------------------------------
-!
-! The DART state vector (control vector) will consist of:  S, T, U, V, PSURF
-! (Salinity, Temperature, U velocity, V velocity, Sea Surface Height).
-! S, T are 3D arrays, located at cell centers.  U,V are at grid cell corners.
-! PSURF is a 2D field (X,Y only).  The Z direction is downward.
-!
-! FIXME: proposed change 1: we put SSH first, then T,U,V, then S, then
-!                           any optional tracers, since SSH is the only 2D
-!                           field; all tracers are 3D.  this simplifies the
-!                           mapping to and from the vars to state vector.
-!
-! FIXME: proposed change 2: we make this completely namelist driven,
-!                           both contents and order of vars.  this should
-!                           wait until restart files are in netcdf format,
-!                           to avoid problems with incompatible namelist
-!                           and IC files.  it also complicates the mapping
-!                           to and from the vars to state vector.
-!------------------------------------------------------------------
-
-integer, parameter :: n3dfields = 4
-integer, parameter :: n2dfields = 1
-integer, parameter :: nfields   = n3dfields + n2dfields
-
-! (the absoft compiler likes them to all be the same length during declaration)
-! we trim the blanks off before use anyway, so ...
-character(len=128) :: progvarnames(nfields) = (/'SALT ','TEMP ','UVEL ','VVEL ','PSURF'/)
-
-integer, parameter :: S_index     = 1
-integer, parameter :: T_index     = 2
-integer, parameter :: U_index     = 3
-integer, parameter :: V_index     = 4
-integer, parameter :: PSURF_index = 5
-
-integer :: start_index(nfields)
-
-! Grid parameters - the values will be read from a
-! standard POP namelist and filled in here.
-
-! nx, ny and nz are the size of the dipole (or irregular) grids. 
-integer :: Nx=-1, Ny=-1, Nz=-1    ! grid counts for each field
-
-! locations of cell centers (C) and edges (G) for each axis.
-real(r8), allocatable :: ZC(:), ZG(:)
-
-! These arrays store the longitude and latitude of the lower left corner of
-! each of the dipole u quadrilaterals and t quadrilaterals.
-real(r8), allocatable :: ULAT(:,:), ULON(:,:), TLAT(:,:), TLON(:,:)
-
-! integer, lowest valid cell number in the vertical
-integer, allocatable  :: KMT(:, :), KMU(:, :)
-! real, depth of lowest valid cell (0 = land).  use only if KMT/KMU not avail.
-real(r8), allocatable :: HT(:,:), HU(:,:)
-
-! compute pressure based on depth - can do once upfront.
-real(r8), allocatable :: pressure(:)
-
-real(r8)        :: endTime
-real(r8)        :: ocean_dynamics_timestep = 900.0_r4
-integer         :: timestepcount = 0
-type(time_type) :: model_time, model_timestep
-
-integer :: model_size    ! the state vector length
-
-
-INTERFACE vector_to_prog_var
-      MODULE PROCEDURE vector_to_2d_prog_var
-      MODULE PROCEDURE vector_to_3d_prog_var
-END INTERFACE
-
-
-!------------------------------------------------
-
-! The regular grid used for dipole interpolation divides the sphere into
-! a set of regularly spaced lon-lat boxes. The number of boxes in
-! longitude and latitude are set by num_reg_x and num_reg_y. Making the
-! number of regular boxes smaller decreases the computation required for
-! doing each interpolation but increases the static storage requirements
-! and the initialization computation (which seems to be pretty small).
-integer, parameter :: num_reg_x = 90, num_reg_y = 90
-
-! The max_reg_list_num controls the size of temporary storage used for
-! initializing the regular grid. Four arrays
-! of size num_reg_x*num_reg_y*max_reg_list_num are needed. The initialization
-! fails and returns an error if max_reg_list_num is too small. A value of
-! 30 is sufficient for the x3 POP grid with 180 regular lon and lat boxes 
-! and a value of 80 is sufficient for for the x1 grid.
-integer, parameter :: max_reg_list_num = 80
-
-! The dipole interpolation keeps a list of how many and which dipole quads
-! overlap each regular lon-lat box. The number for the u and t grids are 
-! stored in u_dipole_num and t_dipole_num. The allocatable arrays
-! u_dipole_lon(lat)_list and t_dipole_lon(lat)_list list the longitude 
-! and latitude indices for the overlapping dipole quads. The entry in
-! u_dipole_start and t_dipole_start for a given regular lon-lat box indicates
-! where the list of dipole quads begins in the u_dipole_lon(lat)_list and
-! t_dipole_lon(lat)_list arrays.
-
-integer :: u_dipole_start(num_reg_x, num_reg_y)
-integer :: u_dipole_num  (num_reg_x, num_reg_y) = 0
-integer :: t_dipole_start(num_reg_x, num_reg_y)
-integer :: t_dipole_num  (num_reg_x, num_reg_y) = 0
-integer, allocatable :: u_dipole_lon_list(:), t_dipole_lon_list(:)
-integer, allocatable :: u_dipole_lat_list(:), t_dipole_lat_list(:)
-
-! Need to check for pole quads: for now we are not interpolating in them
-integer :: pole_x, t_pole_y, u_pole_y
-
-
-
-! Have a global variable saying whether this is dipole or regular lon-lat grid
-! This should be initialized static_init_model. Code to do this is below.
-logical :: dipole_grid
-
-contains
-
-!------------------------------------------------------------------
-!------------------------------------------------------------------
-
-subroutine pop_static_init_model()
-
-! Called to do one time initialization of the model. In this case,
-! it reads in the grid information.
-
-integer :: iunit, io
-integer :: ss, dd
-
-! The Plan:
-!
-!   read in the grid sizes from the horiz grid file and the vert grid file
-!   horiz is netcdf, vert is ascii
-!  
-!   allocate space, and read in actual grid values
-!
-!   figure out model timestep.  FIXME: from where?
-!
-!   Compute the model size.
-!
-!   set the index numbers where the field types change
-!
-
-if ( module_initialized ) return ! only need to do this once.
-
-! Print module information to log file and stdout.
-call register_module(source, revision, revdate)
-
-! Since this routine calls other routines that could call this routine
-! we'll say we've been initialized pretty dang early.
-module_initialized = .true.
-
-! Read the DART namelist for this model
-call find_namelist_in_file('input.nml', 'pop_model_nml', iunit)
-read(iunit, nml = pop_model_nml, iostat = io)
-call check_namelist_read(iunit, io, 'pop_model_nml')
-
-! Record the namelist values used for the run
-call error_handler(E_MSG,'static_init_model','pop_model_nml values are',' ',' ',' ')
-if (do_output()) write(logfileunit, nml=pop_model_nml)
-if (do_output()) write(     *     , nml=pop_model_nml)
-
-
-! Set the time step ... causes POP namelists to be read.
-! Ensures model_timestep is multiple of 'ocean_dynamics_timestep'
-
-model_timestep = set_model_time_step()
-
-call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
-
-write(msgstring,*)'assimilation period is ',dd,' days ',ss,' seconds'
-call error_handler(E_MSG,'static_init_model',msgstring,source,revision,revdate)
-
-
-! get data dimensions, then allocate space, then open the files
-! and actually fill in the arrays.
-
-call get_horiz_grid_dims(Nx, Ny)
-call get_vert_grid_dim(Nz)
-
-! Allocate space for grid variables. 
-allocate(ULAT(Nx,Ny), ULON(Nx,Ny), TLAT(Nx,Ny), TLON(Nx,Ny))
-allocate( KMT(Nx,Ny),  KMU(Nx,Ny))
-allocate(  HT(Nx,Ny),   HU(Nx,Ny))
-allocate(     ZC(Nz),      ZG(Nz))
-
-
-! Fill them in.
-! horiz grid initializes ULAT/LON, TLAT/LON as well.
-! kmt initializes HT/HU if present in input file.
-call read_horiz_grid(Nx, Ny, ULAT, ULON, TLAT, TLON)
-call read_topography(Nx, Ny,  KMT,  KMU)
-call read_vert_grid( Nz, ZC, ZG)
-
-if (debug > 2) call write_grid_netcdf() ! DEBUG only
-if (debug > 2) call write_grid_interptest() ! DEBUG only
-
-
-! compute the offsets into the state vector for the start of each
-! different variable type.
-
-! record where in the state vector the data type changes
-! from one type to another, by computing the starting
-! index for each block of data.
-start_index(S_index)     = 1
-start_index(T_index)     = start_index(S_index) + (Nx * Ny * Nz)
-start_index(U_index)     = start_index(T_index) + (Nx * Ny * Nz)
-start_index(V_index)     = start_index(U_index) + (Nx * Ny * Nz)
-start_index(PSURF_index) = start_index(V_index) + (Nx * Ny * Nz)
-
-! in spite of the staggering, all grids are the same size
-! and offset by half a grid cell.  4 are 3D and 1 is 2D.
-!  e.g. S,T,U,V = 256 x 225 x 70
-!  e.g. PSURF = 256 x 225
-
-if (do_output()) write(logfileunit, *) 'Using grid : Nx, Ny, Nz = ', &
-                                                     Nx, Ny, Nz
-if (do_output()) write(     *     , *) 'Using grid : Nx, Ny, Nz = ', &
-                                                     Nx, Ny, Nz
-
-model_size = (n3dfields * (Nx * Ny * Nz)) + (n2dfields * (Nx * Ny))
-if (do_output()) write(*,*) 'model_size = ', model_size
-
-! initialize the pressure array - pressure in bars
-allocate(pressure(Nz))
-call dpth2pres(Nz, ZC, pressure)
-
-! Initialize the interpolation routines
-call init_interp()
-
-end subroutine pop_static_init_model
-
-!------------------------------------------------------------
-
-subroutine init_interp()
-
-! Initializes data structures needed for POP interpolation for
-! either dipole or irregular grid.
-! This should be called at static_init_model time to avoid 
-! having all this temporary storage in the middle of a run.
-
-integer :: i
-
-! Determine whether this is a irregular lon-lat grid or a dipole.
-! Do this by seeing if the lons have the same values at both
-! the first and last latitude row; this is not the case for dipole. 
-
-dipole_grid = .false.
-do i = 1, nx
-   if(ulon(i, 1) /= ulon(i, ny)) then
-      dipole_grid = .true.
-      call init_dipole_interp()
-      return
-   endif
-enddo
-
-end subroutine init_interp
-
-!------------------------------------------------------------
-
-subroutine init_dipole_interp()
-
-! Build the data structure for interpolation for a dipole grid.
-
-! Need a temporary data structure to build this.
-! These arrays keep a list of the x and y indices of dipole quads 
-! that potentially overlap the regular boxes. Need one for the u 
-! and one for the t grid.
-integer :: ureg_list_lon(num_reg_x, num_reg_y, max_reg_list_num)
-integer :: ureg_list_lat(num_reg_x, num_reg_y, max_reg_list_num)
-integer :: treg_list_lon(num_reg_x, num_reg_y, max_reg_list_num)
-integer :: treg_list_lat(num_reg_x, num_reg_y, max_reg_list_num)
-
-
-real(r8) :: u_c_lons(4), u_c_lats(4), t_c_lons(4), t_c_lats(4), pole_row_lon
-integer  :: i, j, k, pindex
-integer  :: reg_lon_ind(2), reg_lat_ind(2), u_total, t_total, u_index, t_index
-logical  :: is_pole
-
-! Begin by finding the quad that contains the pole for the dipole t_grid. 
-! To do this locate the u quad with the pole on its right boundary. This is on
-! the row that is opposite the shifted pole and exactly follows a lon circle.
-pole_x = nx / 2;
-! Search for the row at which the longitude flips over the pole
-pole_row_lon = ulon(pole_x, 1);
-do i = 1, ny
-   pindex = i
-   if(ulon(pole_x, i) /= pole_row_lon) exit
-enddo
-
-! Pole boxes for u have indices pole_x or pole_x-1 and index - 1;
-! (it's right before the flip).
-u_pole_y = pindex - 1;
-
-! Locate the T dipole quad that contains the pole.
-! We know it is in either the same lat quad as the u pole edge or one higher.
-! Figure out if the pole is more or less than halfway along
-! the u quad edge to find the right one.
-if(ulat(pole_x, u_pole_y) > ulat(pole_x, u_pole_y + 1)) then
-   t_pole_y = u_pole_y;
-else
-   t_pole_y = u_pole_y + 1;
-endif
-
-! Loop through each of the dipole grid quads 
-do i = 1, nx
-   ! There's no wraparound in y, one box less than grid boundaries
-   do j = 1, ny - 1
-      ! Is this the pole quad for the T grid?
-      is_pole = (i == pole_x .and. j == t_pole_y)
-      
-      ! Set up array of lons and lats for the corners of these u and t quads
-      call get_quad_corners(ulon, i, j, u_c_lons)
-      call get_quad_corners(ulat, i, j, u_c_lats)
-      call get_quad_corners(tlon, i, j, t_c_lons)
-      call get_quad_corners(tlat, i, j, t_c_lats)
-
-      ! Get list of regular boxes that cover this u dipole quad
-      ! false indicates that for the u grid there's nothing special about pole
-      call reg_box_overlap(u_c_lons, u_c_lats, .false., reg_lon_ind, reg_lat_ind)         
-
-      ! Update the temporary data structures for the u quad 
-      call update_reg_list(u_dipole_num, ureg_list_lon, &
-         ureg_list_lat, reg_lon_ind, reg_lat_ind, i, j)
-
-      ! Repeat for t dipole quads
-      call reg_box_overlap(t_c_lons, t_c_lats, is_pole, reg_lon_ind, reg_lat_ind)         
-      call update_reg_list(t_dipole_num, treg_list_lon, &
-         treg_list_lat, reg_lon_ind, reg_lat_ind, i, j)
-   enddo
-enddo
-
-if (do_output()) write(*,*)'to determine (minimum) max_reg_list_num values for new grids ...'
-if (do_output()) write(*,*)'u_dipole_num is ',maxval(u_dipole_num)
-if (do_output()) write(*,*)'t_dipole_num is ',maxval(t_dipole_num)
-
-! Invert the temporary data structure. The total number of entries will be 
-! the sum of the number of dipole cells for each regular cell. 
-u_total = sum(u_dipole_num)
-t_total = sum(t_dipole_num)
-
-! Allocate storage for the final structures in module storage
-allocate(u_dipole_lon_list(u_total), u_dipole_lat_list(u_total))
-allocate(t_dipole_lon_list(t_total), t_dipole_lat_list(t_total))
-
-! Fill up the long list by traversing the temporary structure. Need indices 
-! to keep track of where to put the next entry.
-u_index = 1
-t_index = 1
-
-! Loop through each regular grid box
-do i = 1, num_reg_x
-   do j = 1, num_reg_y
-
-      ! The list for this regular box starts at the current indices.
-      u_dipole_start(i, j) = u_index
-      t_dipole_start(i, j) = t_index
-
-      ! Copy all the close dipole quads for regular u box(i, j)
-      do k = 1, u_dipole_num(i, j)
-         u_dipole_lon_list(u_index) = ureg_list_lon(i, j, k) 
-         u_dipole_lat_list(u_index) = ureg_list_lat(i, j, k) 
-         u_index = u_index + 1
-      enddo
-      
-      ! Copy all the close dipoles for regular t box (i, j)
-      do k = 1, t_dipole_num(i, j)
-         t_dipole_lon_list(t_index) = treg_list_lon(i, j, k) 
-         t_dipole_lat_list(t_index) = treg_list_lat(i, j, k) 
-         t_index = t_index + 1
-      enddo
-
-   enddo
-enddo
-
-! Confirm that the indices come out okay as debug
-if(u_index /= u_total + 1) then
-   msgstring = 'Storage indices did not balance for U grid: : contact DART developers'
-   call error_handler(E_ERR, 'init_dipole_interp', msgstring, source, revision, revdate)
-endif
-if(t_index /= t_total + 1) then
-   msgstring = 'Storage indices did not balance for T grid: : contact DART developers'
-   call error_handler(E_ERR, 'init_dipole_interp', msgstring, source, revision, revdate)
-endif
-
-end subroutine init_dipole_interp
-
-!------------------------------------------------------------
-
-subroutine get_reg_box_indices(lon, lat, x_ind, y_ind)
- real(r8), intent(in)  :: lon, lat
- integer,  intent(out) :: x_ind, y_ind
-
-! Given a longitude and latitude in degrees returns the index of the regular
-! lon-lat box that contains the point.
-
-call get_reg_lon_box(lon, x_ind)
-call get_reg_lat_box(lat, y_ind)
-
-end subroutine get_reg_box_indices
-
-!------------------------------------------------------------
-
-subroutine get_reg_lon_box(lon, x_ind)
- real(r8), intent(in)  :: lon
- integer,  intent(out) :: x_ind
-
-! Determine which regular longitude box a longitude is in.
-
-x_ind = int(num_reg_x * lon / 360.0_r8) + 1
-
-! Watch out for exactly at top; assume all lats and lons in legal range
-if(lon == 360.0_r8) x_ind = num_reg_x
-
-end subroutine get_reg_lon_box
-
-!------------------------------------------------------------
-
-subroutine get_reg_lat_box(lat, y_ind)
- real(r8), intent(in)  :: lat
- integer,  intent(out) :: y_ind
-
-! Determine which regular latitude box a latitude is in.
-
-y_ind = int(num_reg_y * (lat + 90.0_r8) / 180.0_r8) + 1
-
-! Watch out for exactly at top; assume all lats and lons in legal range
-if(lat == 90.0_r8)  y_ind = num_reg_y
-
-end subroutine get_reg_lat_box
-
-!------------------------------------------------------------
-
-subroutine reg_box_overlap(x_corners, y_corners, is_pole, reg_lon_ind, reg_lat_ind)
- real(r8), intent(in)  :: x_corners(4), y_corners(4)
- logical,  intent(in)  :: is_pole
- integer,  intent(out) :: reg_lon_ind(2), reg_lat_ind(2)
-
-! Find a set of regular lat lon boxes that covers all of the area covered by 
-! a dipole grid qaud whose corners are given by the dimension four x_corners 
-! and y_corners arrays.  The two dimensional arrays reg_lon_ind and reg_lat_ind
-! return the first and last indices of the regular boxes in latitude and
-! longitude respectively. These indices may wraparound for reg_lon_ind.  
-! A special computation is needed for a dipole quad that has the true north 
-! pole in its interior. The logical is_pole is set to true if this is the case.
-! This can only happen for the t grid.  If the longitude boxes overlap 0
-! degrees, the indices returned are adjusted by adding the total number of
-! boxes to the second index (e.g. the indices might be 88 and 93 for a case
-! with 90 longitude boxes).
-
-real(r8) :: lat_min, lat_max, lon_min, lon_max
-integer  :: i
-
-!  A quad containing the pole is fundamentally different
-if(is_pole) then
-   ! Need all longitude boxes
-   reg_lon_ind(1) = 1
-   reg_lon_ind(2) = num_reg_x
-   ! Need to cover from lowest latitude to top box
-   lat_min = minval(y_corners)
-   reg_lat_ind(1) = int(num_reg_y * (lat_min + 90.0_r8) / 180.0_r8) + 1
-   call get_reg_lat_box(lat_min, reg_lat_ind(1))
-   reg_lat_ind(2) = num_reg_y
-else
-   ! All other quads do not contain pole (pole could be on edge but no problem)
-   ! This is specific to the dipole POP grids that do not go to the south pole
-   ! Finding the range of latitudes is cake
-   lat_min = minval(y_corners)
-   lat_max = maxval(y_corners)
-
-   ! Figure out the indices of the regular boxes for min and max lats
-   call get_reg_lat_box(lat_min, reg_lat_ind(1))
-   call get_reg_lat_box(lat_max, reg_lat_ind(2))
-
-   ! Lons are much trickier. Need to make sure to wraparound the
-   ! right way. There is no guarantee on direction of lons in the
-   ! high latitude dipole rows.
-   ! All longitudes for non-pole rows have to be within 180 degrees
-   ! of one another. 
-   lon_min = minval(x_corners)
-   lon_max = maxval(x_corners)
-   if((lon_max - lon_min) > 180.0_r8) then
-      ! If the max longitude value is more than 180 
-      ! degrees larger than the min, then there must be wraparound.
-      ! Then, find the smallest value > 180 and the largest < 180 to get range.
-      lon_min = 360.0_r8
-      lon_max = 0.0_r8
-      do i=1, 4
-         if(x_corners(i) > 180.0_r8 .and. x_corners(i) < lon_min) lon_min = x_corners(i)
-         if(x_corners(i) < 180.0_r8 .and. x_corners(i) > lon_max) lon_max = x_corners(i)
-      enddo
-   endif
-
-   ! Get the indices for the extreme longitudes
-   call get_reg_lon_box(lon_min, reg_lon_ind(1))
-   call get_reg_lon_box(lon_max, reg_lon_ind(2))
-
-   ! Watch for wraparound again; make sure that second index is greater than first
-   if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
-endif
-
-end subroutine reg_box_overlap
-
-!------------------------------------------------------------
-
-subroutine get_quad_corners(x, i, j, corners)
- real(r8), intent(in)  :: x(:, :)
- integer,  intent(in)  :: i, j
- real(r8), intent(out) :: corners(4)
-
-! Grabs the corners for a given quadrilateral from the global array of lower
-! right corners. Note that corners go counterclockwise around the quad.
-
-integer :: ip1
-
-! Have to worry about wrapping in longitude but not in latitude
-ip1 = i + 1
-if(ip1 > nx) ip1 = 1
-
-corners(1) = x(i,   j  ) 
-corners(2) = x(ip1, j  )
-corners(3) = x(ip1, j+1)
-corners(4) = x(i,   j+1)
-
-end subroutine get_quad_corners
-
-!------------------------------------------------------------
-
-subroutine update_reg_list(reg_list_num, reg_list_lon, reg_list_lat, &
-                           reg_lon_ind, reg_lat_ind, dipole_lon_index, dipole_lat_index)
-
- integer, intent(inout) :: reg_list_num(:, :), reg_list_lon(:, :, :), reg_list_lat(:, :, :)
- integer, intent(inout) :: reg_lon_ind(2), reg_lat_ind(2)
- integer, intent(in)    :: dipole_lon_index, dipole_lat_index
- 
-! Updates the data structure listing dipole quads that are in a given regular box
-integer :: ind_x, index_x, ind_y
-
-! Loop through indices for each possible regular cell
-! Have to watch for wraparound in longitude
-if(reg_lon_ind(2) < reg_lon_ind(1)) reg_lon_ind(2) = reg_lon_ind(2) + num_reg_x
-
-do ind_x = reg_lon_ind(1), reg_lon_ind(2)
-   ! Inside loop, need to go back to wraparound indices to find right box
-   index_x = ind_x
-   if(index_x > num_reg_x) index_x = index_x - num_reg_x
-   
-   do ind_y = reg_lat_ind(1), reg_lat_ind(2)
-      ! Make sure the list storage isn't full
-      if(reg_list_num(index_x, ind_y) >= max_reg_list_num) then
-         write(msgstring,*) 'max_reg_list_num (',max_reg_list_num,') is too small ... increase'
-         call error_handler(E_ERR, 'update_reg_list', msgstring, source, revision, revdate)
-      endif
-
-      ! Increment the count
-      reg_list_num(index_x, ind_y) = reg_list_num(index_x, ind_y) + 1
-      ! Store this quad in the list for this regular box
-      reg_list_lon(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lon_index
-      reg_list_lat(index_x, ind_y, reg_list_num(index_x, ind_y)) = dipole_lat_index
-   enddo
-enddo
-
-end subroutine update_reg_list
-
-!------------------------------------------------------------
-
-subroutine pop_init_conditions(x)
- real(r8), intent(out) :: x(:)
-
-! Returns a model state vector, x, that is some sort of appropriate
-! initial condition for starting up a long integration of the model.
-! At present, this is only used if the namelist parameter 
-! start_from_restart is set to .false. in the program perfect_model_obs.
-
-character(len=128) :: msgstring2, msgstring3
-
-msgstring2 = "cannot run perfect_model_obs with 'start_from_restart = .false.' "
-msgstring3 = 'use pop_to_dart to generate an initial state'
-call error_handler(E_ERR,'init_conditions', &
-                  'WARNING!!  POP model has no built-in default state', &
-                  source, revision, revdate, &
-                  text2=msgstring2, text3=msgstring3)
-
-! this code never reached - just here to avoid compiler warnings
-! about an intent(out) variable not being set to a value.
-x = 0.0_r8
-
-end subroutine pop_init_conditions
-
-!------------------------------------------------------------------
-
-subroutine pop_adv_1step(x, time)
- real(r8),        intent(inout) :: x(:)
- type(time_type), intent(in)    :: time
-
-! If the model could be called as a subroutine, does a single
-! timestep advance.  POP cannot be called this way, so fatal error
-! if this routine is called.
-
-call error_handler(E_ERR,'adv_1step', &
-                  'POP model cannot be called as a subroutine; async cannot = 0', &
-                  source, revision, revdate)
-
-end subroutine pop_adv_1step
-
-!------------------------------------------------------------------
-
-function pop_get_model_size()
- integer :: pop_get_model_size
-
-! Returns the size of the model as an integer. Required for all
-! applications.
-
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-pop_get_model_size = model_size
-
-end function pop_get_model_size
-
-!------------------------------------------------------------------
-
-subroutine pop_init_time(time)
- type(time_type), intent(out) :: time
-
-! Companion interface to init_conditions. Returns a time that is
-! appropriate for starting up a long integration of the model.
-! At present, this is only used if the namelist parameter 
-! start_from_restart is set to .false. in the program perfect_model_obs.
-
-character(len=128) :: msgstring2, msgstring3
-
-msgstring2 = "cannot run perfect_model_obs with 'start_from_restart = .false.' "
-msgstring3 = 'use pop_to_dart to generate an initial state which contains a timestamp'
-call error_handler(E_ERR,'init_time', &
-                  'WARNING!!  POP model has no built-in default time', &
-                  source, revision, revdate, &
-                  text2=msgstring2, text3=msgstring3)
-
-! this code never reached - just here to avoid compiler warnings
-! about an intent(out) variable not being set to a value.
-time = set_time(0,0)
-
-end subroutine pop_init_time
-
-!------------------------------------------------------------------
-
-subroutine pop_model_interpolate(x, location, obs_type, interp_val, istatus)
- real(r8),            intent(in) :: x(:)
- type(location_type), intent(in) :: location
- integer,             intent(in) :: obs_type
- real(r8),           intent(out) :: interp_val
- integer,            intent(out) :: istatus
-
-! Model interpolate will interpolate any state variable (S, T, U, V, PSURF) to
-! the given location given a state vector. The type of the variable being
-! interpolated is obs_type since normally this is used to find the expected
-! value of an observation at some location. The interpolated value is 
-! returned in interp_val and istatus is 0 for success.
-
-! Local storage
-real(r8)       :: loc_array(3), llon, llat, lheight
-integer        :: base_offset, offset, ind
-integer        :: hgt_bot, hgt_top
-real(r8)       :: hgt_fract
-real(r8)       :: top_val, bot_val
-integer        :: hstatus
-logical        :: convert_to_ssh
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! print data min/max values
-if (debug > 2) call print_ranges(x)
-
-! Let's assume failure.  Set return val to missing, then the code can
-! just set istatus to something indicating why it failed, and return.
-! If the interpolation is good, the interp_val will be set to the 
-! good value, and the last line here sets istatus to 0.
-! make any error codes set here be in the 10s
-
-interp_val = MISSING_R8     ! the DART bad value flag
-istatus = 99                ! unknown error
-
-! Get the individual locations values
-loc_array = get_location(location)
-llon    = loc_array(1)
-llat    = loc_array(2)
-lheight = loc_array(3)
-
-if (debug > 1) print *, 'requesting interpolation of ', obs_type, ' at ', llon, llat, lheight
-
-if( vert_is_height(location) ) then
-   ! Nothing to do 
-elseif ( vert_is_surface(location) ) then
-   ! Nothing to do 
-elseif (vert_is_level(location)) then
-   ! convert the level index to an actual depth 
-   ind = nint(loc_array(3))
-   if ( (ind < 1) .or. (ind > size(zc)) ) then 
-      istatus = 11
-      return
-   else
-      lheight = zc(ind)
-   endif
-else   ! if pressure or undefined, we don't know what to do
-   istatus = 17
-   return
-endif
-
-! kind (in-situ) temperature is a combination of potential temp,
-! salinity, and pressure based on depth.  call a routine that
-! interpolates all three, does the conversion, and returns the
-! sensible/in-situ temperature.
-if(obs_type == KIND_TEMPERATURE) then
-   ! we know how to interpolate this from potential temp,
-   ! salinity, and pressure based on depth.
-   call compute_temperature(x, llon, llat, lheight, interp_val, istatus)
-   if (debug > 1) print *, 'interp val, istatus = ', interp_val, istatus
-   return
-endif
-
-
-! The following kinds are either in the state vector (so you
-! can simply interpolate to find the value) or they are a simple
-! transformation of something in the state vector.
-
-convert_to_ssh = .FALSE.
-
-if(obs_type == KIND_SALINITY) then
-   base_offset = start_index(S_index)
-else if(obs_type == KIND_POTENTIAL_TEMPERATURE) then
-   base_offset = start_index(T_index)
-else if(obs_type == KIND_U_CURRENT_COMPONENT) then
-   base_offset = start_index(U_index)
-else if(obs_type == KIND_V_CURRENT_COMPONENT) then
-   base_offset = start_index(V_index)
-else if(obs_type == KIND_SEA_SURFACE_PRESSURE) then
-   base_offset = start_index(PSURF_index)
-else if(obs_type == KIND_SEA_SURFACE_HEIGHT) then
-   base_offset = start_index(PSURF_index) ! simple linear transform of PSURF
-   convert_to_ssh = .TRUE.
-else
-   ! Not a legal type for interpolation, return istatus error
-   istatus = 15
-   return
-endif
-
-! For Sea Surface Height or Pressure don't need the vertical coordinate
-! SSP needs to be converted to a SSH if height is required.
-if( vert_is_surface(location) ) then
-   call lon_lat_interpolate(x(base_offset:), llon, llat, obs_type, 1, interp_val, istatus)
-   if (convert_to_ssh .and. (istatus == 0)) then
-      interp_val = interp_val / 980.6_r8   ! POP uses CGS units
-   endif
-   if (debug > 1) print *, 'interp val, istatus = ', interp_val, istatus
-   return
-endif
-
-! Get the bounding vertical levels and the fraction between bottom and top
-call height_bounds(lheight, Nz, ZC, hgt_bot, hgt_top, hgt_fract, hstatus)
-if(hstatus /= 0) then
-   istatus = 12
-   return
-endif
-
-! do a 2d interpolation for the value at the bottom level, then again for
-! the top level, then do a linear interpolation in the vertical to get the
-! final value.  this sets both interp_val and istatus.
-call do_interp(x, base_offset, hgt_bot, hgt_top, hgt_fract, &
-               llon, llat, obs_type, interp_val, istatus)
-if (debug > 1) print *, 'interp val, istatus = ', interp_val, istatus
-
-end subroutine pop_model_interpolate
-
-!------------------------------------------------------------------
-
-! Three different types of grids are used here. The POP dipole 
-! grid is referred to as a dipole grid and each region is
-! referred to as a quad, short for quadrilateral. 
-! The longitude latitude rectangular grid with possibly irregular
-! spacing in latitude used for some POP applications and testing
-! is referred to as the irregular grid and each region is 
-! called a box.
-! Finally, a regularly spaced longitude latitude grid is used
-! as a computational tool for interpolating from the dipole
-! grid. This is referred to as the regular grid and each region
-! is called a box. 
-! All grids are referenced by the index of the lower left corner
-! of the quad or box.
-
-! The dipole grid is assumed to be global for all applications.
-! The irregular grid is also assumed to be global east
-! west for all applications.
-
-!------------------------------------------------------------------
-
-subroutine lon_lat_interpolate(x, lon, lat, var_type, height, interp_val, istatus)
- real(r8), intent(in) :: x(:)
- real(r8), intent(in) :: lon, lat
- integer,  intent(in) :: var_type, height
- real(r8), intent(out) :: interp_val
- integer,  intent(out) :: istatus
-
-! Subroutine to interpolate to a lon lat location given the state vector 
-! for that level, x. This works just on one horizontal slice.
-! NOTE: Using array sections to pass in the x array may be inefficient on some
-! compiler/platform setups. Might want to pass in the entire array with a base
-! offset value instead of the section if this is an issue.
-! This routine works for either the dipole or a regular lat-lon grid.
-! Successful interpolation returns istatus=0.
-
-! Local storage
-integer  :: lat_bot, lat_top, lon_bot, lon_top, num_inds, start_ind
-integer  :: x_ind, y_ind
-real(r8) :: p(4), x_corners(4), y_corners(4), xbot, xtop
-real(r8) :: lon_fract, lat_fract
-logical  :: masked
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! Succesful return has istatus of 0
-istatus = 0
-
-! Get the lower left corner for either grid type
-if(dipole_grid) then
-   ! Figure out which of the regular grid boxes this is in
-   call get_reg_box_indices(lon, lat, x_ind, y_ind)
-
-   ! Is this on the U or T grid?
-   if(is_on_ugrid(var_type)) then
-      ! On U grid
-      num_inds =  u_dipole_num  (x_ind, y_ind)
-      start_ind = u_dipole_start(x_ind, y_ind)
-
-      ! If there are no quads overlapping, can't do interpolation
-      if(num_inds == 0) then
-         istatus = 1
-         return
-      endif
-
-      ! Search the list of quads to see if (lon, lat) is in one
-      call get_dipole_quad(lon, lat, ulon, ulat, num_inds, start_ind, &
-         u_dipole_lon_list, u_dipole_lat_list, lon_bot, lat_bot, istatus)
-      ! Fail on bad istatus return
-      if(istatus /= 0) return
-
-      ! Getting corners for accurate interpolation
-      call get_quad_corners(ulon, lon_bot, lat_bot, x_corners)
-      call get_quad_corners(ulat, lon_bot, lat_bot, y_corners)
-
-      ! Fail if point is in one of the U boxes that go through the
-      ! pole (this could be fixed up if necessary)
-      if(lat_bot == u_pole_y .and. (lon_bot == pole_x -1 .or. &
-         lon_bot == pole_x)) then
-         istatus = 4
-         return
-      endif
-
-   else
-      ! On T grid
-      num_inds =  t_dipole_num  (x_ind, y_ind)
-      start_ind = t_dipole_start(x_ind, y_ind)
-      call get_dipole_quad(lon, lat, tlon, tlat, num_inds, start_ind, &
-         t_dipole_lon_list, t_dipole_lat_list, lon_bot, lat_bot, istatus)
-      ! Fail on bad istatus return
-      if(istatus /= 0) return
-
-      ! Fail if point is in T box that covers pole
-      if(lon_bot == pole_x .and. lat_bot == t_pole_y) then
-         istatus = 5
-         return
-      endif
-
-      ! Getting corners for accurate interpolation
-      call get_quad_corners(tlon, lon_bot, lat_bot, x_corners)
-      call get_quad_corners(tlat, lon_bot, lat_bot, y_corners)
-   endif
-
-else
-   ! This is an irregular grid
-   ! U and V are on velocity grid
-   if (is_on_ugrid(var_type)) then
-      ! Get the corner indices and the fraction of the distance between
-      call get_irreg_box(lon, lat, ulon, ulat, &
-         lon_bot, lat_bot, lon_fract, lat_fract, istatus)
-   else
-      ! Eta, T and S are on the T grid
-      ! Get the corner indices
-      call get_irreg_box(lon, lat, tlon, tlat, &
-         lon_bot, lat_bot, lon_fract, lat_fract, istatus)
-   endif
-
-   ! Return passing through error status
-   if(istatus /= 0) return
-
-endif
-
-! Find the indices to get the values for interpolating
-lat_top = lat_bot + 1
-if(lat_top > ny) then
-   istatus = 2
-   return
-endif
-
-! Watch for wraparound in longitude
-lon_top = lon_bot + 1
-if(lon_top > nx) lon_top = 1
-
-! Get the values at the four corners of the box or quad
-! Corners go around counterclockwise from lower left
-p(1) = get_val(lon_bot, lat_bot, nx, x, var_type, height, masked)
-if(masked) then
-   istatus = 3
-   return
-endif
-
-p(2) = get_val(lon_top, lat_bot, nx, x, var_type, height, masked)
-if(masked) then
-   istatus = 3
-   return
-endif
-
-p(3) = get_val(lon_top, lat_top, nx, x, var_type, height, masked)
-if(masked) then
-   istatus = 3
-   return
-endif
-
-p(4) = get_val(lon_bot, lat_top, nx, x, var_type, height, masked)
-if(masked) then
-   istatus = 3
-   return
-endif
-
-! Full bilinear interpolation for quads
-if(dipole_grid) then
-   call quad_bilinear_interp(lon, lat, x_corners, y_corners, p, interp_val)
-else
-   ! Rectangular biliear interpolation
-   xbot = p(1) + lon_fract * (p(2) - p(1))
-   xtop = p(4) + lon_fract * (p(3) - p(4))
-   ! Now interpolate in latitude
-   interp_val = xbot + lat_fract * (xtop - xbot)
-endif
-
-end subroutine lon_lat_interpolate
-
-!------------------------------------------------------------
-
-function get_val(lon_index, lat_index, nlon, x, var_type, height, masked)
- integer,     intent(in) :: lon_index, lat_index, nlon, var_type, height
- real(r8),    intent(in) :: x(:)
- logical,    intent(out) :: masked
- real(r8)                :: get_val
-
-! Returns the value from a single level array given the lat and lon indices
-! 'masked' returns true if this is NOT a valid grid location (e.g. land, or
-! below the ocean floor in shallower areas).
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! check the land/ocean bottom map and return if not valid water cell.
-if(is_dry_land(var_type, lon_index, lat_index, height)) then
-   masked = .true.
-   get_val = MISSING_R8
-   return
-endif
-
-! Layout has lons varying most rapidly
-get_val = x((lat_index - 1) * nlon + lon_index)
-
-! this is a valid ocean water cell, not land or below ocean floor
-masked = .false.
-
-end function get_val
-
-!------------------------------------------------------------
-
-subroutine get_irreg_box(lon, lat, lon_array, lat_array, &
-                         found_x, found_y, lon_fract, lat_fract, istatus)
-
- real(r8),   intent(in) :: lon, lat
- real(r8),   intent(in) :: lon_array(nx, ny), lat_array(nx, ny)
- real(r8),  intent(out) :: lon_fract, lat_fract
- integer,   intent(out) :: found_x, found_y, istatus
-
-! Given a longitude and latitude of a point (lon and lat) and the
-! longitudes and latitudes of the lower left corner of the regular grid
-! boxes, gets the indices of the grid box that contains the point and
-! the fractions along each directrion for interpolation.
-
-! Local storage
-integer  :: lat_status, lon_top, lat_top
-
-! Succesful return has istatus of 0
-istatus = 0
-
-! Get latitude box boundaries
-call lat_bounds(lat, ny, lat_array, found_y, lat_top, lat_fract, lat_status)
-
-! Check for error on the latitude interpolation
-if(lat_status /= 0) then
-   istatus = 1
-   return
-endif
-
-! Find out what longitude box and fraction
-call lon_bounds(lon, nx, lon_array, found_x, lon_top, lon_fract)
-
-end subroutine get_irreg_box
-
-!------------------------------------------------------------
-
-subroutine lon_bounds(lon, nlons, lon_array, bot, top, fract)
- real(r8),    intent(in) :: lon
- integer,     intent(in) :: nlons
- real(r8),    intent(in) :: lon_array(:, :)
- integer,    intent(out) :: bot, top
- real(r8),   intent(out) :: fract
-
-! Given a longitude lon, the array of longitudes for grid boundaries, and the
-! number of longitudes in the grid, returns the indices of the longitude
-! below and above the location longitude and the fraction of the distance
-! between. It is assumed that the longitude wraps around for a global grid. 
-! Since longitude grids are going to be regularly spaced, this could be made more efficient.
-! Algorithm fails for a silly grid that has only two longitudes separated by 180 degrees.
-
-! Local storage
-integer  :: i
-real(r8) :: dist_bot, dist_top
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! This is inefficient, someone could clean it up since longitudes are regularly spaced
-! But note that they don't have to start at 0
-do i = 2, nlons
-   dist_bot = lon_dist(lon, lon_array(i - 1, 1))
-   dist_top = lon_dist(lon, lon_array(i, 1))
-   if(dist_bot <= 0 .and. dist_top > 0) then
-      bot = i - 1
-      top = i
-      fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
-      return
-   endif
-enddo
-
-! Falling off the end means it's in between; wraparound
-bot = nlons
-top = 1
-dist_bot = lon_dist(lon, lon_array(bot, 1))
-dist_top = lon_dist(lon, lon_array(top, 1)) 
-fract = abs(dist_bot) / (abs(dist_bot) + dist_top)
-
-end subroutine lon_bounds
-
-!-------------------------------------------------------------
-
-subroutine lat_bounds(lat, nlats, lat_array, bot, top, fract, istatus)
- real(r8),   intent(in) :: lat
- integer,    intent(in) :: nlats
- real(r8),   intent(in) :: lat_array(:, :)
- integer,   intent(out) :: bot, top
- real(r8),  intent(out) :: fract
- integer,   intent(out) :: istatus
-
-! Given a latitude lat, the array of latitudes for grid boundaries, and the
-! number of latitudes in the grid, returns the indices of the latitude
-! below and above the location latitude and the fraction of the distance
-! between. istatus is returned as 0 unless the location latitude is 
-! south of the southernmost grid point (1 returned) or north of the 
-! northernmost (2 returned). If one really had lots of polar obs would 
-! want to worry about interpolating around poles.
-
-! Local storage
-integer :: i
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! Success should return 0, failure a positive number.
-istatus = 0
-
-! Check for too far south or north
-if(lat < lat_array(1, 1)) then
-   istatus = 1
-   return
-else if(lat > lat_array(1, nlats)) then
-   istatus = 2
-   return
-endif
-
-! In the middle, search through
-do i = 2, nlats
-   if(lat <= lat_array(1, i)) then
-      bot = i - 1
-      top = i
-      fract = (lat - lat_array(1, bot)) / (lat_array(1, top) - lat_array(1, bot))
-      return
-   endif
-enddo
-
-! Shouldn't get here. Might want to fail really hard through error handler
-istatus = 40
-
-end subroutine lat_bounds
-
-!------------------------------------------------------------------
-
-function lon_dist(lon1, lon2)
- real(r8), intent(in) :: lon1, lon2
- real(r8)             :: lon_dist
-
-! Returns the smallest signed distance between lon1 and lon2 on the sphere
-! If lon1 is less than 180 degrees east of lon2 the distance is negative
-! If lon1 is less than 180 degrees west of lon2 the distance is positive
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-lon_dist = lon2 - lon1
-if(lon_dist >= -180.0_r8 .and. lon_dist <= 180.0_r8) then
-   return
-else if(lon_dist < -180.0_r8) then
-   lon_dist = lon_dist + 360.0_r8
-else
-   lon_dist = lon_dist - 360.0_r8
-endif
-
-end function lon_dist
-
-!------------------------------------------------------------
-
-subroutine get_dipole_quad(lon, lat, qlons, qlats, num_inds, start_ind, &
-                           x_inds, y_inds, found_x, found_y, istatus)
-
- real(r8), intent(in)  :: lon, lat, qlons(:, :), qlats(:, :)
- integer,  intent(in)  :: num_inds, start_ind, x_inds(:), y_inds(:)
- integer,  intent(out) :: found_x, found_y, istatus
-
-! Given the lon and lat of a point, and a list of the
-! indices of the quads that might contain a point at (lon, lat), determines
-! which quad contains the point.  istatus is returned as 0 if all went 
-! well and 1 if the point was not found to be in any of the quads.
-
-integer :: i, my_index
-real(r8) :: x_corners(4), y_corners(4)
-
-istatus = 0
-
-! Loop through all the quads and see if the point is inside
-do i = 1, num_inds
-   my_index = start_ind + i - 1
-   call get_quad_corners(qlons, x_inds(my_index), y_inds(my_index), x_corners)
-   call get_quad_corners(qlats, x_inds(my_index), y_inds(my_index), y_corners)
-
-   ! Ssearch in this individual quad
-   if(in_quad(lon, lat, x_corners, y_corners)) then
-      found_x = x_inds(my_index)
-      found_y = y_inds(my_index)
-      return
-   endif
-enddo
-
-! Falling off the end means search failed, return istatus 1
-istatus = 1
-
-end subroutine get_dipole_quad
-
-!------------------------------------------------------------
-
-function in_quad(lon, lat, x_corners, y_corners)
- real(r8), intent(in)  :: lon, lat, x_corners(4), y_corners(4)
- logical               :: in_quad
-
-! Return in_quad true if the point (lon, lat) is in the quad with 
-! the given corners.
-
-! Do this by line tracing in latitude for now. For non-pole point, want a vertical
-! line from the lon, lat point to intersect a side of the quad both above
-! and below the point.
-
-real(r8) :: x(2), y(2)
-logical  :: cant_be_in_box, in_box
-integer  :: intercepts_above(4), intercepts_below(4), i
-integer  :: num_above, num_below
-
-! Default answer is point is not in quad
-in_quad = .false.
-
-! Loop through the sides and compute intercept (if any) with a vertical line
-! from the point. This line has equation x=lon.
-do i = 1, 4
-   ! Load up the sides endpoints
-   if(i <= 3) then
-      x(1:2) = x_corners(i:i+1)
-      y(1:2) = y_corners(i:i+1)
-   else
-      x(1) = x_corners(4)
-      x(2) = x_corners(1)
-      y(1) = y_corners(4)
-      y(2) = y_corners(1)
-   endif
-
-   ! Check to see how a vertical line from the point is related to this side
-   call line_intercept(x, y, lon, lat, cant_be_in_box, in_box, intercepts_above(i), &
-      intercepts_below(i))
-
-   ! If cant_be_in_box is true, can return right away
-   if(cant_be_in_box) then
-      in_quad = .false.
-      return
-   ! Return true if on a side
-   else if(in_box) then
-      in_quad = .true.
-      return
-   endif
-
-enddo
-
-! See if the line intercepted a side of the quad both above and below
-num_above = sum(intercepts_above)
-num_below = sum(intercepts_below)
-
-if(num_above > 0 .and. num_below > 0) then
-   in_quad = .true.
-endif
-
-end function in_quad
-
-!------------------------------------------------------------
-
-subroutine line_intercept(side_x_in, side_y, x_point_in, y_point, &
-                          cant_be_in_box, in_box, intercept_above, intercept_below)
-
- real(r8), intent(in)  :: side_x_in(2), side_y(2), x_point_in, y_point
- logical,  intent(out) :: cant_be_in_box, in_box
- integer,  intent(out) :: intercept_above, intercept_below
-
-! Find the intercept of a vertical line from point (x_point, y_point) and
-! a line segment with endpoints side_x and side_y.
-! For a given side have endpoints (side_x1, side_y1) and (side_x2, side_y2)
-! so equation of segment is y = side_y1 + m(x-side_x1) for y 
-! between side_y1 and side_y2.
-! Intersection of vertical line and line containing side 
-! occurs at y = side_y1 + m(x_point - side_x1); need this
-! y to be between side_y1 and side_y2.
-! If the vertical line is colinear with the side but the point is not on the side, return
-! cant_be_in_box as true. If the point is on the side, return in_box true.
-! If the intersection of the vertical line and the side occurs at a point above
-! the given point, return 1 for intercept_above. If the intersection occurs 
-! below, return 1 for intercept_below. If the vertical line does not intersect
-! the segment, return false and 0 for all intent out arguments.
-
-! WARNING: CERTAINLY A PROBLEM FOR THE POLE BOX!!! POLE BOX COULD
-! HAVE SIDES THAT ARE LONGER THAN 180. For now pole boxes are excluded.
-
-! This can probably be made much cleaner and more efficient.
-
-real(r8) :: slope, y_intercept, side_x(2), x_point
-
-! May have to adjust the longitude intent in values, so copy
-side_x = side_x_in
-x_point = x_point_in
-
-! See if the side wraps around in longitude
-if(maxval(side_x) - minval(side_x) > 180.0_r8) then
-   if(side_x(1) < 180.0_r8)  side_x(1) =  side_x(1) + 360.0_r8
-   if(side_x(2) < 180.0_r8)  side_x(2) =  side_x(2) + 360.0_r8
-   if(x_point < 180.0_r8) x_point = x_point + 360.0_r8
-endif
-
-! Initialize the default returns 
-cant_be_in_box   = .false.
-in_box           = .false.
-intercept_above = 0
-intercept_below = 0
-
-! First easy check, if x_point is not between endpoints of segment doesn't intersect
-if(x_point < minval(side_x) .or. x_point > maxval(side_x)) return
-
-! Otherwise line must intersect the segment
-
-! First subblock, slope is undefined
-if(side_x(2) == side_x(1)) then
-   ! The line is colinear with the side
-   ! If y_point is between endpoints then point is on this side
-   if(y_point <= maxval(side_y) .and. y_point >= minval(side_y)) then
-      in_box = .true.
-      return
-   ! If not on side but colinear with side, point cant be in quad
-   else
-      cant_be_in_box = .true.
-      return
-   endif
-
-else
-
-   ! Second possibility; slope is defined
-   ! FIXME: watch out for numerical instability.
-   ! near-zero x's and large side_y's may cause overflow
-   slope = (side_y(2) - side_y(1)) / (side_x(2) - side_x(1))
-
-   ! Intercept of vertical line through is at x_point and...
-   y_intercept = side_y(1) + slope * (x_point - side_x(1))
-
-   ! Intersects the segment, is it above, below, or at the point
-   if(y_intercept == y_point) then
-      in_box = .true.
-      return
-   else if(y_intercept > y_point) then
-      intercept_above = 1
-      return
-   else
-      intercept_below = 1
-      return
-   endif
-endif
-
-end subroutine line_intercept
-
-!------------------------------------------------------------
-
-subroutine quad_bilinear_interp(lon_in, lat, x_corners_in, y_corners, &
-                                p, interp_val)
-
- real(r8),  intent(in) :: lon_in, lat, x_corners_in(4), y_corners(4), p(4)
- real(r8), intent(out) :: interp_val
-
-! Given a longitude and latitude (lon_in, lat), the longitude and
-! latitude of the 4 corners of a quadrilateral and the values at the
-! four corners, interpolates to (lon_in, lat) which is assumed to
-! be in the quad. This is done by bilinear interpolation, fitting
-! a function of the form a + bx + cy + dxy to the four points and 
-! then evaluating this function at (lon, lat). The fit is done by
-! solving the 4x4 system of equations for a, b, c, and d. The system
-! is reduced to a 3x3 by eliminating a from the first three equations
-! and then solving the 3x3 before back substituting. There is concern
-! about the numerical stability of this implementation. Implementation
-! checks showed accuracy to seven decimal places on all tests.
-
-integer :: i
-real(r8) :: m(3, 3), v(3), r(3), a, x_corners(4), lon
-! real(r8) :: lon_mean
-
-! Watch out for wraparound on x_corners.
-lon = lon_in
-x_corners = x_corners_in
-
-! See if the side wraps around in longitude. If the corners longitudes
-! wrap around 360, then the corners and the point to interpolate to
-! must be adjusted to be in the range from 180 to 540 degrees.
-if(maxval(x_corners) - minval(x_corners) > 180.0_r8) then
-   if(lon < 180.0_r8) lon = lon + 360.0_r8
-   do i = 1, 4
-      if(x_corners(i) < 180.0_r8) x_corners(i) = x_corners(i) + 360.0_r8
-   enddo
-endif
-
-
-!*******
-! Problems with extremes in polar cell interpolation can be reduced
-! by this block, but it is not clear that it is needed for actual
-! ocean grid data
-! Find the mean longitude of corners and remove
-!!!lon_mean = sum(x_corners) / 4.0_r8
-!!!x_corners = x_corners - lon_mean
-!!!lon = lon - lon_mean
-! Multiply everybody by the cos of the latitude
-!!!do i = 1, 4
-   !!!x_corners(i) = x_corners(i) * cos(y_corners(i) * deg2rad)
-!!!enddo
-!!!lon = lon * cos(lat * deg2rad)
-
-!*******
-
-
-! Fit a surface and interpolate; solve for 3x3 matrix
-do i = 1, 3
-   ! Eliminate a from the first 3 equations
-   m(i, 1) = x_corners(i) - x_corners(i + 1)
-   m(i, 2) = y_corners(i) - y_corners(i + 1)
-   m(i, 3) = x_corners(i)*y_corners(i) - x_corners(i + 1)*y_corners(i + 1)
-   v(i) = p(i) - p(i + 1)
-enddo
-
-! Solve the matrix for b, c and d
-call mat3x3(m, v, r)
-
-! r contains b, c, and d; solve for a
-a = p(4) - r(1) * x_corners(4) - r(2) * y_corners(4) - &
-   r(3) * x_corners(4)*y_corners(4)
-
-
-!----------------- Implementation test block
-! When interpolating on dipole x3 never exceeded 1e-9 error in this test
-!!!do i = 1, 4
-   !!!interp_val = a + r(1)*x_corners(i) + r(2)*y_corners(i)+ r(3)*x_corners(i)*y_corners(i)
-   !!!if(abs(interp_val - p(i)) > 1e-9) then
-      !!!write(*, *) 'large interp residual ', interp_val - p(i)
-   !!!endif
-!!!enddo
-!----------------- Implementation test block
-
-
-! Now do the interpolation
-interp_val = a + r(1)*lon + r(2)*lat + r(3)*lon*lat
-
-!********
-! Avoid exceeding maxima or minima as stopgap for poles problem
-! When doing bilinear interpolation in quadrangle, can get interpolated
-! values that are outside the range of the corner values
-if(interp_val > maxval(p)) then 
-   interp_val = maxval(p)
-else if(interp_val < minval(p)) then
-   interp_val = minval(p)
-endif
-!********
-
-end subroutine quad_bilinear_interp
-
-!------------------------------------------------------------
-
-subroutine mat3x3(m, v, r)
- real(r8),  intent(in) :: m(3, 3), v(3)
- real(r8), intent(out) :: r(3)
-
-! Solves rank 3 linear system mr = v for r
-! using Cramer's rule. This isn't the best choice
-! for speed or numerical stability so might want to replace
-! this at some point.
-
-real(r8) :: m_sub(3, 3), numer, denom
-integer  :: i
-
-! Compute the denominator, det(m)
-denom = deter3(m)
-
-! Loop to compute the numerator for each component of r
-do i = 1, 3
-   m_sub = m
-   m_sub(:, i) = v   
-   numer = deter3(m_sub)
-   r(i) = numer / denom
-enddo
-
-end subroutine mat3x3
-
-!------------------------------------------------------------
-
-function deter3(m)
- real(r8), intent(in) :: m(3, 3)
- real(r8)             :: deter3
-
-! Computes determinant of 3x3 matrix m
-
-deter3 = m(1,1)*m(2,2)*m(3,3) + m(1,2)*m(2,3)*m(3,1) + &
-         m(1,3)*m(2,1)*m(3,2) - m(3,1)*m(2,2)*m(1,3) - &
-         m(1,1)*m(2,3)*m(3,2) - m(3,3)*m(2,1)*m(1,2)
-
-end function deter3
-
-!------------------------------------------------------------
-
-subroutine height_bounds(lheight, nheights, hgt_array, bot, top, fract, istatus)
- real(r8),    intent(in) :: lheight
- integer,     intent(in) :: nheights
- real(r8),    intent(in) :: hgt_array(nheights)
- integer,    intent(out) :: bot, top
- real(r8),   intent(out) :: fract
- integer,    intent(out) :: istatus
-
-! Local variables
-integer   :: i
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! Succesful istatus is 0
-! Make any failure here return istatus in the 20s
-istatus = 0
-
-! The zc array contains the depths of the center of the vertical grid boxes
-! In this case (unlike how we handle the MIT depths), positive is really down.
-! FIXME: in the MIT model, we're given box widths and we compute the centers,
-! and we computed them with larger negative numbers being deeper.  Here,
-! larger positive numbers are deeper.
-
-! It is assumed that the top box is shallow and any observations shallower
-! than the depth of this box's center are just given the value of the
-! top box.
-if(lheight <= hgt_array(1)) then
-   top = 1
-   bot = 2
-   ! NOTE: the fract definition is the relative distance from bottom to top
-   fract = 1.0_r8 
-if (debug > 7) print *, 'above first level in height'
-if (debug > 7) print *, 'hgt_array, top, bot, fract=', hgt_array(1), top, bot, fract
-   return
-endif
-
-! Search through the boxes
-do i = 2, nheights
-   ! If the location is shallower than this entry, it must be in this box
-   if(lheight < hgt_array(i)) then
-      top = i -1
-      bot = i
-      fract = (hgt_array(bot) - lheight) / (hgt_array(bot) - hgt_array(top))
-if (debug > 7) print *, 'i, hgt_array, top, bot, fract=', i, hgt_array(i), top, bot, fract
-      return
-   endif
-enddo
-
-! Falling off the end means the location is lower than the deepest height
-istatus = 20
-
-end subroutine height_bounds
-
-!------------------------------------------------------------------
-
-function pop_get_model_time_step()
- type(time_type) :: pop_get_model_time_step
-
-! Returns the the time step of the model; the smallest increment
-! in time that the model is capable of advancing the state in a given
-! implementation. This interface is required for all applications.
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-pop_get_model_time_step = model_timestep
-
-end function pop_get_model_time_step
-
-!------------------------------------------------------------------
-
-subroutine pop_get_state_meta_data(index_in, location, var_type)
- integer,             intent(in)  :: index_in
- type(location_type), intent(out) :: location
- integer,             intent(out), optional :: var_type
-
-! Given an integer index into the state vector structure, returns the
-! associated location. A second intent(out) optional argument kind
-! can be returned if the model has more than one type of field (for
-! instance temperature and zonal wind component). This interface is
-! required for all filter applications as it is required for computing
-! the distance between observations and state variables.
-
-real(r8) :: lat, lon, depth
-integer :: lon_index, lat_index, depth_index, local_var
-
-call get_state_indices(index_in, lat_index, lon_index, depth_index, local_var)
-
-if (is_on_ugrid(local_var)) then
-   lon = ULON(lon_index, lat_index)
-   lat = ULAT(lon_index, lat_index)
-else
-   lon = TLON(lon_index, lat_index)
-   lat = TLAT(lon_index, lat_index)
-endif
-
-if (local_var == KIND_SEA_SURFACE_HEIGHT) then
-   depth = 0.0_r8
-else
-   depth = ZC(depth_index)
-endif
-
-if (debug > 5) print *, 'lon, lat, depth = ', lon, lat, depth
-
-location = set_location(lon, lat, depth, VERTISHEIGHT)
-
-if (present(var_type)) then
-   var_type = local_var
-   if(is_dry_land(var_type, lon_index, lat_index, depth_index)) then
-      var_type = KIND_DRY_LAND
-   endif
-endif
-
-end subroutine pop_get_state_meta_data
-
-!------------------------------------------------------------------
-
-subroutine get_state_indices(index_in, lat_index, lon_index, depth_index, var_type)
- integer, intent(in)  :: index_in
- integer, intent(out) :: lat_index, lon_index, depth_index
- integer, intent(out) :: var_type
-
-! Given an integer index into the state vector structure, returns the
-! associated array indices for lat, lon, and depth, as well as the type.
-
-integer :: startind, offset
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-if (debug > 5) print *, 'asking for meta data about index ', index_in
-
-call get_state_kind(index_in, var_type, startind, offset)
-
-if (startind == start_index(PSURF_index)) then
-  depth_index = 1
-else
-  depth_index = (offset / (Nx * Ny)) + 1
-endif
-
-lat_index = (offset - ((depth_index-1)*Nx*Ny)) / Nx + 1
-lon_index =  offset - ((depth_index-1)*Nx*Ny) - ((lat_index-1)*Nx) + 1
-
-if (debug > 5) print *, 'lon, lat, depth index = ', lon_index, lat_index, depth_index
-
-end subroutine get_state_indices
-
-!------------------------------------------------------------------
-
-subroutine get_state_kind(index_in, var_type, startind, offset)
- integer, intent(in)  :: index_in
- integer, intent(out) :: var_type, startind, offset
-
-! Given an integer index into the state vector structure, returns the kind,
-! and both the starting offset for this kind, as well as the offset into
-! the block of this kind.
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-if (debug > 5) print *, 'asking for meta data about index ', index_in
-
-if (index_in < start_index(S_index+1)) then
-   var_type = KIND_SALINITY  
-   startind = start_index(S_index)
-else if (index_in < start_index(T_index+1)) then
-   var_type = KIND_POTENTIAL_TEMPERATURE  
-   startind = start_index(T_index)
-else if (index_in < start_index(U_index+1)) then
-   var_type = KIND_U_CURRENT_COMPONENT
-   startind = start_index(U_index)
-else if (index_in < start_index(V_index+1)) then
-   var_type = KIND_V_CURRENT_COMPONENT
-   startind = start_index(V_index)
-else 
-   var_type = KIND_SEA_SURFACE_PRESSURE
-   startind = start_index(PSURF_index)
-endif
-
-! local offset into this var array
-offset = index_in - startind
-
-if (debug > 5) print *, 'var type = ', var_type
-if (debug > 5) print *, 'startind = ', startind
-if (debug > 5) print *, 'offset = ', offset
-
-end subroutine get_state_kind
-
-!------------------------------------------------------------------
-
-subroutine get_state_kind_inc_dry(index_in, var_type)
- integer, intent(in)  :: index_in
- integer, intent(out) :: var_type
-
-! Given an integer index into the state vector structure, returns the
-! type, taking into account the ocean bottom and dry land.
-
-integer :: lon_index, lat_index, depth_index, startind, offset
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-call get_state_kind(index_in, var_type, startind, offset)
-
-if (startind == start_index(PSURF_index)) then
-  depth_index = 1
-else
-  depth_index = (offset / (Nx * Ny)) + 1
-endif
-
-lat_index = (offset - ((depth_index-1)*Nx*Ny)) / Nx + 1
-lon_index =  offset - ((depth_index-1)*Nx*Ny) - ((lat_index-1)*Nx) + 1
-
-! if on land or below ocean floor, replace type with dry land.
-if(is_dry_land(var_type, lon_index, lat_index, depth_index)) then
-   var_type = KIND_DRY_LAND
-endif
-
-end subroutine get_state_kind_inc_dry
-
-!------------------------------------------------------------------
-
-subroutine pop_end_model()
-
-! Shutdown and clean-up.
-
-! assume if one is allocated, they all were.  if no one ever
-! called the init routine, don't try to dealloc something that
-! was never alloc'd.
-if (allocated(ULAT)) deallocate(ULAT, ULON, TLAT, TLON, KMT, KMU, HT, HU)
-if (allocated(ZC))   deallocate(ZC, ZG, pressure)
-
-end subroutine pop_end_model
-
-!------------------------------------------------------------------
-
-function pop_nc_write_model_atts( ncFileID ) result (ierr)
- integer, intent(in)  :: ncFileID      ! netCDF file identifier
- integer              :: ierr          ! return value of function
-
-! TJH -- Writes the model-specific attributes to a netCDF file.
-!     This includes coordinate variables and some metadata, but NOT
-!     the model state vector.
-!
-! assim_model_mod:init_diag_output uses information from the location_mod
-!     to define the location dimension and variable ID. All we need to do
-!     is query, verify, and fill ...
-!
-! Typical sequence for adding new dimensions,variables,attributes:
-! NF90_OPEN             ! open existing netCDF dataset
-!    NF90_redef         ! put into define mode 
-!    NF90_def_dim       ! define additional dimensions (if any)
-!    NF90_def_var       ! define variables: from name, type, and dims
-!    NF90_put_att       ! assign attribute values
-! NF90_ENDDEF           ! end definitions: leave define mode
-!    NF90_put_var       ! provide values for variable
-! NF90_CLOSE            ! close: save updated netCDF dataset
-
-integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
-
-!----------------------------------------------------------------------
-! variables if we just blast out one long state vector
-!----------------------------------------------------------------------
-
-integer :: StateVarDimID   ! netCDF pointer to state variable dimension (model size)
-integer :: MemberDimID     ! netCDF pointer to dimension of ensemble    (ens_size)
-integer :: TimeDimID       ! netCDF pointer to time dimension           (unlimited)
-
-integer :: StateVarVarID   ! netCDF pointer to state variable coordinate array
-integer :: StateVarID      ! netCDF pointer to 3D [state,copy,time] array
-
-!----------------------------------------------------------------------
-! variables if we parse the state vector into prognostic variables.
-!----------------------------------------------------------------------
-
-! for the dimensions and coordinate variables
-integer :: NlonDimID, NlatDimID, NzDimID
-integer :: ulonVarID, ulatVarID, tlonVarID, tlatVarID, ZGVarID, ZCVarID
-integer :: KMTVarID, KMUVarID
-
-! for the prognostic variables
-integer :: SVarID, TVarID, UVarID, VVarID, PSURFVarID 
-
-!----------------------------------------------------------------------
-! variables for the namelist output
-!----------------------------------------------------------------------
-
-character(len=129), allocatable, dimension(:) :: textblock
-integer :: LineLenDimID, nlinesDimID, nmlVarID
-integer :: nlines, linelen
-logical :: has_pop_namelist
-
-!----------------------------------------------------------------------
-! local variables 
-!----------------------------------------------------------------------
-
-! we are going to need these to record the creation date in the netCDF file.
-! This is entirely optional, but nice.
-
-character(len=8)      :: crdate      ! needed by F90 DATE_AND_TIME intrinsic
-character(len=10)     :: crtime      ! needed by F90 DATE_AND_TIME intrinsic
-character(len=5)      :: crzone      ! needed by F90 DATE_AND_TIME intrinsic
-integer, dimension(8) :: values      ! needed by F90 DATE_AND_TIME intrinsic
-character(len=NF90_MAX_NAME) :: str1
-
-integer :: i
-character(len=128)  :: filename
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-ierr = -1 ! assume things go poorly
-
-!--------------------------------------------------------------------
-! we only have a netcdf handle here so we do not know the filename
-! or the fortran unit number.  but construct a string with at least
-! the netcdf handle, so in case of error we can trace back to see
-! which netcdf file is involved.
-!--------------------------------------------------------------------
-
-write(filename,*) 'ncFileID', ncFileID
-
-!-------------------------------------------------------------------------------
-! make sure ncFileID refers to an open netCDF file, 
-! and then put into define mode.
-!-------------------------------------------------------------------------------
-
-call nc_check(nf90_Inquire(ncFileID,nDimensions,nVariables,nAttributes,unlimitedDimID),&
-                                   'nc_write_model_atts', 'inquire '//trim(filename))
-call nc_check(nf90_Redef(ncFileID),'nc_write_model_atts',   'redef '//trim(filename))
-
-!-------------------------------------------------------------------------------
-! We need the dimension ID for the number of copies/ensemble members, and
-! we might as well check to make sure that Time is the Unlimited dimension. 
-! Our job is create the 'model size' dimension.
-!-------------------------------------------------------------------------------
-
-call nc_check(nf90_inq_dimid(ncid=ncFileID, name='NMLlinelen', dimid=LineLenDimID), &
-                           'nc_write_model_atts','inq_dimid NMLlinelen')
-call nc_check(nf90_inq_dimid(ncid=ncFileID, name='copy', dimid=MemberDimID), &
-                           'nc_write_model_atts', 'copy dimid '//trim(filename))
-call nc_check(nf90_inq_dimid(ncid=ncFileID, name='time', dimid=  TimeDimID), &
-                           'nc_write_model_atts', 'time dimid '//trim(filename))
-
-if ( TimeDimID /= unlimitedDimId ) then
-   write(msgstring,*)'Time Dimension ID ',TimeDimID, &
-             ' should equal Unlimited Dimension ID',unlimitedDimID
-   call error_handler(E_ERR,'nc_write_model_atts', msgstring, source, revision, revdate)
-endif
-
-!-------------------------------------------------------------------------------
-! Define the model size / state variable dimension / whatever ...
-!-------------------------------------------------------------------------------
-call nc_check(nf90_def_dim(ncid=ncFileID, name='StateVariable', len=model_size, &
-        dimid = StateVarDimID),'nc_write_model_atts', 'state def_dim '//trim(filename))
-
-!-------------------------------------------------------------------------------
-! Write Global Attributes 
-!-------------------------------------------------------------------------------
-
-call DATE_AND_TIME(crdate,crtime,crzone,values)
-write(str1,'(''YYYY MM DD HH MM SS = '',i4,5(1x,i2.2))') &
-                  values(1), values(2), values(3), values(5), values(6), values(7)
-
-call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'creation_date' ,str1    ), &
-           'nc_write_model_atts', 'creation put '//trim(filename))
-call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model_source'  ,source  ), &
-           'nc_write_model_atts', 'source put '//trim(filename))
-call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model_revision',revision), &
-           'nc_write_model_atts', 'revision put '//trim(filename))
-call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model_revdate' ,revdate ), &
-           'nc_write_model_atts', 'revdate put '//trim(filename))
-call nc_check(nf90_put_att(ncFileID, NF90_GLOBAL, 'model',  'POP' ), &
-           'nc_write_model_atts', 'model put '//trim(filename))
-
-!-------------------------------------------------------------------------------
-! Determine shape of most important namelist
-!-------------------------------------------------------------------------------
-
-call find_textfile_dims('pop_in', nlines, linelen)
-if (nlines > 0) then
-  has_pop_namelist = .true.
-else
-  has_pop_namelist = .false.
-endif
-
-if (debug > 0)    print *, 'pop namelist: nlines, linelen = ', nlines, linelen
-  
-if (has_pop_namelist) then 
-   allocate(textblock(nlines))
-   textblock = ''
-
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='nlines', &
-                 len = nlines, dimid = nlinesDimID), &
-                 'nc_write_model_atts', 'def_dim nlines ')
-
-   call nc_check(nf90_def_var(ncFileID,name='pop_in', xtype=nf90_char,    &
-                 dimids = (/ linelenDimID, nlinesDimID /),  varid=nmlVarID), &
-                 'nc_write_model_atts', 'def_var pop_in')
-   call nc_check(nf90_put_att(ncFileID, nmlVarID, 'long_name',       &
-                 'contents of pop_in namelist'), 'nc_write_model_atts', 'put_att pop_in')
-
-endif
-
-!-------------------------------------------------------------------------------
-! Here is the extensible part. The simplest scenario is to output the state vector,
-! parsing the state vector into model-specific parts is complicated, and you need
-! to know the geometry, the output variables (PS,U,V,T,Q,...) etc. We're skipping
-! complicated part.
-!-------------------------------------------------------------------------------
-
-if ( output_state_vector ) then
-
-   !----------------------------------------------------------------------------
-   ! Create a variable for the state vector
-   !----------------------------------------------------------------------------
-
-  ! Define the state vector coordinate variable and some attributes.
-   call nc_check(nf90_def_var(ncid=ncFileID,name='StateVariable', xtype=nf90_int, &
-                 dimids=StateVarDimID, varid=StateVarVarID), 'nc_write_model_atts', &
-                 'statevariable def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,StateVarVarID,'long_name','State Variable ID'),&
-                 'nc_write_model_atts','statevariable long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, StateVarVarID, 'units','indexical'), &
-                 'nc_write_model_atts', 'statevariable units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,StateVarVarID,'valid_range',(/ 1,model_size /)),&
-                 'nc_write_model_atts', 'statevariable valid_range '//trim(filename))
-
-   ! Define the actual (3D) state vector, which gets filled as time goes on ... 
-   call nc_check(nf90_def_var(ncid=ncFileID, name='state', xtype=nf90_real, &
-                 dimids=(/StateVarDimID,MemberDimID,unlimitedDimID/),varid=StateVarID),&
-                 'nc_write_model_atts','state def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,StateVarID,'long_name','model state or fcopy'),&
-                 'nc_write_model_atts', 'state long_name '//trim(filename))
-
-   ! Leave define mode so we can fill the coordinate variable.
-   call nc_check(nf90_enddef(ncfileID),'nc_write_model_atts','state enddef '//trim(filename))
-
-   ! Fill the state variable coordinate variable
-   call nc_check(nf90_put_var(ncFileID, StateVarVarID, (/ (i,i=1,model_size) /) ), &
-                 'nc_write_model_atts', 'state put_var '//trim(filename))
-
-else
-
-   !----------------------------------------------------------------------------
-   ! We need to output the prognostic variables.
-   !----------------------------------------------------------------------------
-   ! Define the new dimensions IDs
-   !----------------------------------------------------------------------------
-   
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='i', &
-          len = Nx, dimid = NlonDimID),'nc_write_model_atts', 'i def_dim '//trim(filename))
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='j', &
-          len = Ny, dimid = NlatDimID),'nc_write_model_atts', 'j def_dim '//trim(filename))
-   call nc_check(nf90_def_dim(ncid=ncFileID, name='k', &
-          len = Nz, dimid =   NzDimID),'nc_write_model_atts', 'k def_dim '//trim(filename))
-   
-   !----------------------------------------------------------------------------
-   ! Create the (empty) Coordinate Variables and the Attributes
-   !----------------------------------------------------------------------------
-
-
-   ! U,V Grid Longitudes
-   call nc_check(nf90_def_var(ncFileID,name='ULON', xtype=nf90_real, &
-                 dimids=(/ NlonDimID, NlatDimID /), varid=ulonVarID),&
-                 'nc_write_model_atts', 'ULON def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulonVarID, 'long_name', 'longitudes of U,V grid'), &
-                 'nc_write_model_atts', 'ULON long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulonVarID, 'cartesian_axis', 'X'),  &
-                 'nc_write_model_atts', 'ULON cartesian_axis '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulonVarID, 'units', 'degrees_east'), &
-                 'nc_write_model_atts', 'ULON units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulonVarID, 'valid_range', (/ 0.0_r8, 360.0_r8 /)), &
-                 'nc_write_model_atts', 'ULON valid_range '//trim(filename))
-
-   ! U,V Grid Latitudes
-   call nc_check(nf90_def_var(ncFileID,name='ULAT', xtype=nf90_real, &
-                 dimids=(/ NlonDimID, NlatDimID /), varid=ulatVarID),&
-                 'nc_write_model_atts', 'ULAT def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulatVarID, 'long_name', 'latitudes of U,V grid'), &
-                 'nc_write_model_atts', 'ULAT long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulatVarID, 'cartesian_axis', 'Y'),   &
-                 'nc_write_model_atts', 'ULAT cartesian_axis '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulatVarID, 'units', 'degrees_north'),  &
-                 'nc_write_model_atts', 'ULAT units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID,  ulatVarID,'valid_range',(/ -90.0_r8, 90.0_r8 /)), &
-                 'nc_write_model_atts', 'ULAT valid_range '//trim(filename))
-
-   ! S,T,PSURF Grid Longitudes
-   call nc_check(nf90_def_var(ncFileID,name='TLON', xtype=nf90_real, &
-                 dimids=(/ NlonDimID, NlatDimID /), varid=tlonVarID),&
-                 'nc_write_model_atts', 'TLON def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlonVarID, 'long_name', 'longitudes of S,T,... grid'), &
-                 'nc_write_model_atts', 'TLON long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlonVarID, 'cartesian_axis', 'X'),   &
-                 'nc_write_model_atts', 'TLON cartesian_axis '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlonVarID, 'units', 'degrees_east'),  &
-                 'nc_write_model_atts', 'TLON units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlonVarID, 'valid_range', (/ 0.0_r8, 360.0_r8 /)), &
-                 'nc_write_model_atts', 'TLON valid_range '//trim(filename))
-
-
-   ! S,T,PSURF Grid (center) Latitudes
-   call nc_check(nf90_def_var(ncFileID,name='TLAT', xtype=nf90_real, &
-                 dimids= (/ NlonDimID, NlatDimID /), varid=tlatVarID), &
-                 'nc_write_model_atts', 'TLAT def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlatVarID, 'long_name', 'latitudes of S,T, ... grid'), &
-                 'nc_write_model_atts', 'TLAT long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlatVarID, 'cartesian_axis', 'Y'),   &
-                 'nc_write_model_atts', 'TLAT cartesian_axis '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlatVarID, 'units', 'degrees_north'),  &
-                 'nc_write_model_atts', 'TLAT units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, tlatVarID, 'valid_range', (/ -90.0_r8, 90.0_r8 /)), &
-                 'nc_write_model_atts', 'TLAT valid_range '//trim(filename))
-
-   ! Depths
-   call nc_check(nf90_def_var(ncFileID,name='ZG', xtype=nf90_real, &
-                 dimids=NzDimID, varid= ZGVarID), &
-                 'nc_write_model_atts', 'ZG def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZGVarID, 'long_name', 'depth at grid edges'), &
-                 'nc_write_model_atts', 'ZG long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZGVarID, 'cartesian_axis', 'Z'),   &
-                 'nc_write_model_atts', 'ZG cartesian_axis '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZGVarID, 'units', 'meters'),  &
-                 'nc_write_model_atts', 'ZG units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZGVarID, 'positive', 'down'),  &
-                 'nc_write_model_atts', 'ZG units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZGVarID, 'comment', &
-                  'more positive is closer to the center of the earth'),  &
-                 'nc_write_model_atts', 'ZG comment '//trim(filename))
-
-   ! Depths
-   call nc_check(nf90_def_var(ncFileID,name='ZC',xtype=nf90_real,dimids=NzDimID,varid=ZCVarID), &
-                 'nc_write_model_atts', 'ZC def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZCVarID, 'long_name', 'depth at grid centroids'), &
-                 'nc_write_model_atts', 'ZC long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZCVarID, 'cartesian_axis', 'Z'),   &
-                 'nc_write_model_atts', 'ZC cartesian_axis '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZCVarID, 'units', 'meters'),  &
-                 'nc_write_model_atts', 'ZC units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZCVarID, 'positive', 'down'),  &
-                 'nc_write_model_atts', 'ZC units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, ZCVarID, 'comment', &
-                  'more positive is closer to the center of the earth'),  &
-                 'nc_write_model_atts', 'ZC comment '//trim(filename))
-
-   ! Depth mask
-   call nc_check(nf90_def_var(ncFileID,name='KMT',xtype=nf90_int, &
-                 dimids= (/ NlonDimID, NlatDimID /), varid=KMTVarID), &
-                 'nc_write_model_atts', 'KMT def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMTVarID, 'long_name', 'lowest valid depth index at grid centroids'), &
-                 'nc_write_model_atts', 'KMT long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMTVarID, 'units', 'levels'),  &
-                 'nc_write_model_atts', 'KMT units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMTVarID, 'positive', 'down'),  &
-                 'nc_write_model_atts', 'KMT units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMTVarID, 'comment', &
-                  'more positive is closer to the center of the earth'),  &
-                 'nc_write_model_atts', 'KMT comment '//trim(filename))
-
-   ! Depth mask
-   call nc_check(nf90_def_var(ncFileID,name='KMU',xtype=nf90_int, &
-                 dimids= (/ NlonDimID, NlatDimID /), varid=KMUVarID), &
-                 'nc_write_model_atts', 'KMU def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMUVarID, 'long_name', 'lowest valid depth index at grid corners'), &
-                 'nc_write_model_atts', 'KMU long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMUVarID, 'units', 'levels'),  &
-                 'nc_write_model_atts', 'KMU units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMUVarID, 'positive', 'down'),  &
-                 'nc_write_model_atts', 'KMU units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, KMUVarID, 'comment', &
-                  'more positive is closer to the center of the earth'),  &
-                 'nc_write_model_atts', 'KMU comment '//trim(filename))
-
-   !----------------------------------------------------------------------------
-   ! Create the (empty) Prognostic Variables and the Attributes
-   !----------------------------------------------------------------------------
-
-
-   call nc_check(nf90_def_var(ncid=ncFileID, name='SALT', xtype=nf90_real, &
-         dimids = (/NlonDimID,NlatDimID,NzDimID,MemberDimID,unlimitedDimID/),varid=SVarID),&
-         'nc_write_model_atts', 'S def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, SVarID, 'long_name', 'salinity'), &
-         'nc_write_model_atts', 'S long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, SVarID, 'units', 'kg/kg'), &
-         'nc_write_model_atts', 'S units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, SVarID, 'missing_value', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'S missing '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, SVarID, '_FillValue', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'S fill '//trim(filename))
-
-
-   call nc_check(nf90_def_var(ncid=ncFileID, name='TEMP', xtype=nf90_real, &
-         dimids=(/NlonDimID,NlatDimID,NzDimID,MemberDimID,unlimitedDimID/),varid=TVarID),&
-         'nc_write_model_atts', 'T def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, TVarID, 'long_name', 'Potential Temperature'), &
-         'nc_write_model_atts', 'T long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, TVarID, 'units', 'deg C'), &
-         'nc_write_model_atts', 'T units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, TVarID, 'units_long_name', 'degrees celsius'), &
-         'nc_write_model_atts', 'T units_long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, TVarID, 'missing_value', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'T missing '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, TVarID, '_FillValue', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'T fill '//trim(filename))
-
-
-   call nc_check(nf90_def_var(ncid=ncFileID, name='UVEL', xtype=nf90_real, &
-         dimids=(/NlonDimID,NlatDimID,NzDimID,MemberDimID,unlimitedDimID/),varid=UVarID),&
-         'nc_write_model_atts', 'U def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, UVarID, 'long_name', 'U velocity'), &
-         'nc_write_model_atts', 'U long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, UVarID, 'units', 'cm/s'), &
-         'nc_write_model_atts', 'U units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, UVarID, 'units_long_name', 'centimeters per second'), &
-         'nc_write_model_atts', 'U units_long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, UVarID, 'missing_value', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'U missing '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, UVarID, '_FillValue', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'U fill '//trim(filename))
-
-
-   call nc_check(nf90_def_var(ncid=ncFileID, name='VVEL', xtype=nf90_real, &
-         dimids=(/NlonDimID,NlatDimID,NzDimID,MemberDimID,unlimitedDimID/),varid=VVarID),&
-         'nc_write_model_atts', 'V def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, VVarID, 'long_name', 'V Velocity'), &
-         'nc_write_model_atts', 'V long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, VVarID, 'units', 'cm/s'), &
-         'nc_write_model_atts', 'V units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, VVarID, 'units_long_name', 'centimeters per second'), &
-         'nc_write_model_atts', 'V units_long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, VVarID, 'missing_value', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'V missing '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, VVarID, '_FillValue', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'V fill '//trim(filename))
-
-
-   call nc_check(nf90_def_var(ncid=ncFileID, name='PSURF', xtype=nf90_real, &
-         dimids=(/NlonDimID,NlatDimID,MemberDimID,unlimitedDimID/),varid=PSURFVarID), &
-         'nc_write_model_atts', 'PSURF def_var '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, PSURFVarID, 'long_name', 'surface pressure'), &
-         'nc_write_model_atts', 'PSURF long_name '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, PSURFVarID, 'units', 'dyne/cm2'), &
-         'nc_write_model_atts', 'PSURF units '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, PSURFVarID, 'missing_value', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'PSURF missing '//trim(filename))
-   call nc_check(nf90_put_att(ncFileID, PSURFVarID, '_FillValue', NF90_FILL_REAL), &
-         'nc_write_model_atts', 'PSURF fill '//trim(filename))
-
-   ! Finished with dimension/variable definitions, must end 'define' mode to fill.
-
-   call nc_check(nf90_enddef(ncfileID), 'prognostic enddef '//trim(filename))
-
-   !----------------------------------------------------------------------------
-   ! Fill the coordinate variables
-   !----------------------------------------------------------------------------
-
-   call nc_check(nf90_put_var(ncFileID, ulonVarID, ULON ), &
-                'nc_write_model_atts', 'ULON put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, ulatVarID, ULAT ), &
-                'nc_write_model_atts', 'ULAT put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, tlonVarID, TLON ), &
-                'nc_write_model_atts', 'TLON put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, tlatVarID, TLAT ), &
-                'nc_write_model_atts', 'TLAT put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, ZGVarID, ZG ), &
-                'nc_write_model_atts', 'ZG put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, ZCVarID, ZC ), &
-                'nc_write_model_atts', 'ZC put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, KMTVarID, KMT ), &
-                'nc_write_model_atts', 'KMT put_var '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID, KMUVarID, KMU ), &
-                'nc_write_model_atts', 'KMU put_var '//trim(filename))
-
-endif
-
-!-------------------------------------------------------------------------------
-! Fill the variables we can
-!-------------------------------------------------------------------------------
-
-if (has_pop_namelist) then
-   call file_to_text('pop_in', textblock)
-   call nc_check(nf90_put_var(ncFileID, nmlVarID, textblock ), &
-                 'nc_write_model_atts', 'put_var nmlVarID')
-   deallocate(textblock)
-endif
-
-!-------------------------------------------------------------------------------
-! Flush the buffer and leave netCDF file open
-!-------------------------------------------------------------------------------
-call nc_check(nf90_sync(ncFileID), 'nc_write_model_atts', 'atts sync')
-
-ierr = 0 ! If we got here, things went well.
-
-end function pop_nc_write_model_atts
-
-!------------------------------------------------------------------
-
-function pop_nc_write_model_vars( ncFileID, statevec, copyindex, timeindex ) result (ierr)         
- integer,                intent(in) :: ncFileID      ! netCDF file identifier
- real(r8), dimension(:), intent(in) :: statevec
- integer,                intent(in) :: copyindex
- integer,                intent(in) :: timeindex
- integer                            :: ierr          ! return value of function
-
-! TJH 24 Oct 2006 -- Writes the model variables to a netCDF file.
-!
-! TJH 29 Jul 2003 -- for the moment, all errors are fatal, so the
-! return code is always '0 == normal', since the fatal errors stop execution.
-!
-! For the lorenz_96 model, each state variable is at a separate location.
-! that's all the model-specific attributes I can think of ...
-!
-! assim_model_mod:init_diag_output uses information from the location_mod
-!     to define the location dimension and variable ID. All we need to do
-!     is query, verify, and fill ...
-!
-! Typical sequence for adding new dimensions,variables,attributes:
-! NF90_OPEN             ! open existing netCDF dataset
-!    NF90_redef         ! put into define mode
-!    NF90_def_dim       ! define additional dimensions (if any)
-!    NF90_def_var       ! define variables: from name, type, and dims
-!    NF90_put_att       ! assign attribute values
-! NF90_ENDDEF           ! end definitions: leave define mode
-!    NF90_put_var       ! provide values for variable
-! NF90_CLOSE            ! close: save updated netCDF dataset
-
-integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
-integer :: VarID
-
-real(r8), dimension(Nx,Ny,Nz) :: data_3d
-real(r8), dimension(Nx,Ny)    :: data_2d
-character(len=128)  :: filename
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-ierr = -1 ! assume things go poorly
-
-!--------------------------------------------------------------------
-! we only have a netcdf handle here so we do not know the filename
-! or the fortran unit number.  but construct a string with at least
-! the netcdf handle, so in case of error we can trace back to see
-! which netcdf file is involved.
-!--------------------------------------------------------------------
-
-write(filename,*) 'ncFileID', ncFileID
-
-!-------------------------------------------------------------------------------
-! make sure ncFileID refers to an open netCDF file, 
-!-------------------------------------------------------------------------------
-
-call nc_check(nf90_Inquire(ncFileID,nDimensions,nVariables,nAttributes,unlimitedDimID),&
-              'nc_write_model_vars', 'inquire '//trim(filename))
-
-if ( output_state_vector ) then
-
-   call nc_check(NF90_inq_varid(ncFileID, 'state', VarID), &
-                 'nc_write_model_vars', 'state inq_varid '//trim(filename))
-   call nc_check(NF90_put_var(ncFileID,VarID,statevec,start=(/1,copyindex,timeindex/)),&
-                 'nc_write_model_vars', 'state put_var '//trim(filename))
-
-else
-
-   !----------------------------------------------------------------------------
-   ! We need to process the prognostic variables.
-   ! Replace missing values (0.0) with netcdf missing value.
-   ! Staggered grid causes some logistical problems.
-   !----------------------------------------------------------------------------
-
-   call vector_to_prog_var(statevec, S_index, data_3d)
-   where (data_3d == 0.0_r8) data_3d = NF90_FILL_REAL
-   call nc_check(NF90_inq_varid(ncFileID, 'SALT', VarID), &
-                'nc_write_model_vars', 'S inq_varid '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
-                'nc_write_model_vars', 'S put_var '//trim(filename))
-
-   call vector_to_prog_var(statevec, T_index, data_3d)
-   where (data_3d == 0.0_r8) data_3d = NF90_FILL_REAL
-   call nc_check(NF90_inq_varid(ncFileID, 'TEMP', VarID), &
-                'nc_write_model_vars', 'T inq_varid '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
-                'nc_write_model_vars', 'T put_var '//trim(filename))
-
-   call vector_to_prog_var(statevec, U_index, data_3d)
-   where (data_3d == 0.0_r8) data_3d = NF90_FILL_REAL
-   call nc_check(NF90_inq_varid(ncFileID, 'UVEL', VarID), &
-                'nc_write_model_vars', 'U inq_varid '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
-                'nc_write_model_vars', 'U put_var '//trim(filename))
-
-   call vector_to_prog_var(statevec, V_index, data_3d)
-   where (data_3d == 0.0_r8) data_3d = NF90_FILL_REAL
-   call nc_check(NF90_inq_varid(ncFileID, 'VVEL', VarID), &
-                'nc_write_model_vars', 'V inq_varid '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID,VarID,data_3d,start=(/1,1,1,copyindex,timeindex/)),&
-                'nc_write_model_vars', 'V put_var '//trim(filename))
-
-   call vector_to_prog_var(statevec, PSURF_index, data_2d)
-   where (data_2d == 0.0_r8) data_2d = NF90_FILL_REAL
-   call nc_check(NF90_inq_varid(ncFileID, 'PSURF', VarID), &
-                'nc_write_model_vars', 'PSURF inq_varid '//trim(filename))
-   call nc_check(nf90_put_var(ncFileID,VarID,data_2d,start=(/1,1,copyindex,timeindex/)),&
-                'nc_write_model_vars', 'PSURF put_var '//trim(filename))
-
-endif
-
-!-------------------------------------------------------------------------------
-! Flush the buffer and leave netCDF file open
-!-------------------------------------------------------------------------------
-
-call nc_check(nf90_sync(ncFileID), 'nc_write_model_vars', 'sync '//trim(filename))
-
-ierr = 0 ! If we got here, things went well.
-
-end function pop_nc_write_model_vars
-
-!------------------------------------------------------------------
-
-subroutine pop_pert_model_state(state, pert_state, interf_provided)
- real(r8), intent(in)  :: state(:)
- real(r8), intent(out) :: pert_state(:)
- logical,  intent(out) :: interf_provided
-
-! Perturbs a model state for generating initial ensembles.
-! The perturbed state is returned in pert_state.
-! A model may choose to provide a NULL INTERFACE by returning
-! .false. for the interf_provided argument. This indicates to
-! the filter that if it needs to generate perturbed states, it
-! may do so by adding a perturbation to each model state 
-! variable independently. The interf_provided argument
-! should be returned as .true. if the model wants to do its own
-! perturbing of states.
-
-integer :: i, var_type
-logical, save :: random_seq_init = .false.
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-interf_provided = .true.
-
-! Initialize my random number sequence
-if(.not. random_seq_init) then
-   call init_random_seq(random_seq, my_task_id())
-   random_seq_init = .true.
-endif
-
-! only perturb the actual ocean cells; leave the land and
-! ocean floor values alone.
-do i=1,size(state)
-   call get_state_kind_inc_dry(i, var_type)
-   if (var_type /= KIND_DRY_LAND) then
-      pert_state(i) = random_gaussian(random_seq, state(i), &
-                                      model_perturbation_amplitude)
-   else
-      pert_state(i) = state(i)
-   endif
-enddo
-
-
-end subroutine pop_pert_model_state
-
-!------------------------------------------------------------------
-
-subroutine pop_ens_mean_for_model(ens_mean)
- real(r8), intent(in) :: ens_mean(:)
-
-! If needed by the model interface, this is the current mean
-! for all state vector items across all ensembles. It is up to this
-! code to allocate space and save a copy if it is going to be used
-! later on.  For now, we are ignoring it.
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-
-end subroutine pop_ens_mean_for_model
-
-!------------------------------------------------------------------
-
-subroutine print_ranges(x)
- real(r8), intent(in) :: x(:)
-
-! intended for debugging use = print out the data min/max for each
-! field in the state vector, along with the starting and ending 
-! indices for each field.
-
-integer :: s, e
-
-!  S_index     = 1
-!  T_index     = 2
-!  U_index     = 3
-!  V_index     = 4
-!  PSURF_index = 5
-
-s = 1
-e = start_index(T_index)-1
-print *, 'min/max  salinity: ', minval(x(s:e)), maxval(x(s:e))
-print *, 's/e      salinity: ', s, e
-
-s = start_index(T_index)
-e = start_index(U_index)-1
-print *, 'min/max  pot temp: ', minval(x(s:e)), maxval(x(s:e))
-print *, 's/e      pot temp: ', s, e
-
-s = start_index(U_index)
-e = start_index(V_index)-1
-print *, 'min/max U current: ', minval(x(s:e)), maxval(x(s:e))
-print *, 's/e     U current: ', s, e
-
-s = start_index(V_index)
-e = start_index(PSURF_index)-1
-print *, 'min/max V current: ', minval(x(s:e)), maxval(x(s:e))
-print *, 's/e     V current: ', s, e
-
-s = start_index(PSURF_index)
-e = size(x)
-print *, 'min/max Surf Pres: ', minval(x(s:e)), maxval(x(s:e))
-print *, 's/e     Surf Pres: ', s, e
-
-end subroutine print_ranges
-
-!------------------------------------------------------------------
-
-subroutine pop_restart_file_to_sv(filename, state_vector, model_time)
- character(len=*), intent(in)    :: filename 
- real(r8),         intent(inout) :: state_vector(:)
- type(time_type),  intent(out)   :: model_time
-
-! Reads the current time and state variables from a POP restart
-! file and packs them into a dart state vector.
-
-! temp space to hold data while we are reading it
-real(r8) :: data_2d_array(Nx,Ny), data_3d_array(Nx,Ny,Nz)
-integer  :: i, j, k, ivar, indx
-
-integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
-character(len=NF90_MAX_NAME) :: varname 
-integer :: VarID, numdims, dimlen
-integer :: ncid, iyear, imonth, iday, ihour, iminute, isecond
-character(len=256) :: myerrorstring 
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-state_vector = MISSING_R8
-
-! Check that the input file exists ... 
-! Read the time data. 
-! Note from Nancy Norton as pertains time:
-! "The time recorded in the pop2 restart files is the current time,
-! which corresponds to the time of the XXXX_CUR variables.
-!
-! current time is determined from iyear, imonth, iday, and *seconds_this_day*
-!
-! The ihour, iminute, and isecond variables are used for internal
-! model counting purposes, but because isecond is rounded to the nearest
-! integer, it is possible that using ihour,iminute,isecond information
-! on the restart file to determine the exact curtime would give you a 
-! slightly wrong answer."
-!
-! DART only knows about integer number of seconds, so using the rounded one
-! is what we would have to do anyway ... and we already have a set_date routine
-! that takes ihour, iminute, isecond information.
-
-if ( .not. file_exist(filename) ) then
-   write(msgstring,*) 'cannot open file ', trim(filename),' for reading.'
-   call error_handler(E_ERR,'restart_file_to_sv',msgstring,source,revision,revdate)
-endif
-
-call nc_check( nf90_open(trim(filename), NF90_NOWRITE, ncid), &
-                  'restart_file_to_sv', 'open '//trim(filename))
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'iyear'  , iyear), &
-                  'restart_file_to_sv', 'get_att iyear')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'imonth' , imonth), &
-                  'restart_file_to_sv', 'get_att imonth')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'iday'   , iday), &
-                  'restart_file_to_sv', 'get_att iday')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'ihour'  , ihour), &
-                  'restart_file_to_sv', 'get_att ihour')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'iminute', iminute), &
-                  'restart_file_to_sv', 'get_att iminute')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'isecond', isecond), &
-                  'restart_file_to_sv', 'get_att isecond')
-
-! FIXME: we don't allow a real year of 0 - add one for now, but
-! THIS MUST BE FIXED IN ANOTHER WAY!
-if (iyear == 0) then
-  call error_handler(E_MSG, 'restart_file_to_sv', &
-                     'WARNING!!!   year 0 not supported; setting to year 1')
-  iyear = 1
-endif
-
-model_time = set_date(iyear, imonth, iday, ihour, iminute, isecond)
-
-if (do_output()) &
-    call print_time(model_time,'time for restart file '//trim(filename))
-if (do_output()) &
-    call print_date(model_time,'date for restart file '//trim(filename))
-
-! Start counting and filling the state vector one item at a time,
-! repacking the 3d arrays into a single 1d list of numbers.
-! These must be a fixed number and in a fixed order.
-
-indx = 1
-
-! fill SALT, TEMP, UVEL, VVEL in that order
-! The POP restart files have two time steps for each variable,
-! the variables are named SALT_CUR and SALT_OLD ... for example.
-! We are only interested in the CURrent time step.
-
-do ivar=1, n3dfields
-
-   varname = trim(progvarnames(ivar))//'_CUR'
-   myerrorstring = trim(filename)//' '//trim(varname)
-
-   ! Is the netCDF variable the right shape?
-
-   call nc_check(nf90_inq_varid(ncid,   varname, VarID), &
-            'restart_file_to_sv', 'inq_varid '//trim(myerrorstring))
-
-   call nc_check(nf90_inquire_variable(ncid,VarId,dimids=dimIDs,ndims=numdims), &
-            'restart_file_to_sv', 'inquire '//trim(myerrorstring))
-
-   if (numdims /= 3) then
-      write(msgstring,*) trim(myerrorstring),' does not have exactly 3 dimensions'
-      call error_handler(E_ERR,'restart_file_to_sv',msgstring,source,revision,revdate)
-   endif
-
-   do i = 1,numdims
-      write(msgstring,'(''inquire dimension'',i2,A)') i,trim(myerrorstring)
-      call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
-            'restart_file_to_sv', msgstring)
-
-      if (dimlen /= size(data_3d_array,i)) then
-         write(msgstring,*) trim(myerrorstring),'dim/dimlen',i,dimlen,'not',size(data_3d_array,i)
-         call error_handler(E_ERR,'restart_file_to_sv',msgstring,source,revision,revdate)
-      endif
-   enddo   
-
-   ! Actually get the variable and stuff it into the array
-
-   call nc_check(nf90_get_var(ncid, VarID, data_3d_array), 'restart_file_to_sv', &
-                'get_var '//trim(varname))
-
-   do k = 1, Nz   ! size(data_3d_array,3)
-   do j = 1, Ny   ! size(data_3d_array,2)
-   do i = 1, Nx   ! size(data_3d_array,1)
-      state_vector(indx) = data_3d_array(i, j, k)
-      indx = indx + 1
-   enddo
-   enddo
-   enddo
-
-enddo
-
-! and finally, PSURF (and any other 2d fields)
-do ivar=(n3dfields+1), (n3dfields+n2dfields)
-
-   varname = trim(progvarnames(ivar))//'_CUR'
-   myerrorstring = trim(varname)//' '//trim(filename)
-
-   ! Is the netCDF variable the right shape?
-
-   call nc_check(nf90_inq_varid(ncid,   varname, VarID), &
-            'restart_file_to_sv', 'inq_varid '//trim(myerrorstring))
-
-   call nc_check(nf90_inquire_variable(ncid,VarId,dimids=dimIDs,ndims=numdims), &
-            'restart_file_to_sv', 'inquire '//trim(myerrorstring))
-
-   if (numdims /= 2) then
-      write(msgstring,*) trim(myerrorstring),' does not have exactly 2 dimensions'
-      call error_handler(E_ERR,'restart_file_to_sv',msgstring,source,revision,revdate)
-   endif
-
-   do i = 1,numdims
-      write(msgstring,'(''inquire dimension'',i2,A)') i,trim(myerrorstring)
-      call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
-            'restart_file_to_sv', msgstring)
-
-      if (dimlen /= size(data_2d_array,i)) then
-         write(msgstring,*) trim(myerrorstring),'dim/dimlen',i,dimlen,'not',size(data_2d_array,i)
-         call error_handler(E_ERR,'restart_file_to_sv',msgstring,source,revision,revdate)
-      endif
-   enddo   
-
-   ! Actually get the variable and stuff it into the array
-
-   call nc_check(nf90_get_var(ncid, VarID, data_2d_array), 'restart_file_to_sv', &
-                'get_var '//trim(varname))
-
-   do j = 1, Ny   ! size(data_3d_array,2)
-   do i = 1, Nx   ! size(data_3d_array,1)
-      state_vector(indx) = data_2d_array(i, j)
-      indx = indx + 1
-   enddo
-   enddo
-
-enddo
-
-end subroutine pop_restart_file_to_sv
-
-!------------------------------------------------------------------
-
-subroutine pop_sv_to_restart_file(state_vector, filename, statedate)
- real(r8),         intent(in) :: state_vector(:)
- character(len=*), intent(in) :: filename 
- type(time_type),  intent(in) :: statedate
-
-! Writes the current time and state variables from a dart state
-! vector (1d fortran array) into a POP netcdf restart file.
-
-integer :: iyear, imonth, iday, ihour, iminute, isecond
-type(time_type) :: pop_time
-
-! temp space to hold data while we are writing it
-real(r8) :: data_2d_array(Nx,Ny), data_3d_array(Nx,Ny,Nz)
-
-integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
-character(len=NF90_MAX_NAME)          :: varname 
-character(len=256)                    :: myerrorstring 
-
-integer :: i, ivar, ncid, VarID, numdims, dimlen
-
-!----------------------------------------------------------------------
-! Get the show underway
-!----------------------------------------------------------------------
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-! Check that the input file exists. 
-! make sure the time tag in the restart file matches 
-! the current time of the DART state ...
-
-if ( .not. file_exist(filename)) then
-   write(msgstring,*)trim(filename),' does not exist. FATAL error.'
-   call error_handler(E_ERR,'sv_to_restart_file',msgstring,source,revision,revdate) 
-endif
-
-call nc_check( nf90_open(trim(filename), NF90_WRITE, ncid), &
-                  'sv_to_restart_file', 'open '//trim(filename))
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'iyear'  , iyear), &
-                  'sv_to_restart_file', 'get_att iyear')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'imonth' , imonth), &
-                  'sv_to_restart_file', 'get_att imonth')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'iday'   , iday), &
-                  'sv_to_restart_file', 'get_att iday')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'ihour'  , ihour), &
-                  'sv_to_restart_file', 'get_att ihour')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'iminute', iminute), &
-                  'sv_to_restart_file', 'get_att iminute')
-call nc_check( nf90_get_att(ncid, NF90_GLOBAL, 'isecond', isecond), &
-                  'sv_to_restart_file', 'get_att isecond')
-
-pop_time = set_date(iyear, imonth, iday, ihour, iminute, isecond)
-
-if ( pop_time /= statedate ) then
-   call print_time(statedate,'DART current time',logfileunit) 
-   call print_time( pop_time,'POP  current time',logfileunit) 
-   call print_time(statedate,'DART current time') 
-   call print_time( pop_time,'POP  current time') 
-   write(msgstring,*)trim(filename),' current time /= model time. FATAL error.'
-   call error_handler(E_ERR,'sv_to_restart_file',msgstring,source,revision,revdate) 
-endif
-
-if (do_output()) &
-    call print_time(pop_time,'time of restart file '//trim(filename))
-if (do_output()) &
-    call print_date(pop_time,'date of restart file '//trim(filename))
-
-! fill S, T, U, V in that order
-do ivar=1, n3dfields
-
-   varname = trim(progvarnames(ivar))//'_CUR'
-   myerrorstring = trim(filename)//' '//trim(varname)
-
-   ! Is the netCDF variable the right shape?
-   call nc_check(nf90_inq_varid(ncid,   varname, VarID), &
-            'sv_to_restart_file', 'inq_varid '//trim(myerrorstring))
-
-   call nc_check(nf90_inquire_variable(ncid,VarId,dimids=dimIDs,ndims=numdims), &
-            'sv_to_restart_file', 'inquire '//trim(myerrorstring))
-
-   if (numdims /= 3) then
-      write(msgstring,*) trim(myerrorstring),' does not have exactly 3 dimensions'
-      call error_handler(E_ERR,'sv_to_restart_file',msgstring,source,revision,revdate)
-   endif
-
-   do i = 1,numdims
-      write(msgstring,'(''inquire dimension'',i2,A)') i,trim(myerrorstring)
-      call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
-            'sv_to_restart_file', msgstring)
-
-      if (dimlen /= size(data_3d_array,i)) then
-         write(msgstring,*) trim(myerrorstring),'dim/dimlen',i,dimlen,'not',size(data_3d_array,i)
-         call error_handler(E_ERR,'sv_to_restart_file',msgstring,source,revision,revdate)
-      endif
-   enddo
-
-   call vector_to_prog_var(state_vector, ivar, data_3d_array)
-
-   ! Actually stuff it into the netcdf file
-   call nc_check(nf90_put_var(ncid, VarID, data_3d_array), &
-            'sv_to_restart_file', 'put_var '//trim(myerrorstring))
-
-enddo
-
-! and finally, PSURF (and any other 2d fields)
-do ivar=(n3dfields+1), (n3dfields+n2dfields)
-
-   varname = trim(progvarnames(ivar))//'_CUR'
-   myerrorstring = trim(varname)//' '//trim(filename)
-
-   ! Is the netCDF variable the right shape?
-
-   call nc_check(nf90_inq_varid(ncid,   varname, VarID), &
-            'sv_to_restart_file', 'inq_varid '//trim(myerrorstring))
-
-   call nc_check(nf90_inquire_variable(ncid,VarId,dimids=dimIDs,ndims=numdims), &
-            'sv_to_restart_file', 'inquire '//trim(myerrorstring))
-
-   if (numdims /= 2) then
-      write(msgstring,*) trim(myerrorstring),' does not have exactly 2 dimensions'
-      call error_handler(E_ERR,'sv_to_restart_file',msgstring,source,revision,revdate)
-   endif
-
-   do i = 1,numdims
-      write(msgstring,'(''inquire dimension'',i2,A)') i,trim(myerrorstring)
-      call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
-            'sv_to_restart_file', msgstring)
-
-      if (dimlen /= size(data_2d_array,i)) then
-         write(msgstring,*) trim(myerrorstring),'dim/dimlen',i,dimlen,'not',size(data_2d_array,i)
-         call error_handler(E_ERR,'sv_to_restart_file',msgstring,source,revision,revdate)
-      endif
-   enddo
-
-   call vector_to_prog_var(state_vector, ivar, data_2d_array)
-
-   call nc_check(nf90_put_var(ncid, VarID, data_2d_array), &
-            'sv_to_restart_file', 'put_var '//trim(myerrorstring))
-
-enddo
-
-call nc_check(nf90_close(ncid), 'sv_to_restart_file', 'close '//trim(filename))
-
-end subroutine pop_sv_to_restart_file
-
-!------------------------------------------------------------------
-
-subroutine vector_to_2d_prog_var(x, varindex, data_2d_array)
- real(r8), dimension(:),   intent(in)  :: x
- integer,                  intent(in)  :: varindex
- real(r8), dimension(:,:), intent(out) :: data_2d_array
-
-! convert the values from a 1d fortran array, starting at an offset,
-! into a 2d fortran array.  the 2 dims are taken from the array size.
-
-integer :: i,j,ii
-integer :: dim1,dim2
-character(len=128) :: varname
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-dim1 = size(data_2d_array,1)
-dim2 = size(data_2d_array,2)
-
-varname = progvarnames(varindex)
-
-if (dim1 /= Nx) then
-   write(msgstring,*)trim(varname),' 2d array dim 1 ',dim1,' /= ',Nx
-   call error_handler(E_ERR,'vector_to_2d_prog_var',msgstring,source,revision,revdate) 
-endif
-if (dim2 /= Ny) then
-   write(msgstring,*)trim(varname),' 2d array dim 2 ',dim2,' /= ',Ny
-   call error_handler(E_ERR,'vector_to_2d_prog_var',msgstring,source,revision,revdate) 
-endif
-
-ii = start_index(varindex)
-
-do j = 1,Ny   ! latitudes
-do i = 1,Nx   ! longitudes
-   data_2d_array(i,j) = x(ii)
-   ii = ii + 1
-enddo
-enddo
-
-end subroutine vector_to_2d_prog_var
-
-!------------------------------------------------------------------
-
-subroutine vector_to_3d_prog_var(x, varindex, data_3d_array)
- real(r8), dimension(:),     intent(in)  :: x
- integer,                    intent(in)  :: varindex
- real(r8), dimension(:,:,:), intent(out) :: data_3d_array
-
-! convert the values from a 1d fortran array, starting at an offset,
-! into a 3d fortran array.  the 3 dims are taken from the array size.
-
-integer :: i,j,k,ii
-integer :: dim1,dim2,dim3
-character(len=128) :: varname
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-dim1 = size(data_3d_array,1)
-dim2 = size(data_3d_array,2)
-dim3 = size(data_3d_array,3)
-
-varname = progvarnames(varindex)
-
-if (dim1 /= Nx) then
-   write(msgstring,*)trim(varname),' 3d array dim 1 ',dim1,' /= ',Nx
-   call error_handler(E_ERR,'vector_to_3d_prog_var',msgstring,source,revision,revdate) 
-endif
-if (dim2 /= Ny) then
-   write(msgstring,*)trim(varname),' 3d array dim 2 ',dim2,' /= ',Ny
-   call error_handler(E_ERR,'vector_to_3d_prog_var',msgstring,source,revision,revdate) 
-endif
-if (dim3 /= Nz) then
-   write(msgstring,*)trim(varname),' 3d array dim 3 ',dim3,' /= ',Nz
-   call error_handler(E_ERR,'vector_to_3d_prog_var',msgstring,source,revision,revdate) 
-endif
-
-ii = start_index(varindex)
-
-do k = 1,Nz   ! vertical
-do j = 1,Ny   ! latitudes
-do i = 1,Nx   ! longitudes
-   data_3d_array(i,j,k) = x(ii)
-   ii = ii + 1
-enddo
-enddo
-enddo
-
-end subroutine vector_to_3d_prog_var
-
-!------------------------------------------------------------------
-
-subroutine pop_get_gridsize(num_x, num_y, num_z)
- integer, intent(out) :: num_x, num_y, num_z
-
-! public utility routine.
-
-if ( .not. module_initialized ) call pop_static_init_model
-
- num_x = Nx
- num_y = Ny
- num_z = Nz
-
-end subroutine pop_get_gridsize
-
-!------------------------------------------------------------------
-
-function is_dry_land(obs_type, lon_index, lat_index, hgt_index)
- integer, intent(in)  :: obs_type
- integer, intent(in)  :: lon_index, lat_index, hgt_index
- logical              :: is_dry_land
-
-! returns true if this point is below the ocean floor or if it is
-! on land.
-
-logical :: is_ugrid
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-is_dry_land = .FALSE.    ! start out thinking everything is wet.
-
-is_ugrid = is_on_ugrid(obs_type)
-if ((      is_ugrid .and. hgt_index > KMU(lon_index, lat_index)) .or. &
-    (.not. is_ugrid .and. hgt_index > KMT(lon_index, lat_index))) then
-   is_dry_land = .TRUE.
-   return
-endif
-
-end function is_dry_land
-
-!------------------------------------------------------------------
-
-function is_on_ugrid(obs_type)
- integer, intent(in) :: obs_type
- logical             :: is_on_ugrid
-
-!  returns true if U, V -- everything else is on T grid
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-is_on_ugrid = .FALSE.
-
-if ((obs_type == KIND_U_CURRENT_COMPONENT)  .or.  &
-    (obs_type == KIND_V_CURRENT_COMPONENT)) is_on_ugrid = .TRUE.
-
-end function is_on_ugrid
-
-!------------------------------------------------------------------
-
-subroutine write_grid_netcdf()
-
-! Write the grid to a netcdf file for checking.
-
-integer :: ncid, NlonDimID, NlatDimID, NzDimID
-integer :: nlon, nlat, nz
-integer :: ulatVarID, ulonVarID, TLATvarid, TLONvarid
-integer :: ZGvarid, ZCvarid, KMTvarid, KMUvarid
-
-integer :: dimids(2);
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-nlon = size(ULAT,1)
-nlat = size(ULAT,2)
-nz   = size(ZG)
-
-call nc_check(nf90_create('dart_grid.nc', NF90_CLOBBER, ncid),'write_grid_netcdf')
-
-! define dimensions
-
-call nc_check(nf90_def_dim(ncid, 'i', nlon, NlonDimID),'write_grid_netcdf')
-call nc_check(nf90_def_dim(ncid, 'j', nlat, NlatDimID),'write_grid_netcdf')
-call nc_check(nf90_def_dim(ncid, 'k',   nz,   NzDimID),'write_grid_netcdf')
-
-dimids(1) = NlonDimID 
-dimids(2) = NlatDimID 
-
-! define variables
-
-! FIXME: we should add attributes to say what units the grids are in (degrees).
-call nc_check(nf90_def_var(ncid,  'KMT', nf90_int,     dimids,  KMTvarid),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid,  'KMU', nf90_int,     dimids,  KMUvarid),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid, 'ULON', nf90_double,  dimids, ulonVarID),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid, 'ULAT', nf90_double,  dimids, ulatVarID),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid, 'TLON', nf90_double,  dimids, TLONvarid),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid, 'TLAT', nf90_double,  dimids, TLATvarid),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid,   'ZG', nf90_double, NzDimID,   ZGvarid),'write_grid_netcdf')
-call nc_check(nf90_def_var(ncid,   'ZC', nf90_double, NzDimID,   ZCvarid),'write_grid_netcdf')
-
-call nc_check(nf90_put_att(ncid,ulonVarID,'long_name','U,V grid lons'), &
-                                                     'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,ulatVarID,'long_name','U,V grid lats'), &
-                                                     'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,tlonVarID,'long_name','S,T grid lons'), &
-                                                     'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,tlatVarID,'long_name','S,T grid lats'), &
-                                                    'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,ZCVarID,'long_name','vertical grid centers'), &
-                                                    'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,ZCVarID,'units','meters'), &
-                                                    'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,ZGVarID,'long_name','vertical grid bottoms'), &
-                                                    'write_grid_netcdf')
-call nc_check(nf90_put_att(ncid,ZGVarID,'units','meters'), &
-                                                    'write_grid_netcdf')
-
-call nc_check(nf90_enddef(ncid),'write_grid_netcdf')
-
-! fill variables
-
-call nc_check(nf90_put_var(ncid,  KMTvarid,  KMT),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid,  KMUvarid,  KMU),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid, ulatVarID, ULAT),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid, ulonVarID, ULON),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid, TLATvarid, TLAT),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid, TLONvarid, TLON),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid,   ZGvarid,   ZG),'write_grid_netcdf')
-call nc_check(nf90_put_var(ncid,   ZCvarid,   ZC),'write_grid_netcdf')
-
-call nc_check(nf90_close(ncid),'write_grid_netcdf')
-
-end subroutine write_grid_netcdf
-
-!------------------------------------------------------------------
-
-subroutine pop_get_close_obs(gc, base_obs_loc, base_obs_kind, &
-                         obs, obs_kind, num_close, close_ind, dist)
-
- type(get_close_type),              intent(in) :: gc
- type(location_type),               intent(in) :: base_obs_loc
- integer,                           intent(in) :: base_obs_kind
- type(location_type), dimension(:), intent(in) :: obs
- integer,             dimension(:), intent(in) :: obs_kind
- integer,                           intent(out):: num_close
- integer,             dimension(:), intent(out):: close_ind
- real(r8),  optional, dimension(:), intent(out):: dist
-
-! Given a DART location (referred to as "base") and a set of candidate
-! locations & kinds (obs, obs_kind), returns the subset close to the
-! "base", their indices, and their distances to the "base" ...
-
-! For vertical distance computations, general philosophy is to convert all
-! vertical coordinates to a common coordinate. This coordinate type is defined
-! in the namelist with the variable "vert_localization_coord".
-
-integer :: t_ind, k
-
-! Initialize variables to missing status
-
-num_close = 0
-close_ind = -99
-if (present(dist)) dist = 1.0e9   !something big and positive (far away)
-
-! Get all the potentially close obs but no dist (optional argument dist(:)
-! is not present) This way, we are decreasing the number of distance
-! computations that will follow.  This is a horizontal-distance operation and
-! we don't need to have the relevant vertical coordinate information yet 
-! (for obs).
-
-call loc_get_close_obs(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
-                       num_close, close_ind)
-
-! Loop over potentially close subset of obs priors or state variables
-if (present(dist)) then
-do k = 1, num_close
-
-   t_ind = close_ind(k)
-
-   ! if dry land, leave original 1e9 value.  otherwise, compute real dist.
-   if (obs_kind(t_ind) /= KIND_DRY_LAND) then
-      dist(k) = get_dist(base_obs_loc,       obs(t_ind), &
-                         base_obs_kind, obs_kind(t_ind))
-   endif
-
-enddo
-endif
-
-end subroutine pop_get_close_obs
-
-!------------------------------------------------------------------
-
-subroutine write_grid_interptest()
-
-! Write the grid to an ascii file - in a format suitable for
-! subsequent use in the 'test_interpolation()' code.
-! write_grid_interptest is only possible after reading a real POP grid,
-! so static_init_model() must be called to gather the real POP grid.
-
-integer  :: i, j
-real(r8) :: rowmat(Nx,1), colmat(1,Ny), dmat(Nx,Ny)
-real(r8) :: rtlon, rulon, rtlat, rulat, u_val, t_val
-
-if ( .not. module_initialized ) call pop_static_init_model
-
-!----------------------------------------------------------------------
-! Generate a 'Regular' grid with the same rough 'shape' as the dipole grid
-!----------------------------------------------------------------------
-
-open(unit=12, position='rewind', action='write', file='regular_grid_u')
-open(unit=13, position='rewind', action='write', file='regular_grid_t')
-open(unit=14, position='rewind', action='write', file='regular_grid_u_data')
-open(unit=15, position='rewind', action='write', file='regular_grid_t_data')
-
-write(12, *) Nx, Ny
-write(13, *) Nx, Ny
-
-! Have T-grid starting at 0 and U grid offset by half
-do i = 1, Nx
-   rtlon = (i - 1.0_r8) * 360.0_r8 / Nx
-   rulon = (i - 0.5_r8) * 360.0_r8 / Nx
-   do j = 1, Ny
-      rtlat = -90.0_r8 + (j - 1.0_r8) * 180.0_r8 / Ny
-      rulat = -90.0_r8 + (j - 0.5_r8) * 180.0_r8 / Ny
-      write(12, *) i, j, rulon, rulat
-      write(13, *) i, j, rtlon, rtlat
-      ! Now add some wave pattern data 
-      u_val = sin(3.0_r8*(rulat + 11.0_r8)*2.0_r8*PI/360.0_r8) * &
-              sin(4.0_r8*(rulon + 17.0_r8)*2.0_r8*PI/360.0_r8)
-      t_val = sin(3.0_r8*(rtlat + 11.0_r8)*2.0_r8*PI/360.0_r8) * &
-              sin(4.0_r8*(rtlon + 17.0_r8)*2.0_r8*PI/360.0_r8)
-      write(14, *) rulon, rulat, u_val
-      write(15, *) rtlon, rtlat, t_val
-   enddo
-enddo
-
-close(unit=12)
-close(unit=13)
-close(unit=14)
-close(unit=15)
-
-!----------------------------------------------------------------------
-! POP grid (dipole) next
-!----------------------------------------------------------------------
-
-open(unit=12, position='rewind', action='write', file='dipole_grid_u')
-open(unit=13, position='rewind', action='write', file='dipole_grid_t')
-open(unit=14, position='rewind', action='write', file='dipole_grid_u_data')
-open(unit=15, position='rewind', action='write', file='dipole_grid_t_data')
-
-write(12, *) Nx, Ny
-write(13, *) Nx, Ny
-
-rowmat(:,1) = cos(PI * real((/ (i,i=0,Nx-1) /),r8) / Nx);
-colmat(1,:) = sin(PI * real((/ (i,i=0,Ny-1) /),r8) / Ny);
-dmat        = matmul(rowmat,colmat)
-
-do i = 1, Nx
-   do j = 1, Ny
-      write(12, *) i, j, ULON(i,j), ULAT(i,j)
-      write(13, *) i, j, TLON(i,j), TLAT(i,j)
-      write(14, *)       ULON(i,j), ULAT(i,j), dmat(i, j)
-      write(15, *)       TLON(i,j), TLAT(i,j), dmat(i, j)
-   enddo
-enddo
-
-close(unit=12)
-close(unit=13)
-close(unit=14)
-close(unit=15)
-
-end subroutine write_grid_interptest
-
-!------------------------------------------------------------------
-
-subroutine pop_test_interpolation(test_casenum)
- integer, intent(in) :: test_casenum
-
-! rigorous test of the interpolation routine.
-
-! This is storage just used for temporary test driver
-integer :: imain, jmain, index, istatus, nx_temp, ny_temp
-integer :: dnx_temp, dny_temp, height
-real(r8) :: ti, tj
-
-! Storage for testing interpolation to a different grid
-integer :: dnx, dny
-real(r8), allocatable :: dulon(:, :), dulat(:, :), dtlon(:, :), dtlat(:, :)
-real(r8), allocatable :: reg_u_data(:, :), reg_t_data(:, :)
-real(r8), allocatable :: reg_u_x(:), reg_t_x(:), dipole_u(:, :), dipole_t(:, :)
-
-! Test program reads in two grids; one is the interpolation grid
-! The second is a grid to which to interpolate.
-
-! Read of first grid
-! Set of block options for now
-! Lon lat for u on channel 12, v on 13, u data on 14, t data on 15
-
-! Case 1: regular grid to dipole
-! Case 2: dipole to regular grid
-! Case 3: regular grid to regular grid with same grid as dipole in SH
-! Case 4: regular grid with same grid as dipole in SH to regular grid
-! Case 5: regular grid with same grid as dipole in SH to dipole
-! Case 6: dipole to regular grid with same grid as dipole in SH
-
-! do not let the standard init run
-module_initialized = .true.
-
-if(test_casenum == 1 .or. test_casenum == 3) then
-   ! Case 1 or 3: read in from regular grid 
-   open(unit=12, position='rewind', action='read', file='regular_grid_u')
-   open(unit=13, position='rewind', action='read', file='regular_grid_t')
-   open(unit=14, position='rewind', action='read', file='regular_grid_u_data')
-   open(unit=15, position='rewind', action='read', file='regular_grid_t_data')
-
-else if(test_casenum == 4 .or. test_casenum == 5) then
-   ! Case 3 or 4: read regular with same grid as dipole in SH
-   open(unit=12, position='rewind', action='read', file='regular_griddi_u')
-   open(unit=13, position='rewind', action='read', file='regular_griddi_t')
-   open(unit=14, position='rewind', action='read', file='regular_griddi_u_data')
-   open(unit=15, position='rewind', action='read', file='regular_griddi_t_data')
-
-else if(test_casenum == 2 .or. test_casenum == 6) then
-   ! Case 5: read in from dipole grid
-   open(unit=12, position='rewind', action='read', file='dipole_grid_u')
-   open(unit=13, position='rewind', action='read', file='dipole_grid_t')
-   open(unit=14, position='rewind', action='read', file='dipole_grid_u_data')
-   open(unit=15, position='rewind', action='read', file='dipole_grid_t_data')
-endif
-
-! Get the size of the grid from the input u and t files
-read(12, *) nx, ny
-read(13, *) nx_temp, ny_temp
-if(nx /= nx_temp .or. ny /= ny_temp) then
-   write(msgstring,*)'mismatch nx,nx_temp ',nx,nx_temp,' or ny,ny_temp',ny,ny_temp
-   call error_handler(E_ERR,'test_interpolation',msgstring,source,revision,revdate)
-endif
-
-! Allocate stuff for the first grid (the one being interpolated from)
-allocate(ulon(nx, ny), ulat(nx, ny), tlon(nx, ny), tlat(nx, ny))
-allocate( kmt(nx, ny),  kmu(nx, ny))
-allocate(reg_u_data(nx, ny), reg_t_data(nx, ny))
-! The Dart format 1d data arrays
-allocate(reg_u_x(nx*ny), reg_t_x(nx*ny))
-
-do imain = 1, nx
-   do jmain = 1, ny
-      read(12, *) ti, tj, ulon(imain, jmain), ulat(imain, jmain)
-      read(13, *) ti, tj, tlon(imain, jmain), tlat(imain, jmain)
-      read(14, *) ti, tj, reg_u_data(imain, jmain)
-      read(15, *) ti, tj, reg_t_data(imain, jmain)
-   enddo
-enddo
-
-! Load into 1D dart data array
-index = 0
-do jmain = 1, ny
-   do imain = 1, nx
-      index = index + 1
-      reg_u_x(index) = reg_u_data(imain, jmain)
-      reg_t_x(index) = reg_t_data(imain, jmain)
-   enddo
-enddo
-
-! dummy out vertical; let height always = 1 and allow
-! all grid points to be processed.
-kmt = 2
-kmu = 2
-height = 1
-
-! Initialize the interpolation data structure for this grid.
-call init_interp()
-
-! Now read in the points for the output grid
-! Case 1: regular grid to dipole
-! Case 2: dipole to regular grid
-! Case 3: regular grid to regular grid with same grid as dipole in SH
-! Case 4: regular grid with same grid as dipole in SH to regular grid
-! Case 5: regular grid with same grid as dipole in SH to dipole 
-! Case 6: dipole to regular grid with same grid as dipole in SH
-if(test_casenum == 1 .or. test_casenum == 5) then
-   ! Output to dipole grid
-   open(unit=22, position='rewind', action='read',  file='dipole_grid_u')
-   open(unit=23, position='rewind', action='read',  file='dipole_grid_t')
-   open(unit=24, position='rewind', action='write', file='dipole_grid_u_data.out')
-   open(unit=25, position='rewind', action='write', file='dipole_grid_t_data.out')
-
-else if(test_casenum == 2 .or. test_casenum == 4) then
-   ! Output to regular grid
-   open(unit=22, position='rewind', action='read',  file='regular_grid_u')
-   open(unit=23, position='rewind', action='read',  file='regular_grid_t')
-   open(unit=24, position='rewind', action='write', file='regular_grid_u_data.out')
-   open(unit=25, position='rewind', action='write', file='regular_grid_t_data.out')
-
-else if(test_casenum == 3 .or. test_casenum == 6) then
-   ! Output to regular grid with same grid as dipole in SH
-   open(unit=22, position='rewind', action='read',  file='regular_griddi_u')
-   open(unit=23, position='rewind', action='read',  file='regular_griddi_t')
-   open(unit=24, position='rewind', action='write', file='regular_griddi_u_data.out')
-   open(unit=25, position='rewind', action='write', file='regular_griddi_t_data.out')
-endif
-
-read(22, *) dnx, dny
-read(23, *) dnx_temp, dny_temp
-if(dnx /= dnx_temp .or. dny /= dny_temp) then
-   write(msgstring,*)'mismatch dnx,dnx_temp ',dnx,dnx_temp,' or dny,dny_temp',dny,dny_temp
-   call error_handler(E_ERR,'test_interpolation',msgstring,source,revision,revdate)
-endif
-
-allocate(dulon(dnx, dny), dulat(dnx, dny), dtlon(dnx, dny), dtlat(dnx, dny))
-allocate(dipole_u(dnx, dny), dipole_t(dnx, dny))
-
-dipole_u = 0.0_r8   ! just put some dummy values in to make sure they get changed.
-dipole_t = 0.0_r8   ! just put some dummy values in to make sure they get changed.
-
-do imain = 1, dnx
-do jmain = 1, dny
-   read(22, *) ti, tj, dulon(imain, jmain), dulat(imain, jmain)
-   read(23, *) ti, tj, dtlon(imain, jmain), dtlat(imain, jmain)
-enddo
-enddo
-
-do imain = 1, dnx
-   do jmain = 1, dny
-      ! Interpolate U from the first grid to the second grid
-
-      call lon_lat_interpolate(reg_u_x, dulon(imain, jmain), dulat(imain, jmain), &
-         KIND_U_CURRENT_COMPONENT, height, dipole_u(imain, jmain), istatus)
-
-      if ( istatus /= 0 ) then
-         write(msgstring,'(''cell'',i4,i4,1x,f12.8,1x,f12.8,'' U interp failed - code '',i4)') &
-              imain, jmain, dulon(imain, jmain), dulat(imain, jmain), istatus
-         call error_handler(E_MSG,'test_interpolation',msgstring,source,revision,revdate)
-      endif
-
-      write(24, *) dulon(imain, jmain), dulat(imain, jmain), dipole_u(imain, jmain)
-
-      ! Interpolate U from the first grid to the second grid
-
-      call lon_lat_interpolate(reg_t_x, dtlon(imain, jmain), dtlat(imain, jmain), &
-         KIND_POTENTIAL_TEMPERATURE, height, dipole_t(imain, jmain), istatus)
-
-      if ( istatus /= 0 ) then
-         write(msgstring,'(''cell'',i4,i4,1x,f12.8,1x,f12.8,'' T interp failed - code '',i4)') &
-              imain,jmain, dtlon(imain, jmain), dtlat(imain, jmain), istatus
-         call error_handler(E_MSG,'test_interpolation',msgstring,source,revision,revdate)
-      endif
-
-      write(25, *) dtlon(imain, jmain), dtlat(imain, jmain), dipole_t(imain, jmain)
-
-   enddo
-enddo
-
-end subroutine pop_test_interpolation
-
-!------------------------------------------------------------------
-
-subroutine compute_temperature(x, llon, llat, lheight, interp_val, istatus)
- real(r8), intent(in)  :: x(:), llon, llat, lheight
- real(r8), intent(out) :: interp_val
- integer,  intent(out) :: istatus
-
-! use potential temp, depth, and salinity to compute a sensible (in-situ)
-! temperature
-
-integer  :: hstatus, hgt_bot, hgt_top
-real(r8) :: hgt_fract, salinity_val, potential_temp, pres_val
-real(r8) :: pres_bot, pres_top
-
-interp_val = MISSING_R8
-istatus = 99
-
-! Get the bounding vertical levels and the fraction between bottom and top
-call height_bounds(lheight, Nz, ZC, hgt_bot, hgt_top, hgt_fract, hstatus)
-if(hstatus /= 0) then
-   istatus = 12
-   return
-endif
-
-
-! salinity - in msu (kg/kg).  converter will want psu (g/kg).
-call do_interp(x, start_index(S_index), hgt_bot, hgt_top, hgt_fract, llon, llat, &
-               KIND_SALINITY, salinity_val, istatus)
-if(istatus /= 0) return
-if (debug > 8) print *, 'salinity: ', salinity_val
-
-! potential temperature - degrees C.
-call do_interp(x, start_index(T_index), hgt_bot, hgt_top, hgt_fract, llon, llat, &
-               KIND_POTENTIAL_TEMPERATURE, potential_temp, istatus)
-if(istatus /= 0) return
-if (debug > 8) print *, 'potential temp: ', potential_temp
-
-! compute pressure at location between given levels.  these values are in bars;
-! the convert routine wants decibars as pressure input.  hgt_fract is 0 at bottom, 1 at top
-pres_val = pressure(hgt_bot) + hgt_fract * (pressure(hgt_top) - pressure(hgt_bot))
-if (debug > 8) then
-   print *, 'local pressure: ', pres_val
-   print *, 'bot, top, press: ', hgt_bot, pressure(hgt_bot), &
-                                 hgt_top, pressure(hgt_top), pres_val
-endif
-
-! and finally, convert to sensible (in-situ) temperature.
-! potential temp in degrees C, pressure in decibars, salinity in psu or pss (g/kg).
-call insitu_temp(potential_temp, salinity_val*1000.0_r8, pres_val*10.0_r8, interp_val)
-if (debug > 2) print *, 's,pt,pres,t: ', salinity_val, potential_temp, pres_val, interp_val
-
-end subroutine compute_temperature
-
-!------------------------------------------------------------------
-
-subroutine do_interp(x, base_offset, hgt_bot, hgt_top, hgt_fract, &
-                     llon, llat, obs_type, interp_val, istatus)
- real(r8),  intent(in) :: x(:)
- integer,   intent(in) :: base_offset, hgt_bot, hgt_top
- real(r8),  intent(in) :: hgt_fract, llon, llat
- integer,   intent(in) :: obs_type
- real(r8), intent(out) :: interp_val
- integer,  intent(out) :: istatus
- 
-! do a 2d horizontal interpolation for the value at the bottom level, 
-! then again for the top level, then do a linear interpolation in the 
-! vertical to get the final value.
-
-integer  :: offset
-real(r8) :: bot_val, top_val
-
-! Find the base location for the bottom height and interpolate horizontally 
-!  on this level.  Do bottom first in case it is below the ocean floor; can
-!  avoid the second horizontal interpolation.
-offset = base_offset + (hgt_bot - 1) * nx * ny
-if (debug > 6) &
-   print *, 'bot, field, abs offset: ', hgt_bot, base_offset, offset
-
-call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_bot, bot_val, istatus)
-! Failed istatus from interpolate means give up
-if(istatus /= 0) return
-if (debug > 6) &
-   print *, 'bot_val = ', bot_val
-
-! Find the base location for the top height and interpolate horizontally 
-!  on this level.
-offset = base_offset + (hgt_top - 1) * nx * ny
-if (debug > 6) &
-   print *, 'top, field, abs offset: ', hgt_top, base_offset, offset
-
-call lon_lat_interpolate(x(offset:), llon, llat, obs_type, hgt_top, top_val, istatus)
-! Failed istatus from interpolate means give up
-if(istatus /= 0) return
-if (debug > 6) &
-   print *, 'top_val = ', top_val
-
-! Then weight them by the vertical fraction and return
-interp_val = bot_val + hgt_fract * (top_val - bot_val)
-if (debug > 2) print *, 'do_interp: interp val = ',interp_val
-
-end subroutine do_interp
-
-!------------------------------------------------------------------
-
-subroutine insitu_temp(potemp, s, lpres, insitu_t)
- real(r8), intent(in)  :: potemp, s, lpres
- real(r8), intent(out) :: insitu_t
-
-! CODE FROM POP MODEL -
-! nsc 1 nov 2012:  i have taken the original subroutine with call:
-!  subroutine dpotmp(press,temp,s,rp,potemp)
-! and removed the original 'press' argument (setting it to 0.0 below)
-! and renamed temp -> potemp, and potemp -> insitu_t
-! i also reordered the args to be a bit more logical.  now you specify:
-! potential temp, salinity, local pressure in decibars, and you get
-! back in-situ temperature (called sensible temperature in the atmosphere;
-! what a thermometer would measure).  the original (F77 fixed format) code
-! had a computed goto which is deprecated/obsolete.  i replaced it with 
-! a set of 'if() then else if()' lines.  i did try to not alter the original
-! code so much it wasn't recognizable anymore.
-!
-!  aliciak note: rp = 0 and press = local pressure as function of depth
-!  will return potemp given temp.
-!  the trick here that if you make rp = local pressure and press = 0.0, 
-!  and put potemp in the "temp" variable , it will return insitu temp in the 
-!  potemp variable.
-
-! an example figure of the relationship of potential temp and in-situ temp
-! at depth:  http://oceanworld.tamu.edu/resources/ocng_textbook/chapter06/chapter06_05.htm
-! see the 'potential temperature' section (note graph starts at -1000m)
-
-!     title:
-!     *****
-
-!       insitu_temp  -- calculate sensible (in-situ) temperature from 
-!                       local pressure, salinity, and potential temperature
-
-!     purpose:
-!     *******
-
-!       to calculate sensible temperature, taken from a converter that
-!       went from sensible/insitu temperature to potential temperature
-!
-!       ref: N.P. Fofonoff
-!            Deep Sea Research
-!            in press Nov 1976
-
-!     arguments:
-!     **********
-
-!       potemp     -> potential temperature in celsius degrees
-!       s          -> salinity pss 78
-!       lpres      -> local pressure in decibars
-!       insitu_t   <- in-situ (sensible) temperature (deg c)
-
-
-!     local variables:
-!     *********
-
-integer  :: i,j,n
-real(r8) :: dp,p,q,r1,r2,r3,r4,r5,s1,t,x
-
-!     code:
-!     ****
-
-      s1 = s - 35.0_r8
-      p  = 0.0_r8
-      t  = potemp
-
-      dp = lpres - p
-      n  = int (abs(dp)/1000.0_r8) + 1
-      dp = dp/n
-
-      do i=1,n
-         do j=1,4
-
-            r1 = ((-2.1687e-16_r8 * t + 1.8676e-14_r8) * t - 4.6206e-13_r8) * p
-            r2 = (2.7759e-12_r8*t - 1.1351e-10_r8) * s1
-            r3 = ((-5.4481e-14_r8 * t + 8.733e-12_r8) * t - 6.7795e-10_r8) * t
-            r4 = (r1 + (r2 + r3 + 1.8741e-8_r8)) * p + (-4.2393e-8_r8 * t+1.8932e-6_r8) * s1
-            r5 = r4 + ((6.6228e-10_r8 * t-6.836e-8_r8) * t + 8.5258e-6_r8) * t + 3.5803e-5_r8
-
-            x  = dp*r5
-
-            if (j == 1) then
-               t = t + 0.5_r8 * x
-               q = x
-               p = p + 0.5_r8 * dp
-          
-            else if (j == 2) then
-               t = t + 0.29298322_r8 * (x-q)
-               q = 0.58578644_r8 * x + 0.121320344_r8 * q
-   
-            else if (j == 3) then
-               t = t + 1.707106781_r8 * (x-q)
-               q = 3.414213562_r8*x - 4.121320344_r8*q
-               p = p + 0.5_r8*dp
-
-            else ! j must == 4
-               t = t + (x - 2.0_r8 * q) / 6.0_r8
-
-            endif
-   
-         enddo ! j loop
-      enddo ! i loop
-
-      insitu_t = t
-
-if (debug > 2) print *, 'potential temp, salinity, local pressure -> sensible temp'
-if (debug > 2) print *, potemp, s, lpres, insitu_t
-
-!       potemp     -> potential temperature in celsius degrees
-!       s          -> salinity pss 78
-!       lpres      -> local pressure in decibars
-!       insitu_t   <- in-situ (sensible) temperature (deg c)
-
-end subroutine insitu_temp
-
-!------------------------------------------------------------------
-
-subroutine dpth2pres(nd, depth, pressure)
- integer,  intent(in)  :: nd
- real(r8), intent(in)  :: depth(nd)
- real(r8), intent(out) :: pressure(nd)
-
-!  description:
-!  this function computes pressure in bars from depth in meters
-!  using a mean density derived from depth-dependent global 
-!  average temperatures and salinities from levitus 1994, and 
-!  integrating using hydrostatic balance.
-! 
-!  references:
-! 
-!  levitus, s., r. burgett, and t.p. boyer, world ocean atlas 
-!  volume 3: salinity, noaa atlas nesdis 3, us dept. of commerce, 1994.
-! 
-!  levitus, s. and t.p. boyer, world ocean atlas 1994, volume 4:
-!  temperature, noaa atlas nesdis 4, us dept. of commerce, 1994.
-! 
-!  dukowicz, j. k., 2000: reduction of pressure and pressure
-!  gradient errors in ocean simulations, j. phys. oceanogr., submitted.
-
-!  input parameters:
-!  nd     - size of arrays
-!  depth  - depth in meters. no units check is made
-
-!  output parameters:
-!  pressure - pressure in bars 
-
-!  local variables & parameters:
-integer :: n
-real(r8), parameter :: c1 = 1.0_r8
-
-! -----------------------------------------------------------------------
-!  convert depth in meters to pressure in bars
-! -----------------------------------------------------------------------
-
-      do n=1,nd
-         pressure(n) = 0.059808_r8*(exp(-0.025_r8*depth(n)) - c1)  &
-                     + 0.100766_r8*depth(n) + 2.28405e-7_r8*depth(n)**2
-      end do
-
-if (debug > 2 .and. do_output()) then
-   print *, 'depth->pressure conversion table.  cols are: N, depth(m), pressure(bars)'
-   do n=1,nd
-      print *, n, depth(n), pressure(n)
-   enddo
-endif
-
-end subroutine dpth2pres
-
-!------------------------------------------------------------------
-!------------------------------------------------------------------
-
-
-!------------------------------------------------------------------
-! End of model_mod
-!------------------------------------------------------------------
-
-end module pop_model_mod
-
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$

Deleted: DART/branches/rma_trunk/models/POP/pop_to_dart.f90
===================================================================
--- DART/branches/rma_trunk/models/POP/pop_to_dart.f90	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/pop_to_dart.f90	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,107 +0,0 @@
-! DART software - Copyright 2004 - 2013 UCAR. This open source software is
-! provided by UCAR, "as is", without charge, subject to all terms of use at
-! http://www.image.ucar.edu/DAReS/DART/DART_download
-!
-! $Id$
-
-program pop_to_dart
-
-!----------------------------------------------------------------------
-! purpose: interface between POP and DART
-!
-! method: Read POP "restart" files of model state
-!         Reform fields into a DART state vector (control vector).
-!         Write out state vector in "proprietary" format for DART.
-!         The output is a "DART restart file" format.
-! 
-! USAGE:  The POP filename is read from the pop_in namelist
-!         <edit pop_to_dart_output_file in input.nml:pop_to_dart_nml>
-!         pop_to_dart
-!
-! author: Tim Hoar 6/24/09
-!----------------------------------------------------------------------
-
-use        types_mod, only : r8
-use    utilities_mod, only : initialize_utilities, finalize_utilities, &
-                             find_namelist_in_file, check_namelist_read
-use        model_mod, only : restart_file_to_sv, static_init_model, &
-                             get_model_size, get_pop_restart_filename
-use state_vector_io_mod, only : awrite_state_restart, open_restart_write, close_restart
-use time_manager_mod, only : time_type, print_time, print_date
-
-use netcdf
-implicit none
-
-! version controlled file description for error handling, do not edit
-character(len=256), parameter :: source   = &
-   "$URL$"
-character(len=32 ), parameter :: revision = "$Revision$"
-character(len=128), parameter :: revdate  = "$Date$"
-
-!-----------------------------------------------------------------------
-! namelist parameters with default values.
-!-----------------------------------------------------------------------
-
-character (len = 128) :: pop_to_dart_output_file  = 'dart_ics'
-
-namelist /pop_to_dart_nml/ pop_to_dart_output_file
-
-!----------------------------------------------------------------------
-! global storage
-!----------------------------------------------------------------------
-
-integer               :: io, iunit, x_size
-type(time_type)       :: model_time
-real(r8), allocatable :: statevector(:)
-character (len = 128) :: pop_restart_filename = 'no_pop_restart_filename' 
-
-!----------------------------------------------------------------------
-
-call initialize_utilities(progname='pop_to_dart')
-
-!----------------------------------------------------------------------
-! Call model_mod:static_init_model(), which reads the namelists
-! to set calendar type, starting date, deltaT, etc.
-!----------------------------------------------------------------------
-
-call static_init_model()
-
-!----------------------------------------------------------------------
-! Read the namelist to get the input and output filenames.
-!----------------------------------------------------------------------
-
-call find_namelist_in_file("input.nml", "pop_to_dart_nml", iunit)
-read(iunit, nml = pop_to_dart_nml, iostat = io)
-call check_namelist_read(iunit, io, "pop_to_dart_nml") ! closes, too.
-
-call get_pop_restart_filename( pop_restart_filename )
-
-write(*,*)
-write(*,'(''pop_to_dart:converting POP restart file '',A, &
-      &'' to DART file '',A)') &
-       trim(pop_restart_filename), trim(pop_to_dart_output_file)
-
-!----------------------------------------------------------------------
-! Now that we know the names, get to work.
-!----------------------------------------------------------------------
-
-x_size = get_model_size()
-allocate(statevector(x_size))
-call restart_file_to_sv(pop_restart_filename, statevector, model_time) 
-
-iunit = open_restart_write(pop_to_dart_output_file)
-
-call awrite_state_restart(model_time, statevector, iunit)
-call close_restart(iunit)
-
-call print_date(model_time, str='pop_to_dart:POP  model date')
-call print_time(model_time, str='pop_to_dart:DART model time')
-call finalize_utilities('pop_to_dart')
-
-end program pop_to_dart
-
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$

Deleted: DART/branches/rma_trunk/models/POP/pop_to_dart.html
===================================================================
--- DART/branches/rma_trunk/models/POP/pop_to_dart.html	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/pop_to_dart.html	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,223 +0,0 @@
-<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
-"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
-<HTML>
-<HEAD>
-<TITLE>program pop_to_dart</TITLE>
-<link rel="stylesheet" type="text/css" href="../../doc/html/doc.css" />
-<link href="../../doc/images/dart.ico" rel="shortcut icon" />
-</HEAD>
-<BODY>
-<A NAME="TOP"></A>
-
-<H1>PROGRAM <em class=program>pop_to_dart</em></H1>
-
-<table border=0 summary="" cellpadding=5>
-<tr>
-    <td valign=middle>
-    <img src="../../doc/images/Dartboard7.png" alt="DART project logo" height=70 />
-    </td>
-    <td>
-       <P>Jump to <a href="../../index.html">DART Documentation Main Index</a><br />
-          <small><small>version information for this file: <br />
-          <!-- version tag follows, do not edit -->
-          $Id$</small></small>
-       </P></td>
-</tr>
-</table>
-
-<A HREF="#Namelist">NAMELIST</A> /
-<A HREF="#Modules">MODULES</A> /
-<A HREF="#FilesUsed">FILES</A> /
-<A HREF="#References">REFERENCES</A> /
-<A HREF="#Errors">ERRORS</A> /
-<A HREF="#FuturePlans">PLANS</A> /
-<A HREF="#Legalese">TERMS OF USE</A>
-
-<H2>Overview</H2>
-
-<P>
-   <em class=program>pop_to_dart</em> is the program that reads a POP
-   restart file (usually <em class=file>pop.r.nc</em>) and creates
-   a DART output/restart file 
-   (e.g. <em class=file>perfect_ics, filter_ics, ... </em>).
-   <br />
-   <br />
-   The following variables are extracted from the POP netCDF restart 
-   files and are conveyed to DART:
-   <em class=code>SALT_CUR</em>,
-   <em class=code>TEMP_CUR</em>,
-   <em class=code>UVEL_CUR</em>,
-   <em class=code>VVEL_CUR</em>, and
-   <em class=code>PSURF_CUR</em>.
-   <br />
-   <br />
-   Conditions required for successful execution of <em class=program>pop_to_dart</em>:
-</P>
-
-<UL>
-   <LI>a valid <em class=file>input.nml</em> namelist file for DART</LI>
-   <LI>a valid <em class=file>pop_in</em> namelist file for POP</LI>
-   <LI>the POP geometry files mentioned in the <em class=file>pop_in</em> namelist file</LI>
-   <LI>a POP restart file (<em class=file>pop.r.nc</em>).</LI>
-</UL>
-
-<P>
-Since this program is called repeatedly for every ensemble member,
-we have found it convenient to link the POP restart files
-to the default input filename (<em class=file>pop.r.nc</em>).
-The default DART output filename is <em class=file>dart_ics</em> -
-this may be moved or linked as necessary.
-</P>
-
-<!--==================================================================-->
-<!--=================== DESCRIPTION OF A NAMELIST ====================-->
-<!--==================================================================-->
-
-<A NAME="Namelist"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>NAMELIST</H2>
-<P>
-This namelist is read from the file <em class=file>input.nml</em>.
-Namelists start with an ampersand
-'&amp;' and terminate with a slash '/'.
-Character strings that contain a '/' must be
-enclosed in quotes to prevent them from 
-prematurely terminating the namelist.
-</P>
-
-<div class=namelist>
-<pre>
-&amp;pop_to_dart_nml
-   pop_to_dart_output_file = 'dart_ics'  
-/
-</pre>
-</div>
-
-<br />
-<br />
-
-<div>
-<TABLE border=0 cellpadding=10 width=100% summary='namelist description'>
-<THEAD align=left>
-<TR><TH> Item </TH>
-    <TH> Type </TH>
-    <TH> Description </TH> </TR>
-</THEAD>
-
-<TBODY valign=top>
-
-<TR><TD>pop_to_dart_output_file   </TD>
-    <TD>character(len=128)   </TD>
-    <TD>The name of the DART file containing the model state
-derived from the POP restart file.
-</TD></TR>
-
-</TBODY> 
-</TABLE>
-</div>
-
-<br />
-<br />
-
-<!--==================================================================-->
-
-<A NAME="Modules"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>MODULES USED</H2>
-<PRE>
-assim_model_mod
-dart_pop_mod
-location_mod
-model_mod
-null_mpi_utilities_mod
-obs_kind_mod
-random_nr_mod
-random_seq_mod
-time_manager_mod
-types_mod
-utilities_mod
-</PRE>
-
-<!--==================================================================-->
-<!-- Describe the Files Used by this module.                          -->
-<!--==================================================================-->
-
-<A NAME="FilesUsed"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>FILES Read</H2>
-<UL><LI>POP restart file; <em class=file>pop.r.nc</em></LI>
-    <LI>POP namelist file; <em class=file>pop_in</em></LI>
-    <LI>POP geometry files specified in <em class=file>pop_in</em></LI>
-    <LI>DART namelist file; <em class=file>input.nml</em></LI>
-</UL>
-
-<H2>FILES Written</H2>
-<UL><LI>DART initial conditions/restart file; e.g. <em class=file>filter_ic</em></LI>
-</UL>
-
-<!--==================================================================-->
-<!-- Cite references, if need be.                                     -->
-<!--==================================================================-->
-
-<A NAME="References"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>REFERENCES</H2>
-<P>
-none
-</P>
-
-<!--==================================================================-->
-<!-- Describe all the error conditions and codes.                     -->
-<!--==================================================================-->
-
-<A NAME="Errors"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>ERROR CODES and CONDITIONS</H2>
-<P>
-none - all error messages come from modules that have their own documentation.
-</P>
-
-<H2>KNOWN BUGS</H2>
-<P>
-none
-</P>
-
-<!--==================================================================-->
-<!-- Describe Future Plans.                                           -->
-<!--==================================================================-->
-
-<A NAME="FuturePlans"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>FUTURE PLANS</H2>
-<P>
-None.
-</P>
-
-<!--==================================================================-->
-<!-- Legalese & Metadata                                              -->
-<!--==================================================================-->
-
-<A NAME="Legalese"></A>
-<div class="top">[<a href="#">top</a>]</div><hr />
-<H2>Terms of Use</H2>
-
-<P>
-DART software - Copyright 2004 - 2013 UCAR.<br />
-This open source software is provided by UCAR, "as is",<br />
-without charge, subject to all terms of use at<br />
-<a href="http://www.image.ucar.edu/DAReS/DART/DART_download">
-http://www.image.ucar.edu/DAReS/DART/DART_download</a>
-</P>
-
-<TABLE border=0 cellpadding=0 width=100% summary="">
-<TR><TD valign=top>Contact:       </TD><TD> DART core group   </TD></TR>
-<TR><TD valign=top>Revision:      </TD><TD> $Revision$ </TD></TR>
-<TR><TD valign=top>Source:        </TD><TD> $URL$ </TD></TR>
-<TR><TD valign=top>Change Date:   </TD><TD> $Date$ </TD></TR>
-<TR><TD valign=top>Change&nbsp;history:&nbsp;</TD><TD> try "svn&nbsp;log" or "svn&nbsp;diff" </TD></TR>
-</TABLE>
-
-<!--==================================================================-->
-
-</BODY>
-</HTML>

Deleted: DART/branches/rma_trunk/models/POP/pop_to_dart.nml
===================================================================
--- DART/branches/rma_trunk/models/POP/pop_to_dart.nml	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/pop_to_dart.nml	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,3 +0,0 @@
-&pop_to_dart_nml
-   pop_to_dart_output_file = 'dart_ics'  /
-

Copied: DART/branches/rma_trunk/models/POP/shell_scripts/stitch_post.sh (from rev 10736, DART/branches/rma_trunk/models/POP/stitch_post.sh)
===================================================================
--- DART/branches/rma_trunk/models/POP/shell_scripts/stitch_post.sh	                        (rev 0)
+++ DART/branches/rma_trunk/models/POP/shell_scripts/stitch_post.sh	2016-11-11 17:21:36 UTC (rev 10740)
@@ -0,0 +1,61 @@
+#!/bin/bash 
+
+#Aim: stitch together filter output to make diagnostic files
+module load nco
+
+echo -n "starting time "
+date +"%T"
+
+#Get new files, remove old ones
+rm *.nc 
+cp ../Posterior_Diag.nc .
+cp ../Output/*.nc .
+
+# Make copy the record dimension
+ncpdq -O -a copy,time Posterior_Diag.nc Posterior_Diag.nc
+
+export mean='mean_copy_d01.nc'
+export sd='sd_copy_d01.nc'
+export inf_mean='post_inflate_restart01_mean.nc'
+export inf_sd='post_inflate_restart01_sd.nc'
+
+# create time dimension in POP output
+ncecat -O -u time $mean $mean &
+ncecat -O -u time $sd $sd &
+ncecat -O -u time $inf_mean $inf_mean &
+ncecat -O -u time $inf_sd $inf_sd &
+
+wait
+
+# create copy record dimension in POP output
+ncecat -O -u copy $mean $mean &
+ncecat -O -u copy $sd $sd &
+ncecat -O -u copy $inf_mean $inf_mean &
+ncecat -O -u copy $inf_sd $inf_sd &
+
+wait
+
+## rename variables
+for var in 'SALT' 'TEMP' 'UVEL' 'VVEL' 'PSURF' 
+do
+
+   ncrename -v $var'_CUR',$var $mean &
+   ncrename -v $var'_CUR',$var $sd &
+   ncrename -v $var'_CUR',$var $inf_mean &
+   ncrename -v $var'_CUR',$var $inf_sd &
+
+   wait
+ 
+   # concatenate the 4 files into Posterior_Diag
+   # append
+   ncrcat -A -v $var $mean $sd $inf_mean $inf_sd Posterior_Diag.nc  
+   echo -n "done variable " $var "time "
+   date +"%T"   
+done
+
+# switch time and copy dimensions back
+ncpdq -O -a time,copy Posterior_Diag.nc Posterior_Diag.nc
+echo -n "finished " 
+date +"%T"
+
+#mv Posterior_Diag.nc Posterior_Diag.nc.full

Copied: DART/branches/rma_trunk/models/POP/shell_scripts/stitch_prior.sh (from rev 10736, DART/branches/rma_trunk/models/POP/stitch_prior.sh)
===================================================================
--- DART/branches/rma_trunk/models/POP/shell_scripts/stitch_prior.sh	                        (rev 0)
+++ DART/branches/rma_trunk/models/POP/shell_scripts/stitch_prior.sh	2016-11-11 17:21:36 UTC (rev 10740)
@@ -0,0 +1,61 @@
+#!/bin/bash 
+
+#Aim: stitch together filter output to make diagnostic files
+module load nco
+
+echo -n "starting time "
+date +"%T"
+
+#Get new files, remove old ones
+rm *.nc 
+cp ../Prior_Diag.nc .
+cp ../Output/*.nc .
+
+# Make copy the record dimension
+ncpdq -O -a copy,time Prior_Diag.nc Prior_Diag.nc
+
+export mean='prior_diag_mean01.nc'
+export sd='prior_diag_sd01.nc'
+export inf_mean='prior_diag_inf_mean01.nc'
+export inf_sd='prior_diag_inf_sd01.nc'
+
+# create time dimension in POP output
+ncecat -O -u time $mean $mean &
+ncecat -O -u time $sd $sd &
+ncecat -O -u time $inf_mean $inf_mean &
+ncecat -O -u time $inf_sd $inf_sd &
+
+wait
+
+# create copy record dimension in POP output
+ncecat -O -u copy $mean $mean &
+ncecat -O -u copy $sd $sd &
+ncecat -O -u copy $inf_mean $inf_mean &
+ncecat -O -u copy $inf_sd $inf_sd &
+
+wait
+
+## rename variables
+for var in 'SALT' 'TEMP' 'UVEL' 'VVEL' 'PSURF' 
+do
+
+   ncrename -v $var'_CUR',$var $mean &
+   ncrename -v $var'_CUR',$var $sd &
+   ncrename -v $var'_CUR',$var $inf_mean &
+   ncrename -v $var'_CUR',$var $inf_sd &
+
+   wait
+ 
+   # concatenate the 4 files into Prior_Diag
+   # append
+   ncrcat -A -v $var $mean $sd $inf_mean $inf_sd Prior_Diag.nc  
+   echo -n "done variable " $var "time "
+   date +"%T"   
+done
+
+# switch time and copy dimensions back
+ncpdq -O -a time,copy Prior_Diag.nc Prior_Diag.nc
+echo -n "finished " 
+date +"%T"
+
+#mv Prior_Diag.nc Prior_Diag.nc.full

Deleted: DART/branches/rma_trunk/models/POP/stitch_post.sh
===================================================================
--- DART/branches/rma_trunk/models/POP/stitch_post.sh	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/stitch_post.sh	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,61 +0,0 @@
-#!/bin/bash 
-
-#Aim: stitch together filter output to make diagnostic files
-module load nco
-
-echo -n "starting time "
-date +"%T"
-
-#Get new files, remove old ones
-rm *.nc 
-cp ../Posterior_Diag.nc .
-cp ../Output/*.nc .
-
-# Make copy the record dimension
-ncpdq -O -a copy,time Posterior_Diag.nc Posterior_Diag.nc
-
-export mean='mean_copy_d01.nc'
-export sd='sd_copy_d01.nc'
-export inf_mean='post_inflate_restart01_mean.nc'
-export inf_sd='post_inflate_restart01_sd.nc'
-
-# create time dimension in POP output
-ncecat -O -u time $mean $mean &
-ncecat -O -u time $sd $sd &
-ncecat -O -u time $inf_mean $inf_mean &
-ncecat -O -u time $inf_sd $inf_sd &
-
-wait
-
-# create copy record dimension in POP output
-ncecat -O -u copy $mean $mean &
-ncecat -O -u copy $sd $sd &
-ncecat -O -u copy $inf_mean $inf_mean &
-ncecat -O -u copy $inf_sd $inf_sd &
-
-wait
-
-## rename variables
-for var in 'SALT' 'TEMP' 'UVEL' 'VVEL' 'PSURF' 
-do
-
-   ncrename -v $var'_CUR',$var $mean &
-   ncrename -v $var'_CUR',$var $sd &
-   ncrename -v $var'_CUR',$var $inf_mean &
-   ncrename -v $var'_CUR',$var $inf_sd &
-
-   wait
- 
-   # concatenate the 4 files into Posterior_Diag
-   # append
-   ncrcat -A -v $var $mean $sd $inf_mean $inf_sd Posterior_Diag.nc  
-   echo -n "done variable " $var "time "
-   date +"%T"   
-done
-
-# switch time and copy dimensions back
-ncpdq -O -a time,copy Posterior_Diag.nc Posterior_Diag.nc
-echo -n "finished " 
-date +"%T"
-
-#mv Posterior_Diag.nc Posterior_Diag.nc.full

Deleted: DART/branches/rma_trunk/models/POP/stitch_prior.sh
===================================================================
--- DART/branches/rma_trunk/models/POP/stitch_prior.sh	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/stitch_prior.sh	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,61 +0,0 @@
-#!/bin/bash 
-
-#Aim: stitch together filter output to make diagnostic files
-module load nco
-
-echo -n "starting time "
-date +"%T"
-
-#Get new files, remove old ones
-rm *.nc 
-cp ../Prior_Diag.nc .
-cp ../Output/*.nc .
-
-# Make copy the record dimension
-ncpdq -O -a copy,time Prior_Diag.nc Prior_Diag.nc
-
-export mean='prior_diag_mean01.nc'
-export sd='prior_diag_sd01.nc'
-export inf_mean='prior_diag_inf_mean01.nc'
-export inf_sd='prior_diag_inf_sd01.nc'
-
-# create time dimension in POP output
-ncecat -O -u time $mean $mean &
-ncecat -O -u time $sd $sd &
-ncecat -O -u time $inf_mean $inf_mean &
-ncecat -O -u time $inf_sd $inf_sd &
-
-wait
-
-# create copy record dimension in POP output
-ncecat -O -u copy $mean $mean &
-ncecat -O -u copy $sd $sd &
-ncecat -O -u copy $inf_mean $inf_mean &
-ncecat -O -u copy $inf_sd $inf_sd &
-
-wait
-
-## rename variables
-for var in 'SALT' 'TEMP' 'UVEL' 'VVEL' 'PSURF' 
-do
-
-   ncrename -v $var'_CUR',$var $mean &
-   ncrename -v $var'_CUR',$var $sd &
-   ncrename -v $var'_CUR',$var $inf_mean &
-   ncrename -v $var'_CUR',$var $inf_sd &
-
-   wait
- 
-   # concatenate the 4 files into Prior_Diag
-   # append
-   ncrcat -A -v $var $mean $sd $inf_mean $inf_sd Prior_Diag.nc  
-   echo -n "done variable " $var "time "
-   date +"%T"   
-done
-
-# switch time and copy dimensions back
-ncpdq -O -a time,copy Prior_Diag.nc Prior_Diag.nc
-echo -n "finished " 
-date +"%T"
-
-#mv Prior_Diag.nc Prior_Diag.nc.full

Deleted: DART/branches/rma_trunk/models/POP/work/broken_mkmf_model_mod_check
===================================================================
--- DART/branches/rma_trunk/models/POP/work/broken_mkmf_model_mod_check	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/work/broken_mkmf_model_mod_check	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,18 +0,0 @@
-#!/bin/csh
-#
-# DART software - Copyright 2004 - 2013 UCAR. This open source software is
-# provided by UCAR, "as is", without charge, subject to all terms of use at
-# http://www.image.ucar.edu/DAReS/DART/DART_download
-#
-# DART $Id$
-
-../../../mkmf/mkmf -p model_mod_check -t ../../../mkmf/mkmf.template   \
- -a "../../.." -w path_names_model_mod_check
-
-exit $status
-
-# <next few lines under version control, do not edit>
-# $URL$
-# $Revision$
-# $Date$
-

Deleted: DART/branches/rma_trunk/models/POP/work/mkmf_dart_to_pop
===================================================================
--- DART/branches/rma_trunk/models/POP/work/mkmf_dart_to_pop	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/work/mkmf_dart_to_pop	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,18 +0,0 @@
-#!/bin/csh
-#
-# DART software - Copyright 2004 - 2013 UCAR. This open source software is
-# provided by UCAR, "as is", without charge, subject to all terms of use at
-# http://www.image.ucar.edu/DAReS/DART/DART_download
-#
-# DART $Id$
-
-../../../mkmf/mkmf -p dart_to_pop -t ../../../mkmf/mkmf.template   \
- -a "../../.." path_names_dart_to_pop
-
-exit $status
-
-# <next few lines under version control, do not edit>
-# $URL$
-# $Revision$
-# $Date$
-

Deleted: DART/branches/rma_trunk/models/POP/work/mkmf_pop_to_dart
===================================================================
--- DART/branches/rma_trunk/models/POP/work/mkmf_pop_to_dart	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/work/mkmf_pop_to_dart	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,18 +0,0 @@
-#!/bin/csh
-#
-# DART software - Copyright 2004 - 2013 UCAR. This open source software is
-# provided by UCAR, "as is", without charge, subject to all terms of use at
-# http://www.image.ucar.edu/DAReS/DART/DART_download
-#
-# DART $Id$
-
-../../../mkmf/mkmf -p pop_to_dart -t ../../../mkmf/mkmf.template   \
- -a "../../.." path_names_pop_to_dart
-
-exit $status
-
-# <next few lines under version control, do not edit>
-# $URL$
-# $Revision$
-# $Date$
-

Deleted: DART/branches/rma_trunk/models/POP/work/path_names_dart_to_pop
===================================================================
--- DART/branches/rma_trunk/models/POP/work/path_names_dart_to_pop	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/work/path_names_dart_to_pop	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,22 +0,0 @@
-adaptive_inflate/adaptive_inflate_mod.f90
-assim_model/assim_model_mod.f90
-common/types_mod.f90
-distributed/distributed_state_mod.f90
-distributed/null_win_mod.f90
-ensemble_manager/ensemble_manager_mod.f90
-io/copies_on_off_mod.f90
-io/dart_time_io_mod.f90
-io/direct_netcdf_mod.f90
-io/io_filenames_mod.f90
-io/state_structure_mod.f90
-io/state_vector_io_mod.f90
-location/threed_sphere/location_mod.f90
-models/POP/dart_pop_mod.f90
-models/POP/dart_to_pop.f90
-models/POP/model_mod.f90
-mpi_utilities/null_mpi_utilities_mod.f90
-obs_kind/obs_kind_mod.f90
-random_seq/random_seq_mod.f90
-sort/sort_mod.f90
-time_manager/time_manager_mod.f90
-utilities/utilities_mod.f90

Deleted: DART/branches/rma_trunk/models/POP/work/path_names_pop_to_dart
===================================================================
--- DART/branches/rma_trunk/models/POP/work/path_names_pop_to_dart	2016-11-10 23:40:48 UTC (rev 10739)
+++ DART/branches/rma_trunk/models/POP/work/path_names_pop_to_dart	2016-11-11 17:21:36 UTC (rev 10740)
@@ -1,22 +0,0 @@
-adaptive_inflate/adaptive_inflate_mod.f90
-assim_model/assim_model_mod.f90
-common/types_mod.f90
-distributed/distributed_state_mod.f90
-distributed/null_win_mod.f90
-ensemble_manager/ensemble_manager_mod.f90
-io/copies_on_off_mod.f90
-io/dart_time_io_mod.f90
-io/direct_netcdf_mod.f90
-io/io_filenames_mod.f90
-io/state_structure_mod.f90
-io/state_vector_io_mod.f90
-location/threed_sphere/location_mod.f90
-models/POP/dart_pop_mod.f90
-models/POP/model_mod.f90
-models/POP/pop_to_dart.f90
-mpi_utilities/null_mpi_utilities_mod.f90
-obs_kind/obs_kind_mod.f90
-random_seq/random_seq_mod.f90
-sort/sort_mod.f90
-time_manager/time_manager_mod.f90
-utilities/utilities_mod.f90


More information about the Dart-dev mailing list