diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/.gitignore b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/.gitignore
new file mode 100644
index 0000000000..a9ae5ae27e
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/.gitignore
@@ -0,0 +1,10 @@
+GEBCO_08_v1.nc
+IBCAO_V3_500m_RR.grd
+bedmap2.nc
+DATA
+local
+tpxo7_atlas_netcdf.tar.Z
+input.nml
+MIDAS
+fre_nctools
+mosaic*
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/C180_mosaic/create_c180_mosaic_and_rivers.csh b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/C180_mosaic/create_c180_mosaic_and_rivers.csh
new file mode 100755
index 0000000000..9cc1fc1caa
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/C180_mosaic/create_c180_mosaic_and_rivers.csh
@@ -0,0 +1,19 @@
+#!/bin/csh 
+
+module load fre
+
+
+ln -s ../ocean_hgrid.nc .
+ln -s ../topog.nc .
+make_hgrid --grid_type gnomonic_ed --nlon 360  --grid_name C180_grid
+make_solo_mosaic --num_tiles 6 --dir . --mosaic_name C180_mosaic --tile_file C180_grid.tile1.nc,C180_grid.tile2.nc,C180_grid.tile3.nc,C180_grid.tile4.nc,C180_grid.tile5.nc,C180_grid.tile6.nc
+make_solo_mosaic --num_tiles 1 --dir ./ --mosaic_name ocean_mosaic --tile_file ocean_hgrid.nc --periodx 360.
+make_coupler_mosaic --atmos_mosaic C180_mosaic.nc --ocean_mosaic ocean_mosaic.nc --ocean_topog topog.nc --mosaic_name grid_spec
+rm -f hydrology*.nc
+
+#./run.cr_simple_hydrog.csh -f 0. -t 1.e-5 -m ../grid_spec.nc -o . -c compile
+#(mv work/*.nc .;foreach tile ( tile1 tile2 tile3 tile4 tile5 tile6 );mv hydrography.$tile.nc river_data.$tile.nc;end)
+(tar cvhf mosaic.tar C180_grid* C180_mosaic* lake_frac* land_mask* ocean_hgrid* ocean_mask* ocean_mosaic* river_data* topog* grid_spec.nc)
+(mkdir quick_mosaic;cd quick_mosaic;make_quick_mosaic --mosaic_name grid_spec --input_mosaic ../grid_spec.nc)
+(cd quick_mosaic;ln -s ../C180_grid.tile*.nc .;ln -s ../C180_mosaic.nc .;ln -s ../lake_frac* .;ln -s ../river_data* .;ln -s ../topog.nc .)
+(cd quick_mosaic;tar chvf ../quick_mosaic.tar C180_grid* C180_mosaic* lake_frac* land_mask* ocean_mask* river_data* topog* atmos_mosaic* land_mosaic* grid_spec.nc)
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/C180_mosaic/run.cr_simple_hydrog.csh b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/C180_mosaic/run.cr_simple_hydrog.csh
new file mode 100755
index 0000000000..5cd1116950
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/C180_mosaic/run.cr_simple_hydrog.csh
@@ -0,0 +1,286 @@
+#!/bin/csh -f
+#----------------------------------
+#PBS -N cr_simple_hydrog
+#PBS -l size=1
+#PBS -l walltime=02:00:00
+#PBS -r y
+#PBS -j oe
+#PBS -o $HOME
+#PBS -q batch
+#----------------------------------
+
+# =========================================================================
+#   script creates a simple hydrography, with no manual intervention,
+#     simplified lake fields, no elevation
+#    -- run river_regrid tool
+#    -- run post-processor on river_regrid output
+#    -- remove parallel rivers
+#    -- run post-processor again
+#    -- add GLCC waterbod fractions
+# =========================================================================
+
+#set echo on
+#wipeftmp
+
+if (`gfdl_platform` == "hpcs-csc") then
+    source $MODULESHOME/init/csh
+    module purge
+    module load gcp
+    module load netcdf/4.2
+    module load nco/4.3.1 
+    module load ifort/11.1.073
+    module load icc/11.1.073
+    module load mpich2/1.2.1p1
+    alias gmake gmake -j 8
+else
+   echo "ERROR: invalid platform"
+   exit
+endif
+
+set outdir = ""
+set argv = (`getopt  hc:o:f:t:m: $*`)
+
+while ("$argv[1]" != "--")
+    switch ($argv[1])
+        case -h:
+            set help; breaksw
+        case -o:
+            set outdir = $argv[2]; shift argv; breaksw
+        case -c:
+            set compile = $argv[2]; shift argv; breaksw
+        case -f:
+            set min_frac = $argv[2]; shift argv; breaksw
+        case -t:
+            set land_thresh = $argv[2]; shift argv; breaksw
+        case -m:
+            set mosaic_file = $argv[2]; shift argv; breaksw
+    endsw
+    shift argv
+end
+shift argv
+
+# argument error checking
+
+if (! $?min_frac) then
+   echo "ERROR: no argument given for min_frac (river_regrid)."
+   set help
+endif
+
+if (! $?land_thresh) then
+   echo "ERROR: no argument given for land_thresh (river_regrid)."
+   set help
+endif
+
+if (! $?mosaic_file) then
+   echo "ERROR: no argument given for mosaic_file."
+   set help
+endif
+
+if (! $?compile) then
+   echo "NOTE: no argument given for compile: default is to use previously compiled executables "
+   echo " "
+   set compile = " "
+endif
+
+
+# usage message
+if ($?help) then
+   set cmdname = `basename $0`
+   echo
+   echo "USAGE:  $cmdname [-f min_frac] [-t land_thresh] [-m mosaic] [-c compile] [-o outdir]"
+   echo
+   echo "         -f min_frac:    required, river-regrid, min land fraction; set to 0 to retain small "
+   echo "                         land fractions from gridspec"
+   echo "         -t land_thresh: required, river-regrid, land fraction exceeding (1-land_thresh) is set to 1;"
+   echo "                         recommend 1.e-5"
+   echo "         -m mosaic:      required, mosaic file location"
+   echo "         -c compile:     use '-c compile' to compile code (optional)"
+   echo "         -o outdir:      output directory (optional)"
+   echo
+   exit 1
+endif
+
+
+set home_dir         = /home/kap/lmdt/tools/simple_hydrog/
+set fms_code_dir     = /home/kap/lmdt/analysis/diag3_bronx3/shared/code_siena_201211
+
+set riv_regrd        = /home/Zhi.Liang/bin/tools_20130912/river_regrid              # set path to Zhi's river_regrid tool
+set postp_code_path  = $home_dir/code/cp_river_vars.f90
+set rmvpr_code_path  = $home_dir/code/rmv_parallel_rivers.f90
+set lakes_code_path  = $home_dir/code/path_names_lakes_siena201211
+
+set compile_template = $home_dir/code/TEMPLATE_siena_201211_b3
+
+# set path for the disaggregated, extended river network
+set river_network_ext_file = /archive/kap/lmdt/river_network/netcdf/0.5deg/disagg/river_network_mrg_0.5deg_ad3nov_fill_coast_auto1_0.125.nc
+
+# set path for GLCC data (waterbod)
+set glcc_file = /archive/pcm/land_data/cover_lad/gigbp2_0ll.nc
+
+
+
+if (! -d work) mkdir work
+set WORKDIR = `pwd`/work
+cd $WORKDIR
+
+mkdir -p EXEC/{postp,rmvpr,lakes}
+mkdir -p OUTPUT/{river_regrid,post_regrid,rmv_parallel_rivers,post_rmvp}
+mkdir DATA
+
+if ($compile == "compile") then
+# ------------------------------------------------
+#  COMPILE FORTRAN PROGRAMS
+# ------------------------------------------------
+    gcp -v $compile_template $WORKDIR/EXEC/
+    gcp -v $postp_code_path $WORKDIR/EXEC/postp/
+    cd EXEC/postp/
+    /home/fms/bin/mkmf -p cp_rvar -t $WORKDIR/EXEC/$compile_template:t $postp_code_path:t
+    gmake
+    if ($status != 0) then
+        echo compile of post-processor code failed, exiting...
+        exit
+    endif
+    cd $WORKDIR
+
+    gcp -v $rmvpr_code_path $WORKDIR/EXEC/rmvpr/
+    cd EXEC/rmvpr/
+    /home/fms/bin/mkmf -p rmvp -t $WORKDIR/EXEC/$compile_template:t $rmvpr_code_path:t
+    gmake
+    if ($status != 0) then
+        echo compile of remove-parallel-rivers code failed, exiting...
+        exit
+    endif
+    cd $WORKDIR
+
+    gcp -v $lakes_code_path $WORKDIR/EXEC/lakes/
+    cd EXEC/lakes/
+    /home/fms/bin/mkmf -p cplf -t $WORKDIR/EXEC/$compile_template:t -c "-Duse_netCDF" \
+          $lakes_code_path:t -c $fms_code_dir/shared/include /usr/local/include
+    gmake
+    if ($status != 0) then
+        echo compile of lakes code failed, exiting...
+        exit
+    endif
+    cd $WORKDIR
+else
+    gcp -v $home_dir/exec/postp/cp_rvar $WORKDIR/EXEC/postp/
+    gcp -v $home_dir/exec/rmvpr/rmvp    $WORKDIR/EXEC/rmvpr/
+    gcp -v $home_dir/exec/lakes/cplf    $WORKDIR/EXEC/lakes/
+    chmod ugo+x $WORKDIR/EXEC/postp/cp_rvar $WORKDIR/EXEC/rmvpr/rmvp $WORKDIR/EXEC/lakes/cplf
+endif
+
+dmget $river_network_ext_file
+
+
+# ------------------------------------------------
+#  RUN RIVER_REGRID
+# ------------------------------------------------
+echo ""
+echo RUN RIVER-REGRID
+$riv_regrd --mosaic $mosaic_file --river_src $river_network_ext_file --min_frac $min_frac --land_thresh $land_thresh
+mv river_output*nc OUTPUT/river_regrid/
+
+
+# ------------------------------------------------
+#  POST-PROCESS OUTPUT FROM RIVER_REGRID
+# ------------------------------------------------
+set river_input_files = OUTPUT/river_regrid/river_output*nc 
+echo $#river_input_files > fort.5
+foreach file ($river_input_files)
+   echo OUTPUT/river_regrid/$file:t >> fort.5
+end
+echo ""
+echo RUN POST-PROCESSOR
+$WORKDIR/EXEC/postp/cp_rvar < fort.5
+if ($status != 0) then
+    echo ERROR in post-processing river-regrid output, exiting...
+    exit
+else
+    mv river_network*nc OUTPUT/post_regrid/
+    mv out.cp_river_vars OUTPUT/post_regrid/
+endif
+
+
+# ------------------------------------------------
+#  REMOVE PARALLEL RIVERS
+# ------------------------------------------------
+set add_ocean_const = F
+
+set river_input_files = OUTPUT/post_regrid/river_network*nc
+echo $#river_input_files > fort.5
+foreach file ($river_input_files)
+   echo OUTPUT/post_regrid/$file:t >> fort.5
+end
+echo $add_ocean_const >> fort.5
+echo ""
+echo REMOVE PARALLEL RIVERS
+$WORKDIR/EXEC/rmvpr/rmvp < fort.5
+if ($status != 0) then
+    echo ERROR in removal of parallel rivers, exiting...
+    exit
+else
+    mv river_network*nc OUTPUT/rmv_parallel_rivers/
+    mv out.rmv_parallel_rivers OUTPUT/rmv_parallel_rivers/
+endif
+
+
+# ------------------------------------------------
+#  POST-PROCESS OUTPUT FROM REMOVE-PARALLEL-RIVERS
+# ------------------------------------------------
+set river_input_files = OUTPUT/rmv_parallel_rivers/river_network*nc
+echo $#river_input_files > fort.5
+foreach file ($river_input_files)
+   echo OUTPUT/rmv_parallel_rivers/$file:t >> fort.5
+end
+echo ""
+echo RUN POST-PROCESSOR
+$WORKDIR/EXEC/postp/cp_rvar < fort.5
+if ($status != 0) then
+    echo ERROR in post-processing parallel-river output, exiting...
+    exit
+else
+    mv river_network*nc OUTPUT/post_rmvp/
+    mv out.cp_river_vars OUTPUT/post_rmvp/
+endif
+
+
+# --------------------------------------
+#  ADD GLCC WATERBOD DATA
+# --------------------------------------
+set travel_thresh = 2.
+
+dmget $glcc_file
+gcp -v gfdl:$glcc_file   gfdl:$WORKDIR/DATA/
+set river_input_files = OUTPUT/post_rmvp/river_network*nc
+echo $#river_input_files > fort.5
+foreach file ($river_input_files)
+   echo OUTPUT/post_rmvp/$file:t >> fort.5
+end
+echo DATA/$glcc_file:t >> fort.5
+echo $travel_thresh >> fort.5
+touch input.nml
+echo ""
+echo ADD LAKES
+$WORKDIR/EXEC/lakes/cplf < fort.5
+if ($status != 0) then
+    echo ERROR in addition of lake data, exiting...
+    exit
+endif
+
+@ k = 0
+while ($k < $#river_input_files)
+  @ k = $k + 1
+  gcp OUTPUT/post_rmvp/river_network.tile"$k".nc $WORKDIR/hydrography.tile"$k".nc
+end
+set hydro_files = hydrography*.nc
+foreach file ($hydro_files)
+   set tn = `echo "$file" | awk -Ftile '{print $2}'`
+   ncks -A -a -v lake_frac,lake_depth_sill,lake_tau,WaterBod,PWetland,connected_to_next,whole_lake_area,max_slope_to_next \
+     lake_frac.tile"$tn" "$file:t"
+end
+
+if ($outdir != "") then
+    gcp -v hydrography.tile*nc $outdir/
+endif
+
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/FLOR2_grid-mercator.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/FLOR2_grid-mercator.py
new file mode 100644
index 0000000000..3a230ecd76
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/FLOR2_grid-mercator.py
@@ -0,0 +1,260 @@
+
+from scipy.interpolate import interp1d
+from midas.rectgrid import *
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap, cm
+import netCDF4 as nc
+
+
+# In[161]:
+
+lat0=-65.
+lon0=-300.
+lenlat=125.
+lenlon=360.
+nlon=360*2
+nlat=500*2
+mercator=supergrid(nlon,nlat,'mercator','degrees',lat0,lenlat,lon0,lenlon,cyclic_x=True)
+mercator.grid_metrics()
+#get_ipython().magic(u'matplotlib inline')
+phi=mercator.y[:,0]
+dphi=phi[1:]-phi[0:-1]
+plt.plot(phi[1:],2.0*dphi)
+plt.title('Mercator N/S resolution')
+
+
+# In[162]:
+
+#get_ipython().magic(u'matplotlib inline')
+phi=mercator.y[:,0]
+dphi=phi[1:]-phi[0:-1]
+jind=np.where(phi>-20.)[0][0]
+jind=jind+np.mod(jind,2)
+phi=phi[0:jind]
+dphi=dphi[0:jind]
+N=26
+phi_s = phi[-1]
+dphi_s = dphi[-1]
+phi_e = -10.
+dphi_e = 0.2
+nodes = [0,1,N-2,N-1]
+phi_nodes = [phi_s,phi_s+dphi_s,phi_e-dphi_e,phi_e]
+f2=interp1d(nodes,phi_nodes,kind=3)
+jInd2=np.arange(N)
+phi2=f2(jInd2)
+dphi2=phi2[1:]-phi2[0:-1]
+phi=np.concatenate((phi[0:-1],phi2))
+dphi=phi[1:]-phi[0:-1]
+
+N=56
+phi_s = phi[-1]
+phi2=np.linspace(phi_s,0,N)
+dphi2=phi2[1:]-phi2[0:-1]
+PHI=np.concatenate((phi[0:-1],phi2))
+PHI=np.concatenate((PHI[0:-1],-PHI[::-1]))
+
+DPHI=PHI[1:]-PHI[0:-1]
+plt.plot(PHI[1:],2*DPHI)
+
+
+# In[163]:
+
+LAMBDA=np.linspace(lon0,lon0+360.,nlon+1)
+jind=np.where(PHI>-78.)[0][0]
+jind=jind+np.mod(jind,2)
+jind2=np.where(PHI>60.)[0][0]
+jind2=jind2+np.mod(jind2,2)
+
+PHI2=PHI[jind:jind2-1]
+DPHI2=PHI2[1:]-PHI2[0:-1]
+x,y = np.meshgrid(LAMBDA,PHI2)
+gridB = supergrid(xdat=x,ydat=y,axis_units='degrees',cyclic_x=True)
+gridB.grid_metrics()
+gridB.write_nc('flor_gridB.nc')
+
+
+# In[164]:
+
+ny_ncap=155
+phi_tps = PHI2[-1]
+dphi_tps = DPHI2[-1]
+nodes = [0,1,ny_ncap-2,ny_ncap-1]
+phi_tp_nodes = [phi_tps,phi_tps+dphi_tps,89.9,90.]
+ftp=interp1d(nodes,phi_tp_nodes,kind=3)
+jInd2=np.arange(ny_ncap)
+phitp=ftp(jInd2)
+dphitp=phitp[1:]-phitp[0:-1]
+plt.plot(phitp[0:-1],2*dphitp)
+plt.title('Model N/S Arctic resolution (prior to transform)')
+
+
+# In[165]:
+
+lat0_tp = PHI2[-1]
+XTP,YTP=np.meshgrid(LAMBDA,phitp)
+gridC=supergrid(xdat=XTP,ydat=YTP,axis_units='degrees',tripolar_n=True)
+gridC.grid_metrics()
+gridC.write_nc('flor_gridC.nc')
+
+
+### Combine grids
+
+# In[166]:
+
+f=nc.Dataset('flor_gridB.nc')
+g=nc.Dataset('flor_gridC.nc')
+
+y1=f.variables['y'][:]
+y2=g.variables['y'][:]
+y=np.concatenate((y1,y2[1:,:]),axis=0)
+
+dy1=f.variables['dy'][:]
+dy2=g.variables['dy'][:]
+dy=np.concatenate((dy1,dy2),axis=0)
+
+x1=f.variables['x'][:]
+x2=g.variables['x'][:]
+x=np.concatenate((x1,x2[1:,:]),axis=0)
+
+dx1=f.variables['dx'][:]
+dx2=g.variables['dx'][:]
+dx=np.concatenate((dx1,dx2[1:,:]),axis=0)
+
+
+area1=f.variables['area'][:]
+area2=g.variables['area'][:]
+area=np.concatenate((area1,area2),axis=0)
+
+angle_dx1=f.variables['angle_dx'][:,:]
+angle_dx2=g.variables['angle_dx'][:-1,:]
+angle_dx=np.concatenate((angle_dx1,angle_dx2),axis=0)
+
+fout=nc.Dataset('supergrid.nc','w',format='NETCDF3_CLASSIC')
+
+ny=area.shape[0]; nx = area.shape[1]
+nyp=ny+1; nxp=nx+1
+
+fout.createDimension('nyp',nyp)
+fout.createDimension('nxp',nxp)
+fout.createDimension('ny',ny)
+fout.createDimension('nx',nx)
+string=fout.createDimension('string',255)    
+
+tile=fout.createVariable('tile','S1',('string'))
+yv=fout.createVariable('y','f8',('nyp','nxp'))
+xv=fout.createVariable('x','f8',('nyp','nxp'))    
+yv.units='degrees'
+xv.units='degrees'
+yv[:]=y
+xv[:]=x
+
+tile[0:4]='tile1'
+dyv=fout.createVariable('dy','f8',('ny','nxp'))
+dyv.units='meters'
+dyv[:]=dy
+dxv=fout.createVariable('dx','f8',('nyp','nx'))
+dxv.units='meters'
+dxv[:]=dx
+areav=fout.createVariable('area','f8',('ny','nx'))
+areav.units='m2'
+areav[:]=area
+anglev=fout.createVariable('angle_dx','f8',('nyp','nxp'))
+anglev.units='degrees'
+print nxp,nyp
+print angle_dx.shape
+anglev[:]=angle_dx            
+
+fout.sync()
+fout.close()
+
+
+### Load the final supergrid and create a model grid
+
+# In[167]:
+
+sgrid=supergrid(file='supergrid.nc',cyclic_x=True,tripolar_n=True)
+grid=quadmesh(supergrid=sgrid)
+print 'Total number of i grid points = ',grid.im
+print 'Total number of j grid points = ',grid.jm
+print 'Closest grid point to the equator =',grid.lath[np.where(np.abs(grid.lath)<0.2)[0]]
+
+
+# In[168]:
+
+aiso1=grid.dyh[:,0]/grid.dxh[:,0]
+aiso2=1.0/aiso1
+aiso = np.maximum(aiso1,aiso2)
+#plt.plot(grid.lath,aiso1,color='b')
+#plt.plot(grid.lath,aiso2,color='g')
+plt.plot(grid.lath,aiso,color='c')
+
+plt.plot([grid.lath[0],grid.lath[-1]],[1.,1.],'r--')
+plt.xlim(-80,50)
+plt.ylim(0.25,3)
+plt.title('Grid Anisotropy (nd)')
+
+
+# In[169]:
+
+#get_ipython().magic(u'matplotlib inline')
+print grid.x_T.min(),grid.x_T.max()
+plt.pcolormesh(grid.x_T,grid.y_T,grid.dxh*grid.dyh)
+plt.xlim(-300,60)
+plt.ylim(-90,90)
+plt.title('Cell Area')
+plt.colorbar()
+
+
+# In[172]:
+
+wd=6667000.
+ht=6667000.
+m =  Basemap(projection='stere',width=wd,height=ht,lon_0=0.0,lat_ts=70.,lat_0=90.,resolution='l')
+xx=grid.x_T.copy()
+yy=grid.y_T.copy()
+xx[xx<-180]=xx[xx<-180]+360.
+xplt,yplt=m(xx,yy,inverse=False)
+aiso = grid.dyh/grid.dxh
+aiso = np.maximum(aiso,1.0/aiso)
+cf=m.contourf(xplt,yplt,aiso,np.linspace(1.0,3,10))
+#m.contourf(xplt,yplt,aiso)
+#cl=m.contour(xplt,yplt,aiso,np.linspace(1.0,2.0,10),colors='k')
+#plt.clabel(cl)
+m.drawcoastlines()
+plt.title('Tripolar Anisotropy')
+plt.colorbar(cf)
+
+
+# In[171]:
+
+wd=6667000.
+ht=6667000.
+m =  Basemap(projection='stere',width=wd,height=ht,lon_0=0.0,lat_ts=70.,lat_0=90.,resolution='l')
+xx=grid.x_T
+yy=grid.y_T
+xx[xx<-180]=xx[xx<-180]+360.
+xplt,yplt=m(xx,yy,inverse=False)
+cf=m.contourf(xplt,yplt,grid.dxh*grid.dyh/1.e9,np.arange(0.1,3.6,.1),extend='both')
+m.drawcoastlines()
+m.drawparallels(np.arange(50.,90.,10.))
+m.drawmeridians(np.arange(-180.,180.,15.))
+plt.title('Tripolar grid cell Area')
+plt.colorbar(cf)
+
+
+# In[ ]:
+
+ingrid=quadmesh('GEBCO_08_v2.nc',var='depth',simple_grid=True,cyclic=True)
+TOPO=state('GEBCO_08_v2.nc',grid=ingrid,fields=['depth'])
+TOPO.rename_field('depth','topo')
+TOPO.var_dict['topo']['Ztype']='Fixed'
+fnam = 'topog_gebco.nc'
+R=TOPO.subtile('topo',target=grid)
+R.write_nc(fnam,['mean','max','min','std','count'])
+
+
+# In[ ]:
+
+
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/Jenkins.csh b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/Jenkins.csh
new file mode 100755
index 0000000000..69cc4f9ef2
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/Jenkins.csh
@@ -0,0 +1,21 @@
+#!/bin/csh -fx
+
+echo -n Started $0 in ; pwd
+
+# Modules
+source $MODULESHOME/init/csh
+module use -a /home/fms/local/modulefiles
+module load python
+module load netcdf/4.2 intel_compilers
+module load nco/4.3.1
+
+# Installing this file avoids a slow and unreliable file server (that would otherwise be ftp'd from ftp.oce.orst.edu)
+cp -n /archive/gold/datasets/obs/tpxo7_atlas_netcdf.tar.Z .
+# This is the one non-reproducible file needed in this work flow
+cp -n /archive/aja/permanent/edit_topog.nc.v20140629 edit_topog.nc
+
+# Work around for environment problem inside MIDAS
+setenv PYTHONPATH $cwd/local/lib
+
+# Run through the work flow
+make
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/Makefile b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/Makefile
new file mode 100644
index 0000000000..993b9a3341
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/Makefile
@@ -0,0 +1,130 @@
+# Makefile to create supergrid.nc and interpolated_topog.nc
+# To use:
+#   module load python
+#   setenv PYTHONPATH $cwd/MIDAS
+#
+# then
+#   make all
+# or
+#   make supergrid.nc
+#   make interpolated_topog.nc
+
+SHELL=tcsh -f
+
+all: input forcing obs mosaic
+	md5sum -c md5sums.txt
+
+showenv:
+	env
+	-set
+	-module list
+	which python
+	-python --version
+
+input: ocean_hgrid.nc topog.nc basin_codes.nc
+
+forcing: salt_restore.nc seawifs_1998-2006_smoothed_2X.nc tidal_amplitude.nc precip_factor.nc
+
+obs: WOA05_ptemp_salt_monthly.nc WOA05_ptemp_salt_annual.nc
+
+mosaic: C180_mosaic/mosaic.tar 
+
+# Grid
+
+ocean_hgrid.nc: local
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python create_grids.py
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python merge_grids.py
+	./changeChar.py ocean_hgrid.nc tile tile1
+
+# Topography
+interpolated_topog.nc: GEBCO_08_v2.nc ocean_hgrid.nc
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python create_topo.py
+
+
+edit_topog.nc: interpolated_topog.nc
+	./applyTopoEdits.py --orig=interpolated_topog.nc --edits=topo_edits --output=edit_topog.nc
+
+topog.nc: edit_topog.nc interpolated_topog.nc
+	./ice9.py edit_topog.nc --output topog.nc
+
+basin_codes.nc: topog.nc
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python make_basin_mask.py
+	ncatted -h -a flag_meanings,basin,c,c,'1:Southern Ocean, 2:Atlantic Ocean, 3:Pacific Ocean, 4:Arctic Ocean, 5:Indian Ocean, 6:Mediterranean Sea, 7:Black Sea, 8:Hudson Bay, 9:Baltic Sea, 10:Red Sea, 11:Persian Gulf' basin_codes.nc
+	ncatted -h -a flag_values,basin,c,c,'1,2,3,4,5,6,7,8,9,10,11' basin_codes.nc
+
+WOA05_ptemp_salt_annual.nc: WOA05_ptemp_salt_monthly.nc
+	ncra -h -O $< $@
+
+WOA05_ptemp_salt_monthly.nc: ocean_hgrid.nc topog.nc /archive/gold/datasets/obs/WOA05_pottemp_salt.nc interpWOA05.py local
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python interpWOA05.py
+	ncatted -h -a modulo,TIME,c,c,' ' WOA05_ptemp_salt_monthly.nc
+
+salt_restore.nc: ocean_hgrid.nc topog.nc /archive/gold/datasets/obs/WOA05_pottemp_salt.nc interpSaltRestore.py local
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python interpSaltRestore.py
+	ncatted -h -a modulo,TIME,c,c,' ' salt_restore.nc
+
+precip_factor.nc: ocean_hgrid.nc  local
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python create_precip_scale2d.py
+
+seawifs_1998-2006_smoothed_2X.nc: ocean_hgrid.nc topog.nc /archive/gold/datasets/global/siena_201204/INPUT/seawifs_1998-2006_GOLD_smoothed_2X.nc interpCHL.py local
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python interpCHL.py
+	ncatted -h -a modulo,TIME,c,c,' ' seawifs_1998-2006_smoothed_2X.nc
+
+tidal_amplitude.nc: DATA interpTides.py local
+	unlimit stacksize; setenv PYTHONPATH ./local/lib/python; python interpTides.py
+
+C180_mosaic/mosaic.tar: topog.nc
+	(cd C180_mosaic;./create_c180_mosaic_and_rivers.csh)
+
+
+DATA: tpxo7_atlas_netcdf.tar.Z GEBCO_08_v2.nc
+	tar xzvf tpxo7_atlas_netcdf.tar.Z
+	touch DATA
+
+tpxo7_atlas_netcdf.tar.Z:
+#	wget ftp://ftp.oce.orst.edu/dist/tides/Global/tpxo7_atlas_netcdf.tar.Z
+	ln -s  /archive/gold/datasets/obs/tpxo7_atlas_netcdf.tar.Z
+
+MIDAS:
+	git clone https://github.com/mjharriso/MIDAS.git
+
+# Submodule method for obtaining MIDAS (used by frepp_local target) 
+MIDAS/README.md:
+	(cd ../../..; git submodule init tools/python/025gridGeneration/MIDAS)
+	(cd ../../..; git submodule update tools/python/025gridGeneration/MIDAS)
+
+local: MIDAS 
+	-rm -rf $</build/*
+	mkdir -p $@
+	cd $<; make -f Makefile_GFDL INSTALL_PATH=../local
+	touch $@
+
+# This is a target to allow frepp analysis scripts to build MIDAS only once and then clean up
+frepp_local: MIDAS/README.md
+	make local
+	rm -rf MIDAS
+	mkdir MIDAS
+	touch $<
+	touch local
+
+GEBCO_08_v2.nc :
+	ln -s /archive/gold/datasets/topography/$@ .
+
+md5sums.txt: ocean_hgrid.nc  interpolated_topog.nc topog.nc basin_codes.nc salt_restore.nc seawifs_1998-2006_smoothed_2X.nc tidal_amplitude.nc WOA05_ptemp_salt_monthly.nc WOA05_ptemp_salt_annual.nc C180_mosaic/mosaic.tar C180_mosaic/quick_mosaic.tar
+	echo Grids > $@
+	md5sum sgrid*.nc ocean_hgrid.nc >> $@
+	echo >> $@
+	echo Topography >> $@
+	md5sum *topog*.nc >> $@
+	md5sum basin_codes.nc >> $@
+	echo >> $@
+	echo Data >> $@
+	md5sum salt_restore.nc seawifs_1998-2006_smoothed_2X.nc tidal_amplitude.nc >> $@
+	echo >> $@
+	echo Obs >> $@
+	md5sum WOA05_ptemp_salt_{monthly,annual}.nc >> $@
+	echo >> $@
+	echo Mosaic >> $@
+	md5sum C180_mosaic/mosaic.tar >> $@
+	echo Quick Mosaic >> $@
+	md5sum C180_mosaic/quick_mosaic.tar >> $@
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/README.md b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/README.md
new file mode 100644
index 0000000000..42fee42d9d
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/README.md
@@ -0,0 +1,10 @@
+The Makefile in this directory manages the complete work flow to create
+the grids, topography and some oceanographic fields needed to run OM4_025.
+
+The environment must be setup outside of make. At GFDL do:
+
+  source $MODULESHOME/init/csh
+  module use -a /home/fms/local/modulefiles
+  module load netcdf/4.2 intel_compilers
+  module load nco/4.3.1
+  module load python
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/addDimension.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/addDimension.py
new file mode 100755
index 0000000000..7b2bc8a20f
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/addDimension.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+def error(msg,code=9):
+  print 'Error: ' + msg
+  exit(code)
+
+# Imports
+try: import argparse
+except: error('This version of python is not new enough. python 2.7 or newer is required.')
+try: from netCDF4 import Dataset
+except: error('Unable to import netCDF4 module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_netcdf4')
+
+def main():
+
+  # Command line arguments
+  parser = argparse.ArgumentParser(description=
+       'Adds the named dimension to a netcdf file.',
+       epilog='Written by A.Adcroft, 2013.')
+  parser.add_argument('filename', type=str,
+                      help='netcdf file to modify.')
+  parser.add_argument('dimension', type=str,
+                      nargs='?', default='ntiles',
+                      help='Name of dimension to add.')
+  parser.add_argument('value', type=int,
+                      nargs='?', default=1,
+                      help='Value to set dimension length to.')
+
+  optCmdLineArgs = parser.parse_args()
+
+  rg = Dataset(optCmdLineArgs.filename, 'a' );
+  rg.createDimension(optCmdLineArgs.dimension,optCmdLineArgs.value)
+  rg.close()
+  print 'File "%s" updated.'%(optCmdLineArgs.filename)
+
+# Invoke main()
+if __name__ == '__main__': main()
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/applyTopoEdits.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/applyTopoEdits.py
new file mode 100755
index 0000000000..db2cb062c4
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/applyTopoEdits.py
@@ -0,0 +1,119 @@
+#!/usr/bin/env python
+
+import pandas
+import netCDF4 as nc
+import argparse
+
+"""
+Read a list of edits which were provided by and apply them to a topo file.
+"""
+
+parser = argparse.ArgumentParser()
+parser.add_argument('--orig',type=str,help='path to original netCDF file',default=None)
+parser.add_argument('--edits',type=str,help='path to edits file',default=None)
+parser.add_argument('--output',type=str,help='path to output netCDF file',default='topog_edits.nc')
+
+
+
+
+args=parser.parse_args()
+
+
+if args.orig is None:
+
+    print "Need to supply the original topography netCDF file"
+    raise
+
+if args.edits is None:
+
+    print "Need to supply the path to edits "
+    raise
+
+f=nc.Dataset(args.orig)
+depth=f.variables['depth'][:]
+Edits=pandas.read_csv(args.edits,header=None,names=['jEdit','iEdit','New'])
+jEdits=[]
+for a in Edits.jEdit:
+    jEdits.append(int(a.replace('(','')))
+
+iEdits=[]
+for a in Edits.iEdit:
+    iEdits.append(int(a))
+
+New=[]
+for a in Edits.New:
+    New.append(float(a.replace(')','')))
+        
+Orig=[]
+for i,j,new in zip(iEdits,jEdits,New):
+    Orig.append(depth[j,i])
+    print i,j,depth[j,i],new
+
+g=nc.Dataset(args.output,'w',format='NETCDF3_CLASSIC')
+
+dims=[]
+for d in f.dimensions:
+    dimvals=f.variables[d][:]
+    dims.append((d,len(f.dimensions[d]),dimvals))
+
+
+for d in dims:
+    g.createDimension(d[0],d[1])
+    
+g.createDimension('nEdits',None)
+
+dimv=[]
+for d in dims:
+    dimv.append(g.createVariable(d[0],'f8',(d[0])))
+    dimv[-1].units='degrees'
+
+for v,d in zip(dimv,dims):
+    v[:]=d[2]
+    
+for v in f.variables:
+    if v.find('longitude')==-1 and v.find('latitude')==-1 and v.find('count')==-1 and v.find('min')==-1 and v.find('max')==-1:
+        varout=g.createVariable(v,'f4',('latitude','longitude'))
+        try:
+            units=f.variables[v].units
+            varout.units=units
+        except:
+            pass
+        try:
+            standard_name=f.variables[v].standard_name
+            varout.standard_name=standard_name
+        except:
+            pass
+        try:
+            description=f.variables[v].description
+            varout.description=description
+        except:
+            pass
+        try:
+            long_name=f.variables[v].long_name
+            varout.long_name=long_name
+        except:
+            pass
+
+        dat=f.variables[v][:]
+
+        if v.find('depth')>-1:
+            for i,j,d in zip(iEdits,jEdits,New):
+                print 'Modifying Elevation at j,i= ',j,i,' Old= ',dat[j,i],' New= ',-d
+                dat[j,i]=-d
+
+        varout[:]=dat
+
+        
+ivar=g.createVariable('iEdit','i4',('nEdits'))
+jvar=g.createVariable('jEdit','i4',('nEdits'))
+kvar=g.createVariable('zEdit','f4',('nEdits'))
+
+n=0
+for i,j,d in zip(iEdits,jEdits,Orig):
+    ivar[n]=i
+    jvar[n]=j
+    kvar[n]=d
+    n=n+1
+
+g.sync()
+g.close()
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/changeChar.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/changeChar.py
new file mode 100755
index 0000000000..a8c99850dd
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/changeChar.py
@@ -0,0 +1,45 @@
+#!/usr/bin/env python
+
+def error(msg,code=9):
+  print 'Error: ' + msg
+  exit(code)
+
+# Imports
+try: import argparse
+except: error('This version of python is not new enough. python 2.7 or newer is required.')
+try: from netCDF4 import Dataset, stringtochar
+except: error('Unable to import netCDF4 module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_netcdf4')
+try: import numpy as np
+except: error('Unable to import numpy module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_numpy')
+
+def main():
+
+  # Command line arguments
+  parser = argparse.ArgumentParser(description=
+       'Changes the value of a named char variable in a netcdf file.',
+       epilog='Written by A.Adcroft, 2013.')
+  parser.add_argument('filename', type=str,
+                      help='netcdf file to modify.')
+  parser.add_argument('variable', type=str,
+                      help='Name of char variable to change.')
+  parser.add_argument('value', type=str,
+                      help='Contents to change string to.')
+
+  optCmdLineArgs = parser.parse_args()
+
+  rg = Dataset(optCmdLineArgs.filename, 'a' );
+  if optCmdLineArgs.variable in rg.variables:
+    var = rg.variables[optCmdLineArgs.variable]
+    dat = np.empty(1,'S'+repr(len(var)))
+    dat[0] = optCmdLineArgs.value
+    dc = stringtochar(dat)
+    var[:] = dc
+  else: error('"'+optCmdLineArgs.variable+'" was not found in "'+optCmdLineArgs.filename+'".')
+  rg.close()
+  print 'File "%s" updated.'%(optCmdLineArgs.filename)
+
+# Invoke main()
+if __name__ == '__main__': main()
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_grids.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_grids.py
new file mode 100755
index 0000000000..0bb7fb8c31
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_grids.py
@@ -0,0 +1,98 @@
+#!/usr/bin/env python
+
+#============================================================
+# Generate tiles for the northern/southern caps
+# and central mercator grid. For use in Antarctic ice-sheet
+# modeling.
+#
+# 
+# python create_topo.py
+# Output: mercator_supergrid.nc, ncap_supergrid.nc, scap_supergrid.nc
+# These are supergrids (2x grid tracer refinement) containing positions
+# cell lengths, areas and angles.  
+#
+#============================================================
+
+
+
+from midas.rectgrid_gen import *
+from midas.rectgrid import *
+import numpy as np
+from scipy.interpolate import interp1d
+lat0=-65.
+lon0=-300.
+lenlat=125.
+lenlon=360.
+nlon=360*2
+nlat=500*2
+
+mercator=supergrid(nlon,nlat,'mercator','degrees',lat0,lenlat,lon0,lenlon,cyclic_x=True)
+mercator.grid_metrics()
+
+phi=mercator.y[:,0]
+dphi=phi[1:]-phi[0:-1]
+
+phi=mercator.y[:,0]
+dphi=phi[1:]-phi[0:-1]
+jind=np.where(phi>-20.)[0][0]
+jind=jind+np.mod(jind,2)
+phi=phi[0:jind]
+dphi=dphi[0:jind]
+
+N=26
+phi_s = phi[-1]
+dphi_s = dphi[-1]
+phi_e = -10.
+dphi_e = 0.2
+nodes = [0,1,N-2,N-1]
+phi_nodes = [phi_s,phi_s+dphi_s,phi_e-dphi_e,phi_e]
+f2=interp1d(nodes,phi_nodes,kind=3)
+jInd2=np.arange(N)
+phi2=f2(jInd2)
+dphi2=phi2[1:]-phi2[0:-1]
+phi=np.concatenate((phi[0:-1],phi2))
+dphi=phi[1:]-phi[0:-1]
+
+N=56
+phi_s = phi[-1]
+phi2=np.linspace(phi_s,0,N)
+dphi2=phi2[1:]-phi2[0:-1]
+PHI=np.concatenate((phi[0:-1],phi2))
+PHI=np.concatenate((PHI[0:-1],-PHI[::-1]))
+
+DPHI=PHI[1:]-PHI[0:-1]
+
+LAMBDA=np.linspace(lon0,lon0+360.,nlon+1)
+jind=np.where(PHI>-78.)[0][0]
+jind=jind+np.mod(jind,2)
+jind2=np.where(PHI>60.)[0][0]
+jind2=jind2+np.mod(jind2,2)
+
+PHI2=PHI[jind:jind2-1]
+DPHI2=PHI2[1:]-PHI2[0:-1]
+x,y = np.meshgrid(LAMBDA,PHI2)
+gridA = supergrid(xdat=x,ydat=y,axis_units='degrees',cyclic_x=True)
+gridA.grid_metrics()
+gridA.write_nc('sgridA.nc')
+
+
+ny_ncap=155
+phi_tps = PHI2[-1]
+dphi_tps = DPHI2[-1]
+nodes = [0,1,ny_ncap-2,ny_ncap-1]
+phi_tp_nodes = [phi_tps,phi_tps+dphi_tps,89.9,90.]
+ftp=interp1d(nodes,phi_tp_nodes,kind=3)
+jInd2=np.arange(ny_ncap)
+phitp=ftp(jInd2)
+dphitp=phitp[1:]-phitp[0:-1]
+
+lat0_tp = PHI2[-1]
+XTP,YTP=np.meshgrid(LAMBDA,phitp)
+gridB=supergrid(xdat=XTP,ydat=YTP,axis_units='degrees',tripolar_n=True)
+gridB.grid_metrics()
+gridB.write_nc('sgridB.nc')
+
+
+
+
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_precip_scale2d.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_precip_scale2d.py
new file mode 100644
index 0000000000..4349b9b414
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_precip_scale2d.py
@@ -0,0 +1,38 @@
+# IPython log file
+
+from midas.rectgrid import *
+import matplotlib.pyplot as plt
+import netCDF4 as nc
+import numpy as np
+
+f=nc.Dataset('precip_factor.nc','w',format='NETCDF3_CLASSIC')
+x=np.linspace(0.,360.,721)
+y=np.linspace(-90.,90.,361)
+X,Y=np.meshgrid(x,y)
+sgrid=supergrid(xdat=X, ydat=Y, cyclic_x=True)
+grid=quadmesh(supergrid=sgrid)
+f.createDimension('lon',grid.im)
+f.createDimension('lat',grid.jm)
+timed=f.createDimension('time',0)
+timev=f.createVariable('time','f8',('time'))
+lonv=f.createVariable('lon','f8',('lon'))
+latv=f.createVariable('lat','f8',('lat'))
+timev.calendar='julian'
+timev.modulo='yes'
+timev.units='days since 1900-01-01 00:00:00'
+timev.cartesian_axis='T'
+lonv.units='degrees_E'
+lonv.modulo='yes'
+lonv.cartesian_axis='X'
+latv.units='degrees_N'
+latv.cartesian_axis='Y'
+lonv[:]=grid.lonh
+latv[:]=grid.lath
+precip_scale = np.ones(grid.wet.shape)
+precip_scale[grid.y_T>40]=np.cos((grid.y_T[grid.y_T>40.]-40.)*np.pi/180.)
+v=f.createVariable('precip_factor','f8',('time','lat','lon'))
+v[0,:]=precip_scale
+timev[0]=0.
+f.sync()
+f.close()
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_topo.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_topo.py
new file mode 100755
index 0000000000..8c891a32a6
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/create_topo.py
@@ -0,0 +1,43 @@
+#!python
+
+#============================================================
+# Generate tiles for the northern/southern caps
+# and central mercator grid.
+#
+# python create_topo.py
+# Output: mercator_supergrid.nc, ncap_supergrid.nc, scap_supergrid.nc
+# These are supergrids (2x grid tracer refinement) containing positions
+# cell lengths, areas and angles
+#
+# Generate topography for grid tiles using BEDMAP for the Antarctic cap
+# GEBCO 2 minute data for the Mercator grid and either
+# IBCAO or GEBCO for the Northern cap (these files need to be linked to the
+# current directory prior to running this command)
+
+# python create_topo.py --tile ncap|scap|mercator 
+#
+#============================================================
+
+
+
+from midas.rectgrid import *
+from midas.rectgrid_gen import *
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap, cm
+from mpl_toolkits.basemap import interp 
+import argparse
+
+sgrid=supergrid(file='ocean_hgrid.nc',cyclic_x=True,tripolar_n=True)
+grid=quadmesh(supergrid=sgrid)
+
+ingrid=quadmesh('GEBCO_08_v2.nc',var='depth',simple_grid=True,cyclic=True)
+TOPO=state('GEBCO_08_v2.nc',grid=ingrid,fields=['depth'])
+TOPO.rename_field('depth','topo')
+TOPO.var_dict['topo']['Ztype']='Fixed'
+
+fnam='interpolated_topog.nc'   
+R=TOPO.subtile('topo',target=grid)
+R.rename_field('mean','depth')
+R.write_nc(fnam,['depth','max','min','std','count'])
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/editTopo.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/editTopo.py
new file mode 100755
index 0000000000..0478bb4671
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/editTopo.py
@@ -0,0 +1,456 @@
+#!/usr/bin/env python
+
+def error(msg,code=9):
+  print 'Error: ' + msg
+  exit(code)
+
+
+# Imports
+try: import argparse
+except: error('This version of python is not new enough. python 2.7 or newer is required.')
+try: from netCDF4 import Dataset
+except: error('Unable to import netCDF4 module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_netcdf4')
+try: import numpy as np
+except: error('Unable to import numpy module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_numpy')
+try: import matplotlib.pyplot as plt
+except: error('Unable to import matplotlib.pyplot module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_matplotlib')
+from matplotlib.widgets import Button, RadioButtons
+from matplotlib.colors import LinearSegmentedColormap
+import shutil as sh
+
+
+def main():
+
+  # Command line arguments
+  parser = argparse.ArgumentParser(description=
+       '''Point-wise editting of topography.
+          Button 1 assigns the prescribed level to the cell at the mouse pointer.
+          Adjust the prescribed value with buttons on the bottom.
+          Double click button 1 assigns the highest of the nearest ocean points.
+          Right click on a cell resets to the original value.
+          Scroll wheel zooms in and out.
+          Move the "data window" around with the North, South, East and West buttons.
+          Closing the window writes the file to the output file if one is specified with --output.
+        ''',
+       epilog='Written by A.Adcroft, 2013.')
+  parser.add_argument('filename', type=str,
+                      help='netcdf file to read.')
+  parser.add_argument('variable', type=str,
+                      nargs='?', default='depth',
+                      help='Name of variable to edit. Defaults to "depth".')
+  parser.add_argument('--output', type=str,
+                      nargs='?', default=' ',
+                      help='Write an output file. If no output file is specified, creates the file with the "edit_" prepended to the name  of the input file.')
+
+  optCmdLineArgs = parser.parse_args()
+
+  createGUI(optCmdLineArgs.filename, optCmdLineArgs.variable, optCmdLineArgs.output)
+
+
+def createGUI(fileName, variable, outFile):
+
+  # Open netcdf file
+  try: rg=Dataset( fileName, 'r' );
+  except: error('There was a problem opening "'+fileName+'".')
+
+  rgVar = rg.variables[variable] # handle to the variable
+  dims = rgVar.dimensions # tuple of dimensions
+  depth = rgVar[:] # Read the data
+  #depth = depth[0:600,0:600]
+  (nj,ni) = depth.shape
+  print 'Range of input depths: min=',np.amin(depth),'max=',np.amax(depth)
+
+  try:
+    sg=Dataset( 'supergrid.nc', 'r' );
+    lon = sg.variables['x'][:]; lon = lon[0:2*nj+1:2,0:2*ni+1:2]
+    lat = sg.variables['y'][:]; lat = lat[0:2*nj+1:2,0:2*ni+1:2]
+  except:
+    lon, lat = np.meshgrid( np.arange(ni+1), np.arange(nj+1) )
+  fullData = Topography(lon, lat, depth)
+
+  class Container:
+    def __init__(self):
+      self.view = None
+      self.edits = None
+      self.data = None
+      self.quadMesh = None
+      self.ax = None
+      self.syms = None
+      cdict = {'red': ((0.0, 0.0, 0.0), (0.5, 0.7, 0.0), (1.0, 0.9, 0.0)),
+           'green': ((0.0, 0.0, 0.0), (0.5, 0.7, 0.2), (1.0, 1.0, 0.0)),
+            'blue': ((0.0, 0.0, 0.2), (0.5, 1.0, 0.0), (1.0, 0.9, 0.0))}
+      self.cmap = LinearSegmentedColormap('my_colormap',cdict,256)
+      self.clim = 6000
+      self.climLabel = None
+  All = Container()
+  All.view = View(ni,nj)
+  All.edits = Edits()
+
+  # Read edit data, if it exists
+  if 'iEdit' in rg.variables:
+    jEdit = rg.variables['iEdit'][:]; iEdit = rg.variables['jEdit'][:]
+    zEdit = rg.variables['zEdit'][:]
+    for l,i in enumerate(iEdit):
+      All.edits.setVal( fullData.height[iEdit[l],jEdit[l]] )
+      fullData.height[iEdit[l],jEdit[l]] = zEdit[l] # Restore data
+      All.edits.add( iEdit[l],jEdit[l] )
+  All.data = fullData.cloneWindow( (All.view.i0,All.view.j0), (All.view.iw,All.view.jw) )
+  if All.edits.ijz: All.data.applyEdits(fullData, All.edits.ijz)
+
+  # A mask based solely on value of depth
+  #notLand = np.where( depth<0, 1, 0)
+  #wet = ice9it(600,270,depth)
+
+  All.quadMesh = plt.pcolormesh(All.data.longitude,All.data.latitude,All.data.height,cmap=All.cmap,vmin=-All.clim,vmax=All.clim)
+  All.syms = All.edits.plot(fullData)
+  dir(All.syms)
+  All.ax=plt.gca(); All.ax.set_xlim( All.data.xlim ); All.ax.set_ylim( All.data.ylim )
+  All.climLabel = plt.figtext(.97,.97, 'XXXXX', ha='right', va='top')
+  All.climLabel.set_text('clim = $\pm$%i'%(All.clim))
+  All.edits.label = plt.figtext(.97,.03, 'XXXXX', ha='right', va='bottom')
+  All.edits.label.set_text('New depth = %i'%(All.edits.get()))
+  lowerButtons = Buttons()
+  def resetDto0(event): All.edits.setVal(0)
+  lowerButtons.add('Set 0', resetDto0)
+  def resetDto100(event): All.edits.addToVal(100)
+  lowerButtons.add('+100', resetDto100)
+  def resetDto100(event): All.edits.addToVal(30)
+  lowerButtons.add('+30', resetDto100)
+  def resetDto100(event): All.edits.addToVal(10)
+  lowerButtons.add('+10', resetDto100)
+  def resetDto100(event): All.edits.addToVal(3)
+  lowerButtons.add('+3', resetDto100)
+  def resetDto100(event): All.edits.addToVal(1)
+  lowerButtons.add('+1', resetDto100)
+  def resetDto100(event): All.edits.addToVal(-1)
+  lowerButtons.add('-1', resetDto100)
+  def resetDto100(event): All.edits.addToVal(-3)
+  lowerButtons.add('-3', resetDto100)
+  def resetDto100(event): All.edits.addToVal(-10)
+  lowerButtons.add('-10', resetDto100)
+  def resetDto100(event): All.edits.addToVal(-30)
+  lowerButtons.add('-30', resetDto100)
+  def resetDto100(event): All.edits.addToVal(-100)
+  lowerButtons.add('-100', resetDto100)
+  def resetDto100(event): All.edits.addToVal(-500)
+  lowerButtons.add('-500', resetDto100)
+  def undoLast(event):
+    All.edits.pop()
+    All.data = fullData.cloneWindow( (All.view.i0,All.view.j0), (All.view.iw,All.view.jw) )
+    All.data.applyEdits(fullData, All.edits.ijz)
+    All.quadMesh.set_array(All.data.height.ravel())
+    All.edits.updatePlot(fullData,All.syms)
+    plt.draw()
+  lowerButtons.add('Undo', undoLast)
+  upperButtons = Buttons(bottom=1-.0615)
+  def colorScale(event):
+    Levs = [50, 200, 1000, 6000]
+    i = Levs.index(All.clim)
+    if event=='+clim': i = min(i+1, len(Levs)-1)
+    elif event==' -clim': i = max(i-1, 0)
+    All.clim = Levs[i]
+    #All.quadMesh = plt.pcolormesh(All.data.longitude,All.data.latitude,All.data.height,cmap=All.cmap,vmin=-All.clim,vmax=All.clim)
+    #All.ax.set_xlim( All.data.xlim ); All.ax.set_ylim( All.data.ylim )
+    All.quadMesh.set_clim(vmin=-All.clim, vmax=All.clim)
+    All.climLabel.set_text('clim = $\pm$%i'%(All.clim))
+    plt.draw()
+  def moveVisData(di,dj):
+    All.view.move(di,dj)
+    All.data = fullData.cloneWindow( (All.view.i0,All.view.j0), (All.view.iw,All.view.jw) )
+    All.data.applyEdits(fullData, All.edits.ijz)
+    plt.sca(All.ax); plt.cla()
+    All.quadMesh = plt.pcolormesh(All.data.longitude,All.data.latitude,All.data.height,cmap=All.cmap,vmin=-All.clim,vmax=All.clim)
+    All.ax.set_xlim( All.data.xlim ); All.ax.set_ylim( All.data.ylim )
+    All.syms = All.edits.plot(fullData)
+    plt.draw()
+  def moveWindowLeft(event): moveVisData(-1,0)
+  upperButtons.add('West', moveWindowLeft)
+  def moveWindowRight(event): moveVisData(1,0)
+  upperButtons.add('East', moveWindowRight);
+  def moveWindowDown(event): moveVisData(0,-1)
+  upperButtons.add('South', moveWindowDown)
+  def moveWindowUp(event): moveVisData(0,1)
+  upperButtons.add('North', moveWindowUp)
+  climButtons = Buttons(bottom=1-.0615,left=0.65)
+  def incrCScale(event): colorScale('+clim')
+  climButtons.add('Incr', incrCScale)
+  def incrCScale(event): colorScale(' -clim')
+  climButtons.add('Decr', incrCScale)
+  plt.sca(All.ax)
+  def onClick(event): # Mouse button click
+    if event.inaxes==All.ax and event.button==1 and event.xdata:
+      (i,j) = findPointInMesh(fullData.longitude, fullData.latitude, event.xdata, event.ydata)
+      if not i==None:
+        (I,J) = findPointInMesh(All.data.longitude, All.data.latitude, event.xdata, event.ydata)
+        if event.dblclick:
+          nVal = -99999
+          if All.data.height[I+1,J]<0: nVal = max(nVal, All.data.height[I+1,J])
+          if All.data.height[I-1,J]<0: nVal = max(nVal, All.data.height[I-1,J])
+          if All.data.height[I,J+1]<0: nVal = max(nVal, All.data.height[I,J+1])
+          if All.data.height[I,J-1]<0: nVal = max(nVal, All.data.height[I,J-1])
+          if nVal==-99999: return
+          All.edits.add(i,j,nVal)
+          All.data.height[I,J] = nVal
+        else:
+          All.edits.add(i,j)
+          All.data.height[I,J] = All.edits.get()
+        All.quadMesh.set_array(All.data.height.ravel())
+        All.edits.updatePlot(fullData,All.syms)
+        plt.draw()
+    elif event.inaxes==All.ax and event.button==3 and event.xdata:
+      (i,j) = findPointInMesh(fullData.longitude, fullData.latitude, event.xdata, event.ydata)
+      if not i==None:
+        All.edits.delete(i,j)
+        All.data = fullData.cloneWindow( (All.view.i0,All.view.j0), (All.view.iw,All.view.jw) )
+        All.data.applyEdits(fullData, All.edits.ijz)
+        All.quadMesh.set_array(All.data.height.ravel())
+        All.edits.updatePlot(fullData,All.syms)
+        plt.draw()
+    elif event.inaxes==All.ax and event.button==2 and event.xdata: zoom(event) # Re-center
+  plt.gcf().canvas.mpl_connect('button_press_event', onClick)
+  def zoom(event): # Scroll wheel up/down
+    if event.button == 'up': scale_factor = 1/1.5 # deal with zoom in
+    elif event.button == 'down': scale_factor = 1.5 # deal with zoom out
+    else: scale_factor = 1.0
+    new_xlim, new_ylim = newLims( \
+           All.ax.get_xlim(), All.ax.get_ylim(), (event.xdata,event.ydata), \
+           All.data.xlim, All.data.ylim, scale_factor)
+    if not new_xlim: return # No changein limits
+    All.ax.set_xlim(new_xlim[0], new_xlim[1]); All.ax.set_ylim(new_ylim[0], new_ylim[1])
+    plt.draw() # force re-draw
+  plt.gcf().canvas.mpl_connect('scroll_event', zoom)
+  def statusMesg(x,y):
+    j,i = findPointInMesh(fullData.longitude, fullData.latitude, x, y)
+    if not i==None: return 'lon,lat=%.2f,%.2f  depth(%i,%i)=%.2f'%(x,y,i,j,fullData.height[j,i])
+    else: return 'lon,lat=%.3f,%.3f'%(x,y)
+  All.ax.format_coord = statusMesg
+  plt.show()
+  All.edits.list()
+  if not outFile: outFile = 'edit_'+fileName
+  if not outFile==' ':
+    print 'Creating new file "'+outFile+'"'
+    # Create new netcdf file
+    if not fileName==outFile: sh.copyfile(fileName,outFile)
+    try: rg=Dataset( outFile, 'r+' );
+    except: error('There was a problem opening "'+outFile+'".')
+    rgVar = rg.variables[variable] # handle to the variable
+    dims = rgVar.dimensions # tuple of dimensions
+    rgVar[:] = fullData.height[:,:] # Write the data
+    if All.edits.ijz:
+      print 'Applying %i edits'%(len(All.edits.ijz))
+      if 'nEdits' in rg.dimensions:
+        numEdits = rg.dimensions['nEdits']
+      else: numEdits = rg.createDimension('nEdits', 0)#len(All.edits.ijz))
+      if 'iEdit' in rg.variables: iEd = rg.variables['iEdit']
+      else:
+        iEd = rg.createVariable('iEdit','i4',('nEdits',))
+        iEd.long_name = 'i-index of edited data'
+      if 'jEdit' in rg.variables: jEd = rg.variables['jEdit']
+      else:
+        jEd = rg.createVariable('jEdit','i4',('nEdits',))
+        jEd.long_name = 'j-index of edited data'
+      if 'zEdit' in rg.variables: zEd = rg.variables['zEdit']
+      else:
+        zEd = rg.createVariable('zEdit','f4',('nEdits',))
+        zEd.long_name = 'Original value of edited data'
+        zEd.units = rgVar.units
+      for l,(i,j,z) in enumerate(All.edits.ijz):
+        iEd[l] = j; jEd[l] = i; zEd[l] = rgVar[i,j]; rgVar[i,j] = z
+    rg.close()
+
+
+def ice9it(i,j,depth):
+  # Iterative implementation of "ice 9"
+  wetMask = 0*depth
+  (ni,nj) = wetMask.shape
+  stack = set()
+  stack.add( (i,j) )
+  while stack:
+    (i,j) = stack.pop()
+    if wetMask[i,j] or depth[i,j] >= 0: continue
+    wetMask[i,j] = 1
+    if i>0: stack.add( (i-1,j) )
+    else: stack.add( (ni-1,j) )
+    if i<ni-1: stack.add( (i+1,j) )
+    else: stack.add( (0,j) )
+    if j>0: stack.add( (i,j-1) )
+    if j<nj-1: stack.add( (i,j+1) )
+  return wetMask
+
+
+def findPointInMesh(meshX, meshY, pointX, pointY):
+  def sign(x):
+    if x>0: return 1.0
+    elif x<0: return -1.0
+    else: return 0.
+  def crossProd(u0,v0,u1,v1):
+    return sign( u0*v1 - u1*v0 )
+  def isPointInConvexPolygon(pX, pY, p):
+    u0 = pX[0]-pX[-1]; v0 = pY[0]-pY[-1]
+    u1 = pX[-1] - p[0]; v1 = pY[-1] - p[1]
+    firstSign = crossProd(u0,v0,u1,v1)
+    for n in range(len(pX)-1):
+      u0 = pX[n+1]-pX[n]; v0 = pY[n+1]-pY[n]
+      u1 = pX[n] - p[0]; v1 = pY[n] - p[1]
+      if crossProd(u0,v0,u1,v1)*firstSign<0: return False
+    return True 
+  def recurIJ(mX, mY, p, ij00, ij22):
+    # Unpack indices
+    i0 = ij00[0]; i2 = ij22[0]; j0 = ij00[1]; j2 = ij22[1];
+    # Test bounding box first (bounding box is larger than polygon)
+    xmin=min( np.amin(mX[i0,j0:j2]), np.amin(mX[i2,j0:j2]), np.amin(mX[i0:i2,j0]), np.amin(mX[i0:i2,j2]) )
+    xmax=max( np.amax(mX[i0,j0:j2]), np.amax(mX[i2,j0:j2]), np.amax(mX[i0:i2,j0]), np.amax(mX[i0:i2,j2]) )
+    ymin=min( np.amin(mY[i0,j0:j2]), np.amin(mY[i2,j0:j2]), np.amin(mY[i0:i2,j0]), np.amin(mY[i0:i2,j2]) )
+    ymax=max( np.amax(mY[i0,j0:j2]), np.amax(mY[i2,j0:j2]), np.amax(mY[i0:i2,j0]), np.amax(mY[i0:i2,j2]) )
+    if p[0]<xmin or p[0]>xmax or p[1]<ymin or p[1]>ymax: return None, None
+    if i2>i0+1:
+      i1=int(0.5*(i0+i2))
+      if j2>j0+1: # Four quadrants to test
+        j1=int(0.5*(j0+j2))
+        iAns, jAns = recurIJ(mX, mY, p, (i0,j0), (i1,j1))
+        if iAns==None: iAns, jAns = recurIJ(mX, mY, p, (i1,j1), (i2,j2))
+        if iAns==None: iAns, jAns = recurIJ(mX, mY, p, (i0,j1), (i1,j2))
+        if iAns==None: iAns, jAns = recurIJ(mX, mY, p, (i1,j0), (i2,j1))
+      else: # Two halves, east/west, to test
+        j1=int(0.5*(j0+j2))
+        iAns, jAns = recurIJ(mX, mY, p, (i0,j0), (i1,j2))
+        if iAns==None: iAns, jAns = recurIJ(mX, mY, p, (i1,j0), (i2,j2))
+    else:
+      if j2>j0+1: # Two halves, north/south, to test
+        j1=int(0.5*(j0+j2))
+        iAns, jAns = recurIJ(mX, mY, p, (i0,j0), (i2,j1))
+        if iAns==None: iAns, jAns = recurIJ(mX, mY, p, (i0,j1), (i2,j2))
+      else: # Only one cell left (based on the bounding box)
+        if not isPointInConvexPolygon( \
+              [mX[i0,j0],mX[i0+1,j0],mX[i0+1,j0+1],mX[i0,j0+1]], \
+              [mY[i0,j0],mY[i0+1,j0],mY[i0+1,j0+1],mY[i0,j0+1]], \
+              p): return None, None
+        return i0,j0
+    return iAns, jAns
+  (ni,nj) = meshX.shape; ij00 = [0, 0]; ij22 = [ni-1, nj-1]
+  return recurIJ(meshX, meshY, (pointX, pointY), ij00, ij22)
+
+
+# Calculate a new window by scaling the current window, centering
+# on the cursor if possible.
+def newLims(cur_xlim, cur_ylim, cursor, xlim, ylim, scale_factor):
+  cur_xrange = (cur_xlim[1] - cur_xlim[0])*.5
+  cur_yrange = (cur_ylim[1] - cur_ylim[0])*.5
+  xdata = cursor[0]; ydata = cursor[1]
+  new_xrange = cur_xrange*scale_factor; new_yrange = cur_yrange*scale_factor
+  xdata = min( max( xdata, xlim[0]+new_xrange ), xlim[1]-new_xrange )
+  ydata = min( max( ydata, ylim[0]+new_yrange ), ylim[1]-new_yrange )
+  xL = max( xlim[0], xdata - new_xrange ); xR = min( xlim[1], xdata + new_xrange )
+  yL = max( ylim[0], ydata - new_yrange ); yR = min( ylim[1], ydata + new_yrange )
+  if xL==cur_xlim[0] and xR==cur_xlim[1] and \
+     yL==cur_ylim[0] and yR==cur_ylim[1]: return None, None
+  return (xL, xR), (yL, yR)
+
+
+# Class to handle adding buttons to GUI
+class Buttons:
+  scale = 0.014; space = .01
+  def __init__(self, bottom=.015, left=.015):
+    self.leftEdge = left
+    self.bottomEdge = bottom
+    self.height = .05
+    self.list = []
+  def add(self,label,fn): # fn is callback
+    width = self.scale*len(label)
+    np = [ self.leftEdge, self.bottomEdge, width, self.height ]
+    self.leftEdge = self.leftEdge + width + self.space
+    button = Button(plt.axes(np),label); button.on_clicked( fn )
+    self.list.append( button )
+
+
+# Class to contain edits
+class Edits:
+  def __init__(self):
+    self.newDepth = 0
+    self.ijz = []
+    self.label = None # Handle to text box
+  def setVal(self, newVal):
+    self.newDepth = newVal
+    if self.label: self.label.set_text('New depth = %i'%(self.newDepth))
+  def addToVal(self, increment):
+    self.newDepth += increment
+    if self.label: self.label.set_text('New depth = %i'%(self.newDepth))
+    plt.draw()
+  def get(self): return self.newDepth
+  def delete(self,i,j):
+    for I,J,D in self.ijz:
+      if (i,j)==(I,J): self.ijz.remove((I,J,D))
+  def add(self,i,j,nVal=None):
+    self.delete(i,j)
+    if not nVal==None: self.ijz.append( (i,j,nVal) )
+    else: self.ijz.append( (i,j,self.newDepth) )
+  def pop(self):
+    if self.ijz: self.ijz.pop()
+  def list(self):
+    for a in self.ijz: print a
+  def plot(self,topo):
+    x = []; y= []
+    for i,j,z in self.ijz:
+      tx,ty = topo.cellCoord(j,i)
+      if tx:
+        x.append(tx); y.append(ty)
+    if x:
+      h, = plt.plot(x, y, 'ro', hold=True)
+      return h
+    else: return None
+  def updatePlot(self,topo,h):
+    x = []; y= []
+    for i,j,z in self.ijz:
+      tx,ty = topo.cellCoord(j,i)
+      if tx:
+        x.append(tx); y.append(ty)
+    if x:
+      h.set_xdata(x); h.set_ydata(y)
+
+
+# Class to contain data
+class Topography:
+  def __init__(self, lon, lat, height):
+    self.longitude = lon
+    self.latitude = lat
+    self.height = np.copy( height )
+    self.xlim = ( np.min(lon), np.max(lon) )
+    self.ylim = ( np.min(lat), np.max(lat) )
+  def cloneWindow(self, (i0,j0), (iw,jw)):
+    i1 = i0 + iw; j1 = j0 + jw
+    return Topography( self.longitude[j0:j1+1,i0:i1+1], \
+                       self.latitude[j0:j1+1,i0:i1+1], \
+                       self.height[j0:j1,i0:i1] )
+  def applyEdits(self, origData, ijz):
+    for i,j,z in ijz:
+      x = (origData.longitude[i,j] + origData.longitude[i+1,j+1])/2.
+      y = (origData.latitude[i,j] + origData.latitude[i+1,j+1])/2.
+      (I,J) = findPointInMesh(self.longitude, self.latitude, x, y)
+      if not I==None: self.height[I,J] = z
+  def cellCoord(self,j,i):
+    #ni, nj = self.longitude.shape
+    #if i<0 or j<0 or i>=ni-1 or j>=nj-1: return None, None
+    x = (self.longitude[i,j] + self.longitude[i+1,j+1])/2.
+    y = (self.latitude[i,j] + self.latitude[i+1,j+1])/2.
+    return x,y
+
+# CLass to record the editing window
+class View:
+  def __init__(self, ni, nj):
+    self.ni = ni
+    self.nj = nj
+    self.i0 = 0
+    self.j0 = 0
+    self.iw = min(512,ni)
+    self.jw = min(512,nj)
+  def move(self, di, dj):
+    self.i0 = min( max(0, self.i0+int(di*self.iw/2.)), self.ni-self.iw)
+    self.j0 = min( max(0, self.j0+int(dj*self.jw/2.)), self.nj-self.jw)
+  def geti(self): return (self.i0,self.i0+self.iw)
+  def getj(self): return (self.j0,self.j0+self.jw)
+# Invoke main()
+if __name__ == '__main__': main()
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/generateZinterfaces.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/generateZinterfaces.py
new file mode 100644
index 0000000000..0bb65d23b2
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/generateZinterfaces.py
@@ -0,0 +1,13 @@
+import netCDF4
+import numpy
+
+dz=netCDF4.Dataset('../../../examples/ocean_SIS/MOM6z_SIS_025/INPUT/vgrid_cm4.nc').variables['dz'][:]
+zi=numpy.zeros(76)
+zi[1:]=numpy.cumsum(-dz)
+print zi
+
+Zbot = -netCDF4.Dataset('ocean_topog.nc').variables['depth'][:]
+
+Zi =numpy.zeros((76,1080,1440))
+for k in range(0,76):
+  Zi[k,:,:] = numpy.maximum( Zbot, zi[k] )
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/ice9.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/ice9.py
new file mode 100755
index 0000000000..f801926132
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/ice9.py
@@ -0,0 +1,150 @@
+#!/usr/bin/env python
+
+def error(msg,code=9):
+  print 'Error: ' + msg
+  exit(code)
+
+
+# Imports
+try: import argparse
+except: error('This version of python is not new enough. python 2.7 or newer is required.')
+try: from netCDF4 import Dataset
+except: error('Unable to import netCDF4 module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_netcdf4')
+try: import numpy as np
+except: error('Unable to import numpy module. Check your PYTHONPATH.\n'
+          +'Perhaps try:\n   module load python_numpy')
+import shutil as sh
+
+
+def main():
+
+  # Command line arguments
+  parser = argparse.ArgumentParser(description=
+       'Applies an "ice 9" algorithm to remove detached water from the topography. Also sets land elevation to 0.',
+       epilog='Written by A.Adcroft, 2013.')
+  parser.add_argument('filename', type=str,
+                      help='netcdf file to read.')
+  parser.add_argument('variable', type=str,
+                      nargs='?', default='depth',
+                      help='Name of variable to plot.')
+  parser.add_argument('--output', type=str,
+                      nargs='?', default=' ',
+                      help='name of the output file. If not specified, "iced_" is prepended to the name of the input file.')
+  parser.add_argument('--shallow', type=float,
+                      help='The "shallow" value (+ve, default 1.) to use when calculating the modified_mask. Wet points shallower than this are indicated with mask value of 2.')
+  parser.add_argument('--analyze', action='store_true',
+                      help='Report on impact of round shallow values to zero')
+
+  optCmdLineArgs = parser.parse_args()
+
+  nFileName = optCmdLineArgs.output
+  if nFileName == ' ': nFileName = 'iced_'+optCmdLineArgs.filename
+  shallow = 1
+  if not optCmdLineArgs.shallow==None: shallow = optCmdLineArgs.shallow
+  applyIce9(optCmdLineArgs.filename, nFileName, optCmdLineArgs.variable,
+            0., -40., shallow, optCmdLineArgs.analyze)
+
+def applyIce9(fileName, nFileName, variable, x0, y0, shallow, analyze):
+
+  iRg = Dataset( fileName, 'r' );
+  iDepth = iRg.variables[variable] # handle to the variable
+  depth = iDepth[:] # Read the data
+  print 'Range of input depths: min=',np.amin(depth),'max=',np.amax(depth)
+
+  # Open new netcdf file
+  if fileName==nFileName: error('Output file must be different from the input file')
+  try: rg=Dataset( nFileName, 'w', format='NETCDF3_CLASSIC' );
+  except: error('There was a problem opening "'+nFileName+'".')
+
+  (ny, nx) = depth.shape
+  rg.createDimension('nx',nx)
+  rg.createDimension('ny',ny)
+  rgDepth = rg.createVariable('depth','f4',('ny','nx'))
+  rgDepth.units = iDepth.units
+  rgDepth.standard_name = 'depth below geoid'
+  rgDepth.description = 'Non-negative nominal thickness of the ocean at cell centers'
+  rg.createDimension('ntiles',1)
+
+  # A mask based solely on value of depth
+  notLand = np.where( depth<0, 1, 0)
+  lons=iRg.variables['longitude'][:]
+  lats=iRg.variables['latitude'][:]
+  
+  iind=np.where(lons>-40.)[0][0]
+  jind=np.where(lats>30.)[0][0]
+  print jind,iind
+  print lons[iind],lats[jind]
+  notLand = ice9it(iind,jind,depth)
+
+  rgWet = rg.createVariable('wet','f4',('ny','nx'))
+  rgWet.long_name = 'Wet/dry mask'
+  rgWet.description = 'Values: 1=Ocean, 0=Land'
+  rgWet[:] = notLand # + (1-notLand)*0.3*np.where( depth<0, 1, 0)
+
+  rgDepth[:] = -depth*notLand # Change sign here. Until this point depth has actually been elevation.
+
+  if 'std' in iRg.variables: # Need to copy over list of edits
+    rgH2 = rg.createVariable('h2','f4',('ny','nx'))
+    rgH2.units = iDepth.units+'^2'
+    rgH2.standard_name = 'Variance of sub-grid scale topography'
+    rgH2[:] = iRg.variables['std'][:]**2
+
+  if 'zEdit' in iRg.variables: # Need to copy over list of edits
+    rgMod = rg.createVariable('modified_mask','f4',('ny','nx'))
+    rgMod.long_name = 'Modified mask'
+    rgMod.description = 'Values: 1=Ocean, 0=Land, -1 indicates water points removed by "Ice 9" algorithm. 2 indicates wet points that are shallower than 1m deep.'
+    rgMod[:] = notLand - (1-notLand)*np.where( depth<0, 1, 0) + np.where( (notLand>0) & (depth>-shallow), 1, 0)
+    n = len(iRg.variables['zEdit'])
+    nEd = rg.createDimension('nEdits',n)
+    iEd = rg.createVariable('iEdit','i4',('nEdits',))
+    iEd.long_name = 'i-index of edited data'
+    jEd = rg.createVariable('jEdit','i4',('nEdits',))
+    jEd.long_name = 'j-index of edited data'
+    zEd = rg.createVariable('zEdit','f4',('nEdits',))
+    zEd.long_name = 'Original value of height data'
+    zEd.units = iDepth.units
+    iEd[:] = iRg.variables['iEdit'][:]
+    jEd[:] = iRg.variables['jEdit'][:]
+    zEd[:] = iRg.variables['zEdit'][:]
+
+  rg.close()
+  print 'File "%s" written.'%(nFileName)
+
+  # Analyze the shallow points
+  if analyze:
+    print 'Analyzing...'
+    numNotLand = np.count_nonzero(notLand)
+    print '# of wet points after Ice 9 = %i'%(numNotLand)
+    newDepth = depth*np.where(depth*notLand <= -shallow, 1, 0)
+    numNewWet = np.count_nonzero(newDepth)
+    print '# of wet points deeper than %f = %i'%(-shallow,numNewWet)
+    print '%i - %i = %i fewer points left'%(numNotLand,numNewWet,numNotLand-numNewWet)
+    newWet = ice9it(600,270,newDepth)
+    numNewDeep = np.count_nonzero(newWet)
+    print '# of wet deep points after Ice 9 = %i'%(numNewDeep)
+    print '%i - %i = %i fewer points left'%(numNewWet,numNewDeep,numNewWet-numNewDeep)
+  
+
+def ice9it(i,j,depth):
+  # Iterative implementation of "ice 9"
+  wetMask = 0*depth
+  (nj,ni) = wetMask.shape
+  stack = set()
+  stack.add( (j,i) )
+  while stack:
+    (j,i) = stack.pop()
+    if wetMask[j,i] or depth[j,i] >= 0: continue
+    wetMask[j,i] = 1
+    if i>0: stack.add( (j,i-1) )
+    else: stack.add( (j,ni-1) )
+    if i<ni-1: stack.add( (j,i+1) )
+    else: stack.add( (0,j) )
+    if j>0: stack.add( (j-1,i) )
+    if j<nj-1: stack.add( (j+1,i) )
+    else: stack.add( (j,ni-1-i) )
+  return wetMask
+
+# Invoke main()
+if __name__ == '__main__': main()
+
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpCHL.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpCHL.py
new file mode 100755
index 0000000000..3e72465d51
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpCHL.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+
+from midas.rectgrid import *
+import netCDF4 as nc
+import numpy as np
+
+sgrid=supergrid(file='ocean_hgrid.nc')
+grid=quadmesh(supergrid=sgrid,cyclic=True)
+ogrid=quadmesh('/archive/gold/datasets/global/siena_201204/INPUT/seawifs_1998-2006_GOLD_smoothed_2X.nc',var='CHL_A',cyclic=True)
+O=state('/archive/gold/datasets/global/siena_201204/INPUT/seawifs_1998-2006_GOLD_smoothed_2X.nc',grid=ogrid,fields=['CHL_A'])
+O.var_dict['CHL_A']['Z']=None
+OM=O.horiz_interp('CHL_A',target=grid,method='bilinear')
+OM.rename_field('CHL_A','chl_a')
+OM.var_dict['chl_a']['xax_data']=grid.x_T[0,:]
+OM.var_dict['chl_a']['yax_data']=grid.y_T[:,grid.im/4]
+OM.write_nc('seawifs_1998-2006_smoothed_2X.nc',fields=['chl_a'])
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpSaltRestore.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpSaltRestore.py
new file mode 100755
index 0000000000..280494d309
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpSaltRestore.py
@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+
+from midas.rectgrid import *
+import numpy as np
+import netCDF4 as nc
+
+sgrid=supergrid(file='ocean_hgrid.nc',cyclic_x=True,tripolar_n=True)
+grid=quadmesh(supergrid=sgrid)
+grid.latq=grid.y_T_bounds[:,grid.im/4]
+grid.lath=grid.y_T[:,grid.im/4]
+
+O=state('/archive/gold/datasets/obs/WOA05_pottemp_salt.nc',fields=['SALT'],z_indices=np.arange(0,1),default_calendar='noleap')
+O.grid.cyclic_x=True
+O.var_dict['SALT']['Z'] = None
+print grid.latq.max(),grid.latq.min()
+print O.grid.latq.max(),O.grid.latq.min()
+OM=O.horiz_interp('SALT',target=grid,method='bilinear')
+OM.grid.D=nc.Dataset('topog.nc').variables['depth'][:]
+OM.grid.wet=np.zeros(OM.grid.D.shape)
+OM.grid.wet[OM.grid.D>0.]=1.
+OM.SALT=np.ma.masked_where(OM.SALT<0.,OM.SALT)
+OM.fill_interior('SALT',smooth=True,num_pass=10000)
+OM.mask_where('SALT','grid.D<=0.')
+OM.rename_field('SALT','salt')
+OM.var_dict['salt']['xax_data']=grid.x_T[0,:]
+OM.var_dict['salt']['yax_data']=grid.y_T[:,grid.im/4]
+OM.write_nc('salt_restore.nc',fields=['salt'])
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpTides.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpTides.py
new file mode 100755
index 0000000000..8a17f6985f
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpTides.py
@@ -0,0 +1,125 @@
+#!/usr/bin/env python
+
+from midas.rectgrid import *
+import netCDF4 as nc
+import numpy as np
+
+f=nc.Dataset('DATA/grid_tpxo7_atlas.nc')
+
+lon_u=f.variables['lon_u'][:].T
+lat_u=f.variables['lat_u'][:].T
+lon_v=f.variables['lon_v'][:].T
+lat_v=f.variables['lat_v'][:].T
+
+grid_u=quadmesh(lon=lon_u,lat=lat_u,cyclic=True)
+grid_v=quadmesh(lon=lon_v,lat=lat_v,cyclic=True)
+
+Su=state(grid=grid_u)
+Sv=state(grid=grid_v)
+
+f=nc.Dataset('DATA/u_tpxo7_atlas.nc')
+
+ua=np.zeros((1,8,721,1440))
+va=np.zeros((1,8,721,1440))
+
+for i in np.arange(0,8):
+    tmp=f.variables['ua'][i,:]
+    ua[0,i,:]=tmp.T
+    tmp=f.variables['va'][i,:]
+    va[0,i,:]=tmp.T    
+
+ua[ua==0.0]=-1.e20
+va[va==0.0]=-1.e20
+ua[:,:,-1,:] = ua[:,:,-2,:] # Work around for missing data at N-pole
+va[:,:,-1,:] = va[:,:,-2,:] # Work around for missing data at N-pole
+ua=np.ma.masked_where(ua==-1.e20,ua)
+va=np.ma.masked_where(va==-1.e20,va)
+
+
+vdict_u={}
+vdict_v={}
+
+vdict_u['X']=f.variables['ua'].dimensions[1]
+vdict_u['Y']=f.variables['ua'].dimensions[2]
+vdict_u['Z']=f.variables['ua'].dimensions[0][0:8]
+vdict_u['T']=None
+vdict_u['units']='cm s-1'
+vdict_u['path']='DATA/u_tpxo7_atlas.nc'
+vdict_u['Zdir']=1
+vdict_u['Ztype']='Fixed'
+vdict_u['Zb']=None
+vdict_u['z']=np.arange(0,8)
+vdict_u['z_interfaces']=None
+vdict_u['zunits']='none'
+vdict_u['zax_data']=np.arange(0,8)
+vdict_u['xax_data']=grid_u.lonh
+vdict_u['yax_data']=grid_u.lath
+vdict_u['xunits']='degrees_east'
+vdict_u['yunits']='degrees_north'
+vdict_u['zunits']='component'
+vdict_u['_FillValue']=-1.e20
+vdict_u['missing_value']=-1.e20
+vdict_u['masked']=True
+vdict_v['X']=f.variables['va'].dimensions[1]
+vdict_v['Y']=f.variables['va'].dimensions[2]
+vdict_v['Z']=f.variables['va'].dimensions[0][0:8]
+vdict_v['T']=None
+vdict_v['units']='cm s-1'
+vdict_v['path']='DATA/u_tpxo7_atlas.nc'
+vdict_v['Zdir']=1
+vdict_v['Ztype']='Fixed'
+vdict_v['Zb']=None
+vdict_v['z']=np.arange(0,8)
+vdict_v['z_interfaces']=None
+vdict_v['zunits']='none'
+vdict_v['zax_data']=np.arange(0,8)
+vdict_v['xax_data']=grid_v.lonh
+vdict_v['yax_data']=grid_v.lath
+vdict_v['xunits']='degrees_east'
+vdict_v['yunits']='degrees_north'
+vdict_v['zunits']='component'
+vdict_v['_FillValue']=-1.e20
+vdict_v['missing_value']=-1.e20
+vdict_v['masked']=True
+
+Su.add_field_from_array(ua,'ua',var_dict=vdict_u)
+Sv.add_field_from_array(va,'va',var_dict=vdict_v)
+
+grid=quadmesh(lon=lon_v,lat=lat_u,cyclic=True)
+grid.wet=np.ones((grid.jm,grid.im))
+
+S=Su.horiz_interp('ua',target=grid,method='bilinear')
+S=Sv.horiz_interp('va',target=grid,method='bilinear',PrevState=S)
+
+u2mod = (S.ua**2.0 + S.va**2.0)
+umod=np.sum(u2mod,axis=1)**0.5
+umod=umod[:,np.newaxis,:]
+umod = 1.e-2*umod
+
+vdict=S.var_dict['ua'].copy()
+vdict['units']='m s-1'
+vdict['Z'] = None
+
+S.add_field_from_array(umod,'umod',var_dict=vdict)
+
+sgrid=supergrid(file='ocean_hgrid.nc',cyclic_x=True,tripolar_n=True)
+output_grid = quadmesh(supergrid=sgrid,cyclic=True)
+output_grid.D=nc.Dataset('topog.nc').variables['depth'][:]
+output_grid.wet = np.zeros(output_grid.D.shape)
+output_grid.wet[output_grid.D>0.]=1.
+
+S.fill_interior('umod')
+
+print S.umod.shape
+
+T=S.horiz_interp('umod',target=output_grid,method='bilinear')
+
+T.fill_interior('umod')
+
+T.var_dict['umod']['yax_data'] = output_grid.y_T[:,output_grid.im/4]
+
+#T.umod = np.ma.filled(T.umod,0.0)
+
+T.rename_field('umod','tideamp')
+
+T.write_nc('tidal_amplitude.nc',['tideamp'])
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpWOA05.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpWOA05.py
new file mode 100755
index 0000000000..dde875af3a
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/interpWOA05.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+from midas.rectgrid import *
+import netCDF4 as nc
+import numpy as np
+
+sgrid=supergrid(file='ocean_hgrid.nc',cyclic_x=True,tripolar_n=True)
+grid=quadmesh(supergrid=sgrid)
+grid.lath=grid.y_T[:,grid.im/4]
+grid.latq=grid.y_T_bounds[:,grid.im/4]
+grid.D=nc.Dataset('topog.nc').variables['depth'][:]
+grid.wet=np.zeros(grid.D.shape)
+grid.wet[grid.D>0.]=1
+S=state(grid=grid)
+
+# Model vertical grid
+#dz=nc.Dataset('../../../examples/ocean_SIS/OM4_025/INPUT/vgrid_75_2m.nc').variables['dz'][:]
+#nk = dz.shape[0]
+#zi=np.zeros(nk+1)
+#zi[1:]=np.cumsum(-dz)
+# Analysis vertical  grid
+zi=-nc.Dataset('../../../ice_ocean_SIS/OM4_025/INPUT/vgrid_75_2m.nc').variables['zw'][:]
+nk=zi.shape[0]-1
+
+zb =np.zeros((nk+1,S.grid.jm,S.grid.im))
+for k in range(0,nk+1):
+  zb[k,:]=zi[k]
+  zb[k,:,:] = np.maximum( -S.grid.D, zb[k] )
+
+for n in np.arange(0,12):
+   O=state('/archive/gold/datasets/obs/WOA05_pottemp_salt.nc',fields=['SALT','PTEMP'],time_indices=np.arange(n,n+1),default_calendar='noleap',z_orientation=-1)
+   O.grid.cyclic_x=True
+   O.rename_field('PTEMP','ptemp')
+   O.rename_field('SALT','salt')
+   OM=O.horiz_interp('salt',target=S.grid,method='bilinear')
+   OM=O.horiz_interp('ptemp',target=S.grid,method='bilinear',PrevState=OM)
+   OM.adjust_thickness('ptemp')
+   OM.adjust_thickness('salt')
+   OM.fill_interior('salt',smooth=True,num_pass=10000)
+   OM.fill_interior('ptemp',smooth=True,num_pass=10000)
+
+   OM.remap_ALE(fields=['ptemp','salt'],z_bounds=zb,zbax_data=-zi,method='ppm_h4',bndy_extrapolation=False)
+   OM.rename_field('ptemp_remap','ptemp')
+   OM.rename_field('salt_remap','salt')
+   OM.mask_where('ptemp','grid.wet==0.')
+   OM.mask_where('salt','grid.wet==0.')
+   OM.ptemp=np.ma.masked_where(OM.var_dict['ptemp']['dz'][np.newaxis,:]<1.e-2, OM.ptemp)
+   OM.salt=np.ma.masked_where(OM.var_dict['ptemp']['dz'][np.newaxis,:]<1.e-2, OM.salt)
+
+   if n==0:
+      OM.write_nc('WOA05_ptemp_salt_monthly.nc',['ptemp','salt'],append=False,write_interface_positions=True)
+   else:
+      OM.write_nc('WOA05_ptemp_salt_monthly.nc',['ptemp','salt'],append=True,write_interface_positions=True)
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/make_basin_mask.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/make_basin_mask.py
new file mode 100755
index 0000000000..b73739d352
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/make_basin_mask.py
@@ -0,0 +1,198 @@
+#!/usr/bin/env python
+
+from midas.rectgrid import *
+from midas.rectgrid_gen import *
+import netCDF4
+import numpy
+
+def ice9it(i, j, depth, minD=0.):
+  """
+  Recursive implementation of "ice 9".
+  Returns 1 where depth>minD and is connected to depth[j,i], 0 otherwise.
+  """
+  wetMask = 0*depth      
+
+  (nj,ni) = wetMask.shape
+  stack = set()
+  stack.add( (j,i) )
+  while stack:
+    (j,i) = stack.pop()
+    if wetMask[j,i] or depth[j,i] <= minD: continue
+    wetMask[j,i] = 1
+  
+    if i>0: stack.add( (j,i-1) )
+    else: stack.add( (j,ni-1) ) # Periodic beyond i=0
+
+    if i<ni-1: stack.add( (j,i+1) )
+    else: stack.add( (j,0) ) # Periodic beyond i=ni-1
+    
+    if j>0: stack.add((j-1,i))
+        
+    if j<nj-1: stack.add( (j+1,i) )
+    else: stack.add( (j,ni-1-i) ) # Tri-polar fold beyond j=nj-1
+
+  return wetMask
+
+def ice9(x, y, depth, xy0):
+  ji = nearestJI(x, y, xy0)
+  return ice9it(ji[1], ji[0], depth)
+
+def nearestJI(x, y, (x0, y0)):
+  """
+  Find (j,i) of cell with center nearest to (x0,y0).
+  """
+  return numpy.unravel_index( ((x-x0)**2 + (y-y0)**2).argmin() , x.shape)
+
+def southOf(x, y, xy0, xy1):
+  """
+  Returns 1 for point south/east of the line that passes through xy0-xy1, 0 otherwise.
+  """
+  x0 = xy0[0]; y0 = xy0[1]; x1 = xy1[0]; y1 = xy1[1]
+  dx = x1 - x0; dy = y1 - y0
+  Y = (x-x0)*dy - (y-y0)*dx
+  Y[Y>=0] = 1; Y[Y<=0] = 0
+  return Y
+
+# Rewrite
+print 'Reading grid ...',
+x = netCDF4.Dataset('ocean_hgrid.nc').variables['x'][1::2,1::2] # Cell centers
+y = netCDF4.Dataset('ocean_hgrid.nc').variables['y'][1::2,1::2] # Cell centers
+print 'reading topography ...',
+depth = netCDF4.Dataset('topog.nc').variables['depth'][:]
+print 'done.'
+
+print 'Generating global wet mask ...',
+wet = ice9(x, y, depth, (0,-35)) # All ocean points seeded from South Atlantic
+print 'done.'
+
+code = 0*wet
+
+print 'Finding Cape of Good Hope ...',
+tmp = 1 - wet; tmp[x<-30] = 0
+tmp = ice9(x, y, tmp, (20,-30.))
+yCGH = (tmp*y).min()
+print 'done.', yCGH
+
+print 'Finding Melbourne ...',
+tmp = 1 - wet; tmp[x>-180] = 0
+tmp = ice9(x, y, tmp, (-220,-25.))
+yMel = (tmp*y).min()
+print 'done.', yMel
+
+print 'Processing Persian Gulf ...'
+tmp = wet*( 1-southOf(x, y, (55.,23.), (56.5,27.)) )
+tmp = ice9(x, y, tmp, (53.,25.))
+code[tmp>0] = 11
+wet = wet - tmp # Removed named points
+
+print 'Processing Red Sea ...'
+tmp = wet*( 1-southOf(x, y, (40.,11.), (45.,13.)) )
+tmp = ice9(x, y, tmp, (40.,18.))
+code[tmp>0] = 10
+wet = wet - tmp # Removed named points
+
+print 'Processing Black Sea ...'
+tmp = wet*( 1-southOf(x, y, (26.,42.), (32.,40.)) )
+tmp = ice9(x, y, tmp, (32.,43.))
+code[tmp>0] = 7
+wet = wet - tmp # Removed named points
+
+print 'Processing Mediterranean ...'
+tmp = wet*( southOf(x, y, (-5.7,35.5), (-5.7,36.5)) )
+tmp = ice9(x, y, tmp, (4.,38.))
+code[tmp>0] = 6
+wet = wet - tmp # Removed named points
+
+print 'Processing Baltic ...'
+tmp = wet*( southOf(x, y, (8.6,56.), (8.6,60.)) )
+tmp = ice9(x, y, tmp, (10.,58.))
+code[tmp>0] = 9
+wet = wet - tmp # Removed named points
+
+print 'Processing Hudson Bay ...'
+tmp = wet*( 
+           ( 1-(1-southOf(x, y, (-95.,66.), (-83.5,67.5)))
+              *(1-southOf(x, y, (-83.5,67.5), (-84.,71.))) 
+           )*( 1-southOf(x, y, (-70.,58.), (-70.,65.)) ) )
+tmp = ice9(x, y, tmp, (-85.,60.))
+code[tmp>0] = 8
+wet = wet - tmp # Removed named points
+
+print 'Processing Arctic ...'
+tmp = wet*( 
+          (1-southOf(x, y, (-171.,66.), (-166.,65.5))) * (1-southOf(x, y, (-64.,66.4), (-50.,68.5))) # Lab Sea
+     +    southOf(x, y, (-50.,0.), (-50.,90.)) * (1- southOf(x, y, (0.,65.5), (360.,65.5))  ) # Denmark Strait
+     +    southOf(x, y, (-18.,0.), (-18.,65.)) * (1- southOf(x, y, (0.,64.9), (360.,64.9))  ) # Iceland-Sweden
+     +    southOf(x, y, (20.,0.), (20.,90.)) # Barents Sea
+     +    (1-southOf(x, y, (-280.,55.), (-200.,65.)))
+          )
+tmp = ice9(x, y, tmp, (0.,85.))
+code[tmp>0] = 4
+wet = wet - tmp # Removed named points
+
+print 'Processing Pacific ...'
+tmp = wet*( (1-southOf(x, y, (0.,yMel), (360.,yMel)))
+           -southOf(x, y, (-257,1), (-257,0))*southOf(x, y, (0,3), (1,3))
+           -southOf(x, y, (-254.25,1), (-254.25,0))*southOf(x, y, (0,-5), (1,-5)) 
+           -southOf(x, y, (-243.7,1), (-243.7,0))*southOf(x, y, (0,-8.4), (1,-8.4)) 
+           -southOf(x, y, (-234.5,1), (-234.5,0))*southOf(x, y, (0,-8.9), (1,-8.9)) 
+          )
+tmp = ice9(x, y, tmp, (-150.,0.))
+code[tmp>0] = 3
+wet = wet - tmp # Removed named points
+
+print 'Processing Atlantic ...'
+tmp = wet*(1-southOf(x, y, (0.,yCGH), (360.,yCGH)))
+tmp = ice9(x, y, tmp, (-20.,0.))
+code[tmp>0] = 2
+wet = wet - tmp # Removed named points
+
+print 'Processing Indian ...'
+tmp = wet*(1-southOf(x, y, (0.,yCGH), (360.,yCGH)))
+tmp = ice9(x, y, tmp, (55.,0.))
+code[tmp>0] = 5
+wet = wet - tmp # Removed named points
+
+print 'Processing Southern Ocean ...'
+tmp = ice9(x, y, wet, (0.,-55.))
+code[tmp>0] = 1
+wet = wet - tmp # Removed named points
+
+code[wet>0] = -9
+(j,i) = numpy.unravel_index( wet.argmax(), x.shape)
+if j:
+  print 'There are leftover points unassigned to a basin code'
+  while j:
+    print x[j,i],y[j,i],[j,i]
+    wet[j,i]=0
+    (j,i) = numpy.unravel_index( wet.argmax(), x.shape)
+else: print 'All points assigned a basin code'
+
+sgrid=supergrid(file='ocean_hgrid.nc',cyclic_x=True,tripolar_n=True)
+grid=quadmesh(supergrid=sgrid)
+grid.D=netCDF4.Dataset('topog.nc').variables['depth'][:]
+grid.wet=numpy.zeros(grid.D.shape)
+grid.wet[grid.D>0.]=1.0
+grid.lath=grid.y_T[:,grid.im/4]  # should not be needed
+grid.latq=grid.y_T_bounds[:,grid.im/4+1] # ditto
+
+mask=code
+
+S=state(grid=grid)
+var_dict={}
+var_dict['X']='Longitude'
+var_dict['xax_data']= grid.lonh
+var_dict['xunits']= 'degrees_E'
+var_dict['Y']='Latitude'
+var_dict['yax_data']= grid.lath
+var_dict['yunits']= 'degrees_N'
+var_dict['Z']=None
+var_dict['T']=None
+var_dict['_FillValue']=None
+var_dict['missing_value']=None
+var_dict['flag_values']='1,2,3,4,5,6,7,8,9,10,11'
+var_dict['flag_meanings']='1:Southern Ocean, 2:Atlantic Ocean, 3:Pacific Ocean, 4:Arctic Ocean, 5:Indian Ocean, 6:Mediterranean Sea, 7:Black Sea, 8:Hudson Bay, 9:Baltic Sea, 10:Red Sea, 11:Persian Gulf'
+
+S.add_field_from_array(mask,'basin',var_dict=var_dict)
+
+S.write_nc('basin_codes.nc',['basin'])
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/md5sums.txt b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/md5sums.txt
new file mode 100644
index 0000000000..bf1f9774a1
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/md5sums.txt
@@ -0,0 +1,24 @@
+Grids
+9b65729ae7c31763e610ceb55ee13d10  sgridA.nc
+770370f885f74e5a3d248a5ba626fdbc  sgridB.nc
+a3b5fc9ecc173c0377866e888580bd35  ocean_hgrid.nc
+
+Topography
+e58b478085f3e7b3df6d7a1cafd2a1f0  edit_topog.nc
+66c17af456a3305e6131889d6d2a6923  interpolated_topog.nc
+ca830bb2f054c89aaa3d3dad52a4de6c  topog.nc
+d3df7dd668eda3e16f131147bda6e54b  basin_codes.nc
+
+Data
+fee749b92d6bd7bcdb569e15c5d60d31  salt_restore.nc
+ae3d8ec17033428adc881a2f2af79e1e  seawifs_1998-2006_smoothed_2X.nc
+d0dcef39a65d71b1a7d6e1ebec9b06d1  tidal_amplitude.nc
+
+Obs
+708c256e3a663ca93af44742693c1d94  WOA05_ptemp_salt_monthly.nc
+e5978a20dcbe8d805a2aebc391fcbd19  WOA05_ptemp_salt_annual.nc
+
+Mosaic
+355e70e71d77feb34342059894fda960  C180_mosaic/mosaic.tar
+Quick Mosaic
+617876e77f1e8465fcebaf4b703a9bc9  C180_mosaic/quick_mosaic.tar
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/merge_grids.py b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/merge_grids.py
new file mode 100755
index 0000000000..2c81fdadd8
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/merge_grids.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+import netCDF4 as nc
+import numpy as np
+
+f=nc.Dataset('sgridA.nc')
+g=nc.Dataset('sgridB.nc')
+y1=f.variables['y'][:]
+y2=g.variables['y'][:]
+y=np.concatenate((y1,y2[1:,:]),axis=0)
+
+dy1=f.variables['dy'][:]
+dy2=g.variables['dy'][:]
+dy=np.concatenate((dy1,dy2),axis=0)
+
+x1=f.variables['x'][:]
+x2=g.variables['x'][:]
+x=np.concatenate((x1,x2[1:,:]),axis=0)
+
+dx1=f.variables['dx'][:]
+dx2=g.variables['dx'][:]
+dx=np.concatenate((dx1,dx2[1:,:]),axis=0)
+
+
+area1=f.variables['area'][:]
+area2=g.variables['area'][:]
+area=np.concatenate((area1,area2),axis=0)
+
+angle_dx1=f.variables['angle_dx'][:,:]
+angle_dx2=g.variables['angle_dx'][:-1,:]
+angle_dx=np.concatenate((angle_dx1,angle_dx2),axis=0)
+
+fout=nc.Dataset('ocean_hgrid.nc','w',format='NETCDF3_CLASSIC')
+
+ny=area.shape[0]; nx = area.shape[1]
+nyp=ny+1; nxp=nx+1
+
+print 'ny,nx= ',ny,nx
+
+nyp=fout.createDimension('nyp',nyp)
+nxp=fout.createDimension('nxp',nxp)
+ny=fout.createDimension('ny',ny)
+nx=fout.createDimension('nx',nx)
+string=fout.createDimension('string',255)    
+
+tile=fout.createVariable('tile','S1',('string'))
+yv=fout.createVariable('y','f8',('nyp','nxp'))
+xv=fout.createVariable('x','f8',('nyp','nxp'))    
+yv.units='degrees'
+xv.units='degrees'
+yv[:]=y
+xv[:]=x
+
+tile[0:4]='tile1'
+dyv=fout.createVariable('dy','f8',('ny','nxp'))
+dyv.units='meters'
+dyv[:]=dy
+dxv=fout.createVariable('dx','f8',('nyp','nx'))
+dxv.units='meters'
+dxv[:]=dx
+areav=fout.createVariable('area','f8',('ny','nx'))
+areav.units='m2'
+areav[:]=area
+anglev=fout.createVariable('angle_dx','f8',('nyp','nxp'))
+anglev.units='degrees'
+anglev[:]=angle_dx            
+
+fout.sync()
+fout.close()
diff --git a/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/topo_edits b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/topo_edits
new file mode 100644
index 0000000000..87e3a9c06d
--- /dev/null
+++ b/land_ice_ocean_LM3_SIS2/OM_360x320_C180/preprocessing/topo_edits
@@ -0,0 +1,47 @@
+(196, 219, 130.0)
+(206, 294, -0)
+(206, 295, -0)
+(206, 296, -0)
+(206, 298, -0)
+(206, 297, -0)
+(208, 297, 100.0)
+(208, 296, 100.0)
+(208, 295, 100.0)
+(213, 326, 130.0)
+(212, 326, 130.0)
+(181, 342, 130.0)
+(181, 343, 130.0)
+(196, 357, 130.0)
+(195, 357, 130.0)
+(196, 355, 130.0)
+(219, 336, 130.0)
+(228, 300, 130.0)
+(152, 43, 160.0)
+(154, 42, 160.0)
+(125, 55, 190.0)
+(124, 62, -0)
+(124, 63, -0)
+(148, 60, -0)
+(147, 60, -0)
+(146, 60, -0)
+(231, 81, 30.0)
+(230, 81, 30.0)
+(229, 81, 30.0)
+(215, 80, 30.0)
+(191, 218, 0)
+(292, 249, 200)
+(293, 251, 200)
+(294, 251, 200)
+(207, 293, 600)
+(207, 294, 600)
+(256, 273, 600)
+(256, 272, 600)
+(255, 272, 600)
+(180, 343, 150)
+(182, 342, 150)
+(183, 341, 150)
+(196, 356, 150)
+(236, 311, 30)
+(234, 311, 30)
+(235, 311, 30)
+(234, 312, 30)