From 97988bbd61942388ce31c497bd2bd6a6c075dccd Mon Sep 17 00:00:00 2001
From: Kaitlin Alexander <kaalexander@users.noreply.github.com>
Date: Wed, 18 Nov 2015 14:41:05 +1100
Subject: [PATCH] Made initialisation files more modular

---
 calc_z.pyc               | Bin 0 -> 984 bytes
 file_guide.txt           |  27 ++-
 romscice_ini.py          | 368 ++++++++++++++++++++-------------------
 romscice_ini_iceshelf.py |  56 +++---
 romscice_ini_woa.py      | 303 ++++++++++++++++----------------
 woa_netcdf.py            | 182 +++++++++----------
 6 files changed, 484 insertions(+), 452 deletions(-)
 create mode 100644 calc_z.pyc

diff --git a/calc_z.pyc b/calc_z.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..bf018bf428b12cff8784f747e1b65a4c0b672e03
GIT binary patch
literal 984
zcmb7B&2G~`5T13?Bu(?z#ID1ILy?fxLqiT!94ZP5=Nu?M5)z6W$HY;cIAk}4M!ToN
zbMR_B0um1Z->egIX0<z>-*3M06n@vMKl(4<MYKH*zAx~k4>7?NloGv}IP{I^gP=))
zu7e_Cvp+b2fTmcSMHR>syfWSZPu{|ha#OMLkfkn@XR4$SaM*zK?;Qs1nj_M<W$NlG
zxu(Xjx=6Z4?!BNR&dfcqV01lGCk2oLPT?r0<qjwwr<lerrYqEe2OoV%p6{}rSHbRd
zUDG6uTPGW7f<uE`-?T_$yaH)K!0Z|rU<AGnfwei_nMG^1Nr6n6ber@seJXH`>Co4R
zREhmCUwwQ_#${9uY?mIm8lhUQ`s^-|?-|c#c;xPT>kcRJLWsVD!A?P&%g-(Gp5wr2
zao_9U0RMkDg~0iBHRA<st)Evjb_Tp9&>{*~;bIzNE}|%tiB;rmk)Sn#W4U!~Vhl&L
ziC3Alnk@@>;dcLd;bs=ailNhNF<K>Au7~QMu98@VBWeA)KE}?7lE|VAEp(VJlhF1Z
zBiqI&((*~;Tbq}&8p!p^rVYlKvcr6dXVXO4l!5KBSkFh`ve&~ry0JMo7JKs~00Xi@
zK8DtlNw!kG`Ngx|O&C7C==WCHT#i(hrLq@=X*AlLFYj!_-0Gl6o-#<eBub(xe9>?o
ziBoYVPMnt06iwlY4(z@&40v2SCb%u{UO$An@BiN<Dk#64XOndrKj$?|hOf5&0#;GA
AM*si-

literal 0
HcmV?d00001

diff --git a/file_guide.txt b/file_guide.txt
index a868a08..072a0d9 100644
--- a/file_guide.txt
+++ b/file_guide.txt
@@ -36,12 +36,11 @@ romscice_ini.py: Builds a ROMS initialisation file using ECCO2 data for
                  temperature and salinity; starts with a motionless ocean and
 		 zero sea surface height. Depends on calc_z.py.
 		 To run: First edit user parameters (file paths and grid
-		         parameters) partway through the file (after the
-			 function interp_ecco2roms). Then make sure that you
-			 have scipy version 0.14 or higher (on raijin, this
-			 means switching to python/2.7.6; instructions at the
-			 top of the file). Then open python/ipython and type
-			 "run romscice_ini.py".
+		         parameters) near the bottom of the file. Then make
+			 sure that you have scipy version 0.14 or higher (on
+			 raijin, this means switching to python/2.7.6;
+			 instructions at the top of the file). Then open
+			 python/ipython and type "run romscice_ini.py".
 
 woa_netcdf.py: Converts World Ocean Atlas temperature and salinity data from
                text file (FESOM input format) to NetCDF.
@@ -55,12 +54,12 @@ romscice_ini_woa.py: Builds a ROMS initialisation file using World Ocean Atlas
 		     and uses WOA data assumed to be converted from FESOM input
 		     using woa_netcdf.py
 		     To run: First edit user parameters (file paths and grid
-		             parameters) partway through the file (after the
-			     function interp_woa2roms). Then make sure that you
-			     have scipy version 0.14 or higher (on raijin, this
-			     means switching to python/2.7.6; instructions at
-			     the top of the file). Then open python/ipython and
-			     type "run romscice_ini_woa.py".
+		             parameters) near the bottom of the file (after the
+			     Then make sure that you have scipy version 0.14 or
+			     higher (on raijin, this  means switching to
+			     python/2.7.6; instructions at the top of the file).
+			     Then open python/ipython and type "run
+			     romscice_ini_woa.py".
 
 convert_ecco.job: A batch job which converts 1 year of ECCO2 data into
                   northern lateral bounday conditions for ROMS. Depends on
@@ -154,8 +153,8 @@ romscice_ini_iceshelf.py: Given a ROMS initial condition file from ECCO2 data
 			  it's more realistic than filling the cavities with
 			  constant salinity, and certainly more stable).
 			  To run: Edit the user parameters (paths to 3 files)
-			          near the top of the script, then open python
-				  or ipython and type "run
+			          near the bottom of the script, then open
+				  python or ipython and type "run
 				  romscice_ini_iceshelf.py".
 
 
diff --git a/romscice_ini.py b/romscice_ini.py
index 2eb89ae..7679392 100644
--- a/romscice_ini.py
+++ b/romscice_ini.py
@@ -16,6 +16,176 @@
 from scipy.interpolate import RegularGridInterpolator
 from calc_z import *
 
+
+# Main routine
+def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b, hc, N, nbdry_ecco):
+
+    # Read ECCO2 data and grid
+    print 'Reading ECCO2 data'
+    theta_fid = Dataset(theta_file, 'r')
+    lon_ecco_raw = theta_fid.variables['LONGITUDE_T'][:]
+    lat_ecco = theta_fid.variables['LATITUDE_T'][0:nbdry_ecco]
+    depth_ecco_raw = theta_fid.variables['DEPTH_T'][:]
+    theta_raw = transpose(theta_fid.variables['THETA'][0,:,0:nbdry_ecco,:])
+    theta_fid.close()
+    salt_fid = Dataset(salt_file, 'r')
+    salt_raw = transpose(salt_fid.variables['SALT'][0,:,0:nbdry_ecco,:])
+    salt_fid.close()
+
+    # The ECCO2 longitude axis doesn't wrap around; there is a gap between
+    # almost-180W and almost-180E, and the ROMS grid has points in this gap.
+    # So copy the last longitude value (mod 360) to the beginning, and the
+    # first longitude value (mod 360) to the end.
+    lon_ecco = zeros(size(lon_ecco_raw)+2)
+    lon_ecco[0] = lon_ecco_raw[-1]-360
+    lon_ecco[1:-1] = lon_ecco_raw
+    lon_ecco[-1] = lon_ecco_raw[0]+360
+
+    # The shallowest ECCO2 depth value is 5 m, but ROMS needs 0 m. So add the
+    # index depth = 0 m to the beginning. Later we will just copy the 5 m values
+    # for theta and salt into this index. Similarly, add the index depth = 6000 m
+    # to the end.
+    depth_ecco = zeros(size(depth_ecco_raw)+2)
+    depth_ecco[0] = 0.0
+    depth_ecco[1:-1] = depth_ecco_raw
+    depth_ecco[-1] = 6000.0
+
+    # Copy the theta and salt values to the new longitude and depth indices,
+    # making sure to preserve the mask.
+    theta = ma.array(zeros((size(lon_ecco), size(lat_ecco), size(depth_ecco))))
+    theta[1:-1,:,1:-1] = ma.copy(theta_raw)
+    theta[0,:,1:-1] = ma.copy(theta_raw[-1,:,:])
+    theta[-1,:,1:-1] = ma.copy(theta_raw[0,:,:])
+    theta[:,:,0] = ma.copy(theta[:,:,1])
+    theta[:,:,-1] = ma.copy(theta[:,:,-2])
+    salt = ma.array(zeros((size(lon_ecco), size(lat_ecco), size(depth_ecco))))
+    salt[1:-1,:,1:-1] = ma.copy(salt_raw)
+    salt[0,:,1:-1] = ma.copy(salt_raw[-1,:,:])
+    salt[-1,:,1:-1] = ma.copy(salt_raw[0,:,:])
+    salt[:,:,0] = ma.copy(salt[:,:,1])
+    salt[:,:,-1] = ma.copy(salt[:,:,-2])
+
+    # Read ROMS grid
+    print 'Reading ROMS grid'
+    grid_fid = Dataset(grid_file, 'r')
+    lon_roms = grid_fid.variables['lon_rho'][:,:]
+    lat_roms = grid_fid.variables['lat_rho'][:,:]
+    h = grid_fid.variables['h'][:,:]
+    zice = grid_fid.variables['zice'][:,:]
+    mask_rho = grid_fid.variables['mask_rho'][:,:]
+    mask_zice = grid_fid.variables['mask_zice'][:,:]
+    grid_fid.close()
+    num_lon = size(lon_roms, 1)
+    num_lat = size(lon_roms, 0)
+    # Mask h and zice with zeros
+    h = h*mask_rho
+    zice = zice*mask_zice
+
+    # Get a 3D array of ROMS z-coordinates, as well as 1D arrays of s-coordinates
+    # and stretching curves
+    z_roms_3d, sc_r, Cs_r = calc_z(h, zice, lon_roms, lat_roms, theta_s, theta_b, hc, N)
+    # Copy the latitude and longitude values into 3D arrays of the same shape
+    lon_roms_3d = tile(lon_roms, (N,1,1))
+    lat_roms_3d = tile(lat_roms, (N,1,1))
+
+    # Regridding happens here...
+    print 'Interpolating temperature'
+    temp = interp_ecco2roms(theta, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, -0.5)
+    print 'Interpolating salinity'
+    salt = interp_ecco2roms(salt, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, 34.5)
+
+    # Set initial velocities and sea surface height to zero
+    u = zeros((N, num_lat, num_lon-1))
+    v = zeros((N, num_lat-1, num_lon))
+    ubar = zeros((num_lat, num_lon-1))
+    vbar = zeros((num_lat-1, num_lon))
+    zeta = zeros((num_lat, num_lon))
+
+    print 'Writing to NetCDF file'
+    out_fid = Dataset(output_file, 'w')
+    # Define dimensions
+    out_fid.createDimension('xi_u', num_lon-1)
+    out_fid.createDimension('xi_v', num_lon)
+    out_fid.createDimension('xi_rho', num_lon)
+    out_fid.createDimension('eta_u', num_lat)
+    out_fid.createDimension('eta_v', num_lat-1)
+    out_fid.createDimension('eta_rho', num_lat)
+    out_fid.createDimension('s_rho', N)
+    out_fid.createDimension('ocean_time', None)
+    out_fid.createDimension('one', 1);
+    # Define variables and assign values
+    out_fid.createVariable('tstart', 'f8', ('one'))
+    out_fid.variables['tstart'].long_name = 'start processing day'
+    out_fid.variables['tstart'].units = 'day'
+    out_fid.variables['tstart'][:] = 0.0
+    out_fid.createVariable('tend', 'f8', ('one'))
+    out_fid.variables['tend'].long_name = 'end processing day'
+    out_fid.variables['tend'].units = 'day'
+    out_fid.variables['tend'][:] = 0.0
+    out_fid.createVariable('theta_s', 'f8', ('one'))
+    out_fid.variables['theta_s'].long_name = 'S-coordinate surface control parameter'
+    out_fid.variables['theta_s'][:] = theta_s
+    out_fid.createVariable('theta_b', 'f8', ('one'))
+    out_fid.variables['theta_b'].long_name = 'S-coordinate bottom control parameter'
+    out_fid.variables['theta_b'].units = 'nondimensional'
+    out_fid.variables['theta_b'][:] = theta_b
+    out_fid.createVariable('Tcline', 'f8', ('one'))
+    out_fid.variables['Tcline'].long_name = 'S-coordinate surface/bottom layer width'
+    out_fid.variables['Tcline'].units = 'meter'
+    out_fid.variables['Tcline'][:] = Tcline
+    out_fid.createVariable('hc', 'f8', ('one'))
+    out_fid.variables['hc'].long_name = 'S-coordinate parameter, critical depth'
+    out_fid.variables['hc'].units = 'meter'
+    out_fid.variables['hc'][:] = hc
+    out_fid.createVariable('Cs_r', 'f8', ('s_rho'))
+    out_fid.variables['Cs_r'].long_name = 'S-coordinate stretching curves at RHO-points'
+    out_fid.variables['Cs_r'].units = 'nondimensional'
+    out_fid.variables['Cs_r'].valid_min = -1.0
+    out_fid.variables['Cs_r'].valid_max = 0.0
+    out_fid.variables['Cs_r'][:] = Cs_r
+    out_fid.createVariable('ocean_time', 'f8', ('ocean_time'))
+    out_fid.variables['ocean_time'].long_name = 'time since initialization'
+    out_fid.variables['ocean_time'].units = 'seconds'
+    out_fid.variables['ocean_time'][0] = 0.0
+    out_fid.createVariable('u', 'f8', ('ocean_time', 's_rho', 'eta_u', 'xi_u'))
+    out_fid.variables['u'].long_name = 'u-momentum component'
+    out_fid.variables['u'].units = 'meter second-1'
+    out_fid.variables['u'][0,:,:,:] = u
+    out_fid.createVariable('v', 'f8', ('ocean_time', 's_rho', 'eta_v', 'xi_v'))
+    out_fid.variables['v'].long_name = 'v-momentum component'
+    out_fid.variables['v'].units = 'meter second-1'
+    out_fid.variables['v'][0,:,:,:] = v
+    out_fid.createVariable('ubar', 'f8', ('ocean_time', 'eta_u', 'xi_u'))
+    out_fid.variables['ubar'].long_name = 'vertically integrated u-momentum component'
+    out_fid.variables['ubar'].units = 'meter second-1'
+    out_fid.variables['ubar'][0,:,:] = ubar
+    out_fid.createVariable('vbar', 'f8', ('ocean_time', 'eta_v', 'xi_v'))
+    out_fid.variables['vbar'].long_name = 'vertically integrated v-momentum component'
+    out_fid.variables['vbar'].units = 'meter second-1'
+    out_fid.variables['vbar'][0,:,:] = vbar
+    out_fid.createVariable('zeta', 'f8', ('ocean_time', 'eta_rho', 'xi_rho'))
+    out_fid.variables['zeta'].long_name = 'free-surface'
+    out_fid.variables['zeta'].units = 'meter'
+    out_fid.variables['zeta'][0,:,:] = zeta
+    out_fid.createVariable('temp', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
+    out_fid.variables['temp'].long_name = 'potential temperature'
+    out_fid.variables['temp'].units = 'Celsius'
+    out_fid.variables['temp'][0,:,:,:] = temp
+    out_fid.createVariable('salt', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
+    out_fid.variables['salt'].long_name = 'salinity'
+    out_fid.variables['salt'].units = 'PSU'
+    out_fid.variables['salt'][0,:,:,:] = salt
+    out_fid.createVariable('sc_r', 'f8', ('s_rho'))
+    out_fid.variables['sc_r'].long_name = 'S-coordinate at rho-points'
+    out_fid.variables['sc_r'].units = 'nondimensional'
+    out_fid.variables['sc_r'].valid_min = -1.0
+    out_fid.variables['sc_r'].valid_max = 0.0
+    out_fid.variables['sc_r'][:] = sc_r
+    out_fid.close()
+
+    print 'Finished'
+
+
 # Given an array on the ECCO2 grid, fill in the land mask with constant values,
 # and then interpolate to the ROMS grid.
 # Input:
@@ -67,189 +237,31 @@ def interp_ecco2roms (A, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3
     return B
 
 
-# File paths to edit here:
-# Path to ROMS grid file
-grid_file = '/short/m68/kaa561/roms_circumpolar/data/caisom001_OneQuartergrd.nc'
-# Path to ECCO2 files for temperature and salinity in January 1995
-theta_file = '../data/ECCO2/THETA.1440x720x50.199501.nc'
-salt_file = '../data/ECCO2/SALT.1440x720x50.199501.nc'
-# Path to desired output file
-output_file = '../data/caisom001_ini_1995.nc'
-
-# Grid parameters; check grid_file and *.in file to make sure these are correct
-Tcline = 40
-theta_s = 0.9
-theta_b = 4
-hc = 40
-N = 31
-
-# Northernmost index of ECCO2 grid to read (1-based)
-nbdry_ecco = 161
-
-# Read ECCO2 data and grid
-print 'Reading ECCO2 data'
-theta_fid = Dataset(theta_file, 'r')
-lon_ecco_raw = theta_fid.variables['LONGITUDE_T'][:]
-lat_ecco = theta_fid.variables['LATITUDE_T'][0:nbdry_ecco]
-depth_ecco_raw = theta_fid.variables['DEPTH_T'][:]
-theta_raw = transpose(theta_fid.variables['THETA'][0,:,0:nbdry_ecco,:])
-theta_fid.close()
-salt_fid = Dataset(salt_file, 'r')
-salt_raw = transpose(salt_fid.variables['SALT'][0,:,0:nbdry_ecco,:])
-salt_fid.close()
+# User parameters
+if __name__ == "__main__":
 
-# The ECCO2 longitude axis doesn't wrap around; there is a gap between
-# almost-180W and almost-180E, and the ROMS grid has points in this gap.
-# So copy the last longitude value (mod 360) to the beginning, and the
-# first longitude value (mod 360) to the end.
-lon_ecco = zeros(size(lon_ecco_raw)+2)
-lon_ecco[0] = lon_ecco_raw[-1]-360
-lon_ecco[1:-1] = lon_ecco_raw
-lon_ecco[-1] = lon_ecco_raw[0]+360
+    # File paths to edit here:
+    # Path to ROMS grid file
+    grid_file = '/short/m68/kaa561/roms_circumpolar/data/caisom001_OneQuartergrd.nc'
+    # Path to ECCO2 files for temperature and salinity in January 1995
+    theta_file = '../data/ECCO2/THETA.1440x720x50.199501.nc'
+    salt_file = '../data/ECCO2/SALT.1440x720x50.199501.nc'
+    # Path to desired output file
+    output_file = '../data/caisom001_ini_1995.nc'
 
-# The shallowest ECCO2 depth value is 5 m, but ROMS needs 0 m. So add the
-# index depth = 0 m to the beginning. Later we will just copy the 5 m values
-# for theta and salt into this index. Similarly, add the index depth = 6000 m
-# to the end.
-depth_ecco = zeros(size(depth_ecco_raw)+2)
-depth_ecco[0] = 0.0
-depth_ecco[1:-1] = depth_ecco_raw
-depth_ecco[-1] = 6000.0
+    # Grid parameters; check grid_file and *.in file to make sure these are correct
+    Tcline = 40
+    theta_s = 0.9
+    theta_b = 4
+    hc = 40
+    N = 31
 
-# Copy the theta and salt values to the new longitude and depth indices,
-# making sure to preserve the mask.
-theta = ma.array(zeros((size(lon_ecco), size(lat_ecco), size(depth_ecco))))
-theta[1:-1,:,1:-1] = ma.copy(theta_raw)
-theta[0,:,1:-1] = ma.copy(theta_raw[-1,:,:])
-theta[-1,:,1:-1] = ma.copy(theta_raw[0,:,:])
-theta[:,:,0] = ma.copy(theta[:,:,1])
-theta[:,:,-1] = ma.copy(theta[:,:,-2])
-salt = ma.array(zeros((size(lon_ecco), size(lat_ecco), size(depth_ecco))))
-salt[1:-1,:,1:-1] = ma.copy(salt_raw)
-salt[0,:,1:-1] = ma.copy(salt_raw[-1,:,:])
-salt[-1,:,1:-1] = ma.copy(salt_raw[0,:,:])
-salt[:,:,0] = ma.copy(salt[:,:,1])
-salt[:,:,-1] = ma.copy(salt[:,:,-2])
+    # Northernmost index of ECCO2 grid to read (1-based)
+    nbdry_ecco = 161
 
-# Read ROMS grid
-print 'Reading ROMS grid'
-grid_fid = Dataset(grid_file, 'r')
-lon_roms = grid_fid.variables['lon_rho'][:,:]
-lat_roms = grid_fid.variables['lat_rho'][:,:]
-h = grid_fid.variables['h'][:,:]
-zice = grid_fid.variables['zice'][:,:]
-mask_rho = grid_fid.variables['mask_rho'][:,:]
-mask_zice = grid_fid.variables['mask_zice'][:,:]
-grid_fid.close()
-num_lon = size(lon_roms, 1)
-num_lat = size(lon_roms, 0)
-# Mask h and zice with zeros
-h = h*mask_rho
-zice = zice*mask_zice
+    run(grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b, hc, N, nbdry_ecco)
 
-# Get a 3D array of ROMS z-coordinates, as well as 1D arrays of s-coordinates
-# and stretching curves
-z_roms_3d, sc_r, Cs_r = calc_z(h, zice, lon_roms, lat_roms, theta_s, theta_b, hc, N)
-# Copy the latitude and longitude values into 3D arrays of the same shape
-lon_roms_3d = tile(lon_roms, (N,1,1))
-lat_roms_3d = tile(lat_roms, (N,1,1))
-
-# Regridding happens here...
-print 'Interpolating temperature'
-temp = interp_ecco2roms(theta, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, -0.5)
-print 'Interpolating salinity'
-salt = interp_ecco2roms(salt, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, 34.5)
-
-# Set initial velocities and sea surface height to zero
-u = zeros((N, num_lat, num_lon-1))
-v = zeros((N, num_lat-1, num_lon))
-ubar = zeros((num_lat, num_lon-1))
-vbar = zeros((num_lat-1, num_lon))
-zeta = zeros((num_lat, num_lon))
-
-print 'Writing to NetCDF file'
-out_fid = Dataset(output_file, 'w')
-# Define dimensions
-out_fid.createDimension('xi_u', num_lon-1)
-out_fid.createDimension('xi_v', num_lon)
-out_fid.createDimension('xi_rho', num_lon)
-out_fid.createDimension('eta_u', num_lat)
-out_fid.createDimension('eta_v', num_lat-1)
-out_fid.createDimension('eta_rho', num_lat)
-out_fid.createDimension('s_rho', N)
-out_fid.createDimension('ocean_time', None)
-out_fid.createDimension('one', 1);
-# Define variables and assign values
-out_fid.createVariable('tstart', 'f8', ('one'))
-out_fid.variables['tstart'].long_name = 'start processing day'
-out_fid.variables['tstart'].units = 'day'
-out_fid.variables['tstart'][:] = 0.0
-out_fid.createVariable('tend', 'f8', ('one'))
-out_fid.variables['tend'].long_name = 'end processing day'
-out_fid.variables['tend'].units = 'day'
-out_fid.variables['tend'][:] = 0.0
-out_fid.createVariable('theta_s', 'f8', ('one'))
-out_fid.variables['theta_s'].long_name = 'S-coordinate surface control parameter'
-out_fid.variables['theta_s'][:] = theta_s
-out_fid.createVariable('theta_b', 'f8', ('one'))
-out_fid.variables['theta_b'].long_name = 'S-coordinate bottom control parameter'
-out_fid.variables['theta_b'].units = 'nondimensional'
-out_fid.variables['theta_b'][:] = theta_b
-out_fid.createVariable('Tcline', 'f8', ('one'))
-out_fid.variables['Tcline'].long_name = 'S-coordinate surface/bottom layer width'
-out_fid.variables['Tcline'].units = 'meter'
-out_fid.variables['Tcline'][:] = Tcline
-out_fid.createVariable('hc', 'f8', ('one'))
-out_fid.variables['hc'].long_name = 'S-coordinate parameter, critical depth'
-out_fid.variables['hc'].units = 'meter'
-out_fid.variables['hc'][:] = hc
-out_fid.createVariable('Cs_r', 'f8', ('s_rho'))
-out_fid.variables['Cs_r'].long_name = 'S-coordinate stretching curves at RHO-points'
-out_fid.variables['Cs_r'].units = 'nondimensional'
-out_fid.variables['Cs_r'].valid_min = -1.0
-out_fid.variables['Cs_r'].valid_max = 0.0
-out_fid.variables['Cs_r'][:] = Cs_r
-out_fid.createVariable('ocean_time', 'f8', ('ocean_time'))
-out_fid.variables['ocean_time'].long_name = 'time since initialization'
-out_fid.variables['ocean_time'].units = 'seconds'
-out_fid.variables['ocean_time'][0] = 0.0
-out_fid.createVariable('u', 'f8', ('ocean_time', 's_rho', 'eta_u', 'xi_u'))
-out_fid.variables['u'].long_name = 'u-momentum component'
-out_fid.variables['u'].units = 'meter second-1'
-out_fid.variables['u'][0,:,:,:] = u
-out_fid.createVariable('v', 'f8', ('ocean_time', 's_rho', 'eta_v', 'xi_v'))
-out_fid.variables['v'].long_name = 'v-momentum component'
-out_fid.variables['v'].units = 'meter second-1'
-out_fid.variables['v'][0,:,:,:] = v
-out_fid.createVariable('ubar', 'f8', ('ocean_time', 'eta_u', 'xi_u'))
-out_fid.variables['ubar'].long_name = 'vertically integrated u-momentum component'
-out_fid.variables['ubar'].units = 'meter second-1'
-out_fid.variables['ubar'][0,:,:] = ubar
-out_fid.createVariable('vbar', 'f8', ('ocean_time', 'eta_v', 'xi_v'))
-out_fid.variables['vbar'].long_name = 'vertically integrated v-momentum component'
-out_fid.variables['vbar'].units = 'meter second-1'
-out_fid.variables['vbar'][0,:,:] = vbar
-out_fid.createVariable('zeta', 'f8', ('ocean_time', 'eta_rho', 'xi_rho'))
-out_fid.variables['zeta'].long_name = 'free-surface'
-out_fid.variables['zeta'].units = 'meter'
-out_fid.variables['zeta'][0,:,:] = zeta
-out_fid.createVariable('temp', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
-out_fid.variables['temp'].long_name = 'potential temperature'
-out_fid.variables['temp'].units = 'Celsius'
-out_fid.variables['temp'][0,:,:,:] = temp
-out_fid.createVariable('salt', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
-out_fid.variables['salt'].long_name = 'salinity'
-out_fid.variables['salt'].units = 'PSU'
-out_fid.variables['salt'][0,:,:,:] = salt
-out_fid.createVariable('sc_r', 'f8', ('s_rho'))
-out_fid.variables['sc_r'].long_name = 'S-coordinate at rho-points'
-out_fid.variables['sc_r'].units = 'nondimensional'
-out_fid.variables['sc_r'].valid_min = -1.0
-out_fid.variables['sc_r'].valid_max = 0.0
-out_fid.variables['sc_r'][:] = sc_r
-out_fid.close()
-
-print 'Finished'
+    
 
         
     
diff --git a/romscice_ini_iceshelf.py b/romscice_ini_iceshelf.py
index 5d3bfee..ba26188 100644
--- a/romscice_ini_iceshelf.py
+++ b/romscice_ini_iceshelf.py
@@ -7,30 +7,38 @@
 # but it's more realistic than filling the cavities with constant salinity, and
 # certainly more stable).
 
-# Paths (other users need to edit these)
-ini_file = '../data/caisom001_ini_1995.nc' # Initial conditions file to edit
-is_ini_file = '../data/caisom001_ini_iceshelf.nc' # Ben's file
-grid_file = '../apps/common/grid/caisom001_OneQuartergrd.nc' # ROMS grid file
-
-# Read zice and make a boolean mask
-grid_id = Dataset(grid_file, 'r')
-mask_zice = grid_id.variables['mask_zice'][:,:]
-zice_index = mask_zice > 0.9 # better for floating-point than ==1
-grid_id.close()
-
-ini_id = Dataset(ini_file, 'a')
-is_id = Dataset(is_ini_file, 'r')
-num_depth = ini_id.variables['salt'].shape[1]
-
-for k in range(num_depth):
-    # For each depth level, replace the salinity values in ice shelf cavities
-    print 'Processing depth level ' + str(k+1) + ' of ' + str(num_depth)
-    ini_salt = ini_id.variables['salt'][0,k,:,:]
-    is_salt = is_id.variables['salt'][0,k,:,:]
-    ini_salt[zice_index] = is_salt[zice_index]
-    ini_id.variables['salt'][0,k,:,:] = ini_salt
-ini_id.close()
-is_id.close()
+def run(ini_file, is_ini_file, grid_file):
 
+    # Read zice and make a boolean mask
+    grid_id = Dataset(grid_file, 'r')
+    mask_zice = grid_id.variables['mask_zice'][:,:]
+    zice_index = mask_zice > 0.9 # better for floating-point than ==1
+    grid_id.close()
+
+    ini_id = Dataset(ini_file, 'a')
+    is_id = Dataset(is_ini_file, 'r')
+    num_depth = ini_id.variables['salt'].shape[1]
+
+    for k in range(num_depth):
+        # For each depth level, replace the salinity values in ice shelf cavities
+        print 'Processing depth level ' + str(k+1) + ' of ' + str(num_depth)
+        ini_salt = ini_id.variables['salt'][0,k,:,:]
+        is_salt = is_id.variables['salt'][0,k,:,:]
+        ini_salt[zice_index] = is_salt[zice_index]
+        ini_id.variables['salt'][0,k,:,:] = ini_salt
+    ini_id.close()
+    is_id.close()
+
+
+# User parameters
+if __name__ == "__main__":    
+
+    ini_file = '../data/caisom001_ini_1995.nc' # Initial conditions file to edit
+    is_ini_file = '../data/caisom001_ini_iceshelf.nc' # Ben's file
+    grid_file = '../apps/common/grid/caisom001_OneQuartergrd.nc' # ROMS grid file
+
+    run(ini_file, is_ini_file, grid_file)
+
+    
 
 
diff --git a/romscice_ini_woa.py b/romscice_ini_woa.py
index ae7171e..7219d0a 100644
--- a/romscice_ini_woa.py
+++ b/romscice_ini_woa.py
@@ -15,6 +15,139 @@
 from scipy.spatial import KDTree
 from calc_z import *
 
+
+# Main routine
+def run (grid_file, woa_file, output_file, Tcline, theta_s, theta_b, hc, N, nbdry_woa):
+
+    # Read WOA data and grid
+    print 'Reading World Ocean Atlas data'
+    woa_id = Dataset(woa_file, 'r')
+    lon_woa = woa_id.variables['longitude'][:]
+    lat_woa = woa_id.variables['latitude'][:nbdry_woa]
+    depth_woa = woa_id.variables['depth'][:]
+    temp_woa = transpose(woa_id.variables['temp'][:,:nbdry_woa,:])
+    salt_woa = transpose(woa_id.variables['salt'][:,:nbdry_woa,:])
+    woa_id.close()
+
+    # Read ROMS grid
+    print 'Reading ROMS grid'
+    grid_id = Dataset(grid_file, 'r')
+    lon_roms = grid_id.variables['lon_rho'][:,:]
+    lat_roms = grid_id.variables['lat_rho'][:,:]
+    h = grid_id.variables['h'][:,:]
+    zice = grid_id.variables['zice'][:,:]
+    mask_rho = grid_id.variables['mask_rho'][:,:]
+    mask_zice = grid_id.variables['mask_zice'][:,:]
+    grid_id.close()
+    num_lon = size(lon_roms, 1)
+    num_lat = size(lon_roms, 0)
+    # Mask h and zice with zeros
+    h = h*mask_rho
+    zice = zice*mask_zice
+
+    # Get 3D array of ROMS z-coordinates, as well as 1D arrays of s-coordinates
+    # and stretching curves
+    z_roms_3d, sc_r, Cs_r = calc_z(h, zice, lon_roms, lat_roms, theta_s, theta_b, hc, N)
+    # Copy the latitude and longitude values into 3D arrays of the same shape
+    lon_roms_3d = tile(lon_roms, (N,1,1))
+    lat_roms_3d = tile(lat_roms, (N,1,1))
+
+    # Regridding happens here
+    print 'Interpolating temperature'
+    temp = interp_woa2roms(temp_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, -0.5)
+    print 'Interpolating salinity'
+    salt = interp_woa2roms(salt_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, 34.5)
+
+    # Set initial velocities and sea surface height to zero
+    u = zeros((N, num_lat, num_lon-1))
+    v = zeros((N, num_lat-1, num_lon))
+    ubar = zeros((num_lat, num_lon-1))
+    vbar = zeros((num_lat-1, num_lon))
+    zeta = zeros((num_lat, num_lon))
+
+    print 'Writing to NetCDF file'
+    out_id = Dataset(output_file, 'w')
+    # Define dimensions
+    out_id.createDimension('xi_u', num_lon-1)
+    out_id.createDimension('xi_v', num_lon)
+    out_id.createDimension('xi_rho', num_lon)
+    out_id.createDimension('eta_u', num_lat)
+    out_id.createDimension('eta_v', num_lat-1)
+    out_id.createDimension('eta_rho', num_lat)
+    out_id.createDimension('s_rho', N)
+    out_id.createDimension('ocean_time', None)
+    out_id.createDimension('one', 1);
+    # Define variables and assign values
+    out_id.createVariable('tstart', 'f8', ('one'))
+    out_id.variables['tstart'].long_name = 'start processing day'
+    out_id.variables['tstart'].units = 'day'
+    out_id.variables['tstart'][:] = 0.0
+    out_id.createVariable('tend', 'f8', ('one'))
+    out_id.variables['tend'].long_name = 'end processing day'
+    out_id.variables['tend'].units = 'day'
+    out_id.variables['tend'][:] = 0.0
+    out_id.createVariable('theta_s', 'f8', ('one'))
+    out_id.variables['theta_s'].long_name = 'S-coordinate surface control parameter'
+    out_id.variables['theta_s'][:] = theta_s
+    out_id.createVariable('theta_b', 'f8', ('one'))
+    out_id.variables['theta_b'].long_name = 'S-coordinate bottom control parameter'
+    out_id.variables['theta_b'].units = 'nondimensional'
+    out_id.variables['theta_b'][:] = theta_b
+    out_id.createVariable('Tcline', 'f8', ('one'))
+    out_id.variables['Tcline'].long_name = 'S-coordinate surface/bottom layer width'
+    out_id.variables['Tcline'].units = 'meter'
+    out_id.variables['Tcline'][:] = Tcline
+    out_id.createVariable('hc', 'f8', ('one'))
+    out_id.variables['hc'].long_name = 'S-coordinate parameter, critical depth'
+    out_id.variables['hc'].units = 'meter'
+    out_id.variables['hc'][:] = hc
+    out_id.createVariable('Cs_r', 'f8', ('s_rho'))
+    out_id.variables['Cs_r'].long_name = 'S-coordinate stretching curves at RHO-points'
+    out_id.variables['Cs_r'].units = 'nondimensional'
+    out_id.variables['Cs_r'].valid_min = -1.0
+    out_id.variables['Cs_r'].valid_max = 0.0
+    out_id.variables['Cs_r'][:] = Cs_r
+    out_id.createVariable('ocean_time', 'f8', ('ocean_time'))
+    out_id.variables['ocean_time'].long_name = 'time since initialization'
+    out_id.variables['ocean_time'].units = 'seconds'
+    out_id.variables['ocean_time'][0] = 0.0
+    out_id.createVariable('u', 'f8', ('ocean_time', 's_rho', 'eta_u', 'xi_u'))
+    out_id.variables['u'].long_name = 'u-momentum component'
+    out_id.variables['u'].units = 'meter second-1'
+    out_id.variables['u'][0,:,:,:] = u
+    out_id.createVariable('v', 'f8', ('ocean_time', 's_rho', 'eta_v', 'xi_v'))
+    out_id.variables['v'].long_name = 'v-momentum component'
+    out_id.variables['v'].units = 'meter second-1'
+    out_id.variables['v'][0,:,:,:] = v
+    out_id.createVariable('ubar', 'f8', ('ocean_time', 'eta_u', 'xi_u'))
+    out_id.variables['ubar'].long_name = 'vertically integrated u-momentum component'
+    out_id.variables['ubar'].units = 'meter second-1'
+    out_id.variables['ubar'][0,:,:] = ubar
+    out_id.createVariable('vbar', 'f8', ('ocean_time', 'eta_v', 'xi_v'))
+    out_id.variables['vbar'].long_name = 'vertically integrated v-momentum component'
+    out_id.variables['vbar'].units = 'meter second-1'
+    out_id.variables['vbar'][0,:,:] = vbar
+    out_id.createVariable('zeta', 'f8', ('ocean_time', 'eta_rho', 'xi_rho'))
+    out_id.variables['zeta'].long_name = 'free-surface'
+    out_id.variables['zeta'].units = 'meter'
+    out_id.variables['zeta'][0,:,:] = zeta
+    out_id.createVariable('temp', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
+    out_id.variables['temp'].long_name = 'potential temperature'
+    out_id.variables['temp'].units = 'Celsius'
+    out_id.variables['temp'][0,:,:,:] = temp
+    out_id.createVariable('salt', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
+    out_id.variables['salt'].long_name = 'salinity'
+    out_id.variables['salt'].units = 'PSU'
+    out_id.variables['salt'][0,:,:,:] = salt
+    out_id.createVariable('sc_r', 'f8', ('s_rho'))
+    out_id.variables['sc_r'].long_name = 'S-coordinate at rho-points'
+    out_id.variables['sc_r'].units = 'nondimensional'
+    out_id.variables['sc_r'].valid_min = -1.0
+    out_id.variables['sc_r'].valid_max = 0.0
+    out_id.variables['sc_r'][:] = sc_r
+    out_id.close()
+
+
 # Given an array on the World Ocean Atlas grid, interpolate onto the ROMS grid,
 # fill the ice shelf cavities with nearest-neighbour values from the same depth
 # level, and fill the land mask with constant values.
@@ -78,150 +211,28 @@ def interp_woa2roms (A, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z
     return B
 
 
-# File paths to edit here:
-# Path to ROMS grid file
-grid_file = '../apps/common/grid/caisom001_OneQuartergrd.nc'
-# Path to World Ocean Atlas NetCDF file (converted from FESOM input using 
-# woa_netcdf.py)
-woa_file = '/short/y99/kaa561/FESOM/woa01_ts.nc'
-# Path to desired output file
-output_file = '../data/caisom001_ini.nc'
-
-# Grid parameters: check grid_file and *.in file to make sure these are correct
-Tcline = 40
-theta_s = 0.9
-theta_b = 4
-hc = 40
-N = 31
-
-# Northernmost index of WOA grid to read (1-based)
-nbdry_woa = 32
-
-# Read WOA data and grid
-print 'Reading World Ocean Atlas data'
-woa_id = Dataset(woa_file, 'r')
-lon_woa = woa_id.variables['longitude'][:]
-lat_woa = woa_id.variables['latitude'][:nbdry_woa]
-depth_woa = woa_id.variables['depth'][:]
-temp_woa = transpose(woa_id.variables['temp'][:,:nbdry_woa,:])
-salt_woa = transpose(woa_id.variables['salt'][:,:nbdry_woa,:])
-woa_id.close()
-
-# Read ROMS grid
-print 'Reading ROMS grid'
-grid_id = Dataset(grid_file, 'r')
-lon_roms = grid_id.variables['lon_rho'][:,:]
-lat_roms = grid_id.variables['lat_rho'][:,:]
-h = grid_id.variables['h'][:,:]
-zice = grid_id.variables['zice'][:,:]
-mask_rho = grid_id.variables['mask_rho'][:,:]
-mask_zice = grid_id.variables['mask_zice'][:,:]
-grid_id.close()
-num_lon = size(lon_roms, 1)
-num_lat = size(lon_roms, 0)
-# Mask h and zice with zeros
-h = h*mask_rho
-zice = zice*mask_zice
-
-# Get 3D array of ROMS z-coordinates, as well as 1D arrays of s-coordinates
-# and stretching curves
-z_roms_3d, sc_r, Cs_r = calc_z(h, zice, lon_roms, lat_roms, theta_s, theta_b, hc, N)
-# Copy the latitude and longitude values into 3D arrays of the same shape
-lon_roms_3d = tile(lon_roms, (N,1,1))
-lat_roms_3d = tile(lat_roms, (N,1,1))
-
-# Regridding happens here
-print 'Interpolating temperature'
-temp = interp_woa2roms(temp_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, -0.5)
-print 'Interpolating salinity'
-salt = interp_woa2roms(salt_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, 34.5)
-
-# Set initial velocities and sea surface height to zero
-u = zeros((N, num_lat, num_lon-1))
-v = zeros((N, num_lat-1, num_lon))
-ubar = zeros((num_lat, num_lon-1))
-vbar = zeros((num_lat-1, num_lon))
-zeta = zeros((num_lat, num_lon))
-
-print 'Writing to NetCDF file'
-out_id = Dataset(output_file, 'w')
-# Define dimensions
-out_id.createDimension('xi_u', num_lon-1)
-out_id.createDimension('xi_v', num_lon)
-out_id.createDimension('xi_rho', num_lon)
-out_id.createDimension('eta_u', num_lat)
-out_id.createDimension('eta_v', num_lat-1)
-out_id.createDimension('eta_rho', num_lat)
-out_id.createDimension('s_rho', N)
-out_id.createDimension('ocean_time', None)
-out_id.createDimension('one', 1);
-# Define variables and assign values
-out_id.createVariable('tstart', 'f8', ('one'))
-out_id.variables['tstart'].long_name = 'start processing day'
-out_id.variables['tstart'].units = 'day'
-out_id.variables['tstart'][:] = 0.0
-out_id.createVariable('tend', 'f8', ('one'))
-out_id.variables['tend'].long_name = 'end processing day'
-out_id.variables['tend'].units = 'day'
-out_id.variables['tend'][:] = 0.0
-out_id.createVariable('theta_s', 'f8', ('one'))
-out_id.variables['theta_s'].long_name = 'S-coordinate surface control parameter'
-out_id.variables['theta_s'][:] = theta_s
-out_id.createVariable('theta_b', 'f8', ('one'))
-out_id.variables['theta_b'].long_name = 'S-coordinate bottom control parameter'
-out_id.variables['theta_b'].units = 'nondimensional'
-out_id.variables['theta_b'][:] = theta_b
-out_id.createVariable('Tcline', 'f8', ('one'))
-out_id.variables['Tcline'].long_name = 'S-coordinate surface/bottom layer width'
-out_id.variables['Tcline'].units = 'meter'
-out_id.variables['Tcline'][:] = Tcline
-out_id.createVariable('hc', 'f8', ('one'))
-out_id.variables['hc'].long_name = 'S-coordinate parameter, critical depth'
-out_id.variables['hc'].units = 'meter'
-out_id.variables['hc'][:] = hc
-out_id.createVariable('Cs_r', 'f8', ('s_rho'))
-out_id.variables['Cs_r'].long_name = 'S-coordinate stretching curves at RHO-points'
-out_id.variables['Cs_r'].units = 'nondimensional'
-out_id.variables['Cs_r'].valid_min = -1.0
-out_id.variables['Cs_r'].valid_max = 0.0
-out_id.variables['Cs_r'][:] = Cs_r
-out_id.createVariable('ocean_time', 'f8', ('ocean_time'))
-out_id.variables['ocean_time'].long_name = 'time since initialization'
-out_id.variables['ocean_time'].units = 'seconds'
-out_id.variables['ocean_time'][0] = 0.0
-out_id.createVariable('u', 'f8', ('ocean_time', 's_rho', 'eta_u', 'xi_u'))
-out_id.variables['u'].long_name = 'u-momentum component'
-out_id.variables['u'].units = 'meter second-1'
-out_id.variables['u'][0,:,:,:] = u
-out_id.createVariable('v', 'f8', ('ocean_time', 's_rho', 'eta_v', 'xi_v'))
-out_id.variables['v'].long_name = 'v-momentum component'
-out_id.variables['v'].units = 'meter second-1'
-out_id.variables['v'][0,:,:,:] = v
-out_id.createVariable('ubar', 'f8', ('ocean_time', 'eta_u', 'xi_u'))
-out_id.variables['ubar'].long_name = 'vertically integrated u-momentum component'
-out_id.variables['ubar'].units = 'meter second-1'
-out_id.variables['ubar'][0,:,:] = ubar
-out_id.createVariable('vbar', 'f8', ('ocean_time', 'eta_v', 'xi_v'))
-out_id.variables['vbar'].long_name = 'vertically integrated v-momentum component'
-out_id.variables['vbar'].units = 'meter second-1'
-out_id.variables['vbar'][0,:,:] = vbar
-out_id.createVariable('zeta', 'f8', ('ocean_time', 'eta_rho', 'xi_rho'))
-out_id.variables['zeta'].long_name = 'free-surface'
-out_id.variables['zeta'].units = 'meter'
-out_id.variables['zeta'][0,:,:] = zeta
-out_id.createVariable('temp', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
-out_id.variables['temp'].long_name = 'potential temperature'
-out_id.variables['temp'].units = 'Celsius'
-out_id.variables['temp'][0,:,:,:] = temp
-out_id.createVariable('salt', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'))
-out_id.variables['salt'].long_name = 'salinity'
-out_id.variables['salt'].units = 'PSU'
-out_id.variables['salt'][0,:,:,:] = salt
-out_id.createVariable('sc_r', 'f8', ('s_rho'))
-out_id.variables['sc_r'].long_name = 'S-coordinate at rho-points'
-out_id.variables['sc_r'].units = 'nondimensional'
-out_id.variables['sc_r'].valid_min = -1.0
-out_id.variables['sc_r'].valid_max = 0.0
-out_id.variables['sc_r'][:] = sc_r
-out_id.close()
+# User parameters
+if __name__ == "__main__":
+
+    # Path to ROMS grid file
+    grid_file = '../apps/common/grid/caisom001_OneQuartergrd.nc'
+    # Path to World Ocean Atlas NetCDF file (converted from FESOM input using 
+    # woa_netcdf.py)
+    woa_file = '/short/y99/kaa561/FESOM/woa01_ts.nc'
+    # Path to desired output file
+    output_file = '../data/caisom001_ini.nc'
+
+    # Grid parameters: check grid_file and *.in file to make sure these are correct
+    Tcline = 40
+    theta_s = 0.9
+    theta_b = 4
+    hc = 40
+    N = 31
+
+    # Northernmost index of WOA grid to read (1-based)
+    nbdry_woa = 32
+
+    run(grid_file, woa_file, output_file, Tcline, theta_s, theta_b, hc, N, nbdry_woa)
+
+    
 
diff --git a/woa_netcdf.py b/woa_netcdf.py
index 0be6046..8599d8f 100644
--- a/woa_netcdf.py
+++ b/woa_netcdf.py
@@ -1,100 +1,102 @@
 from netCDF4 import Dataset
 from numpy import *
 
-# Read the World Ocean Atlas temperature and salinity climatology from a text
-# file (FESOM input format), and save in a NetCDF file for easier use.
+if __name__ == "__main__":
 
-# File paths
-in_file = '/short/y99/kaa561/FESOM/annual_woa01_ts.out'
-out_file = '/short/y99/kaa561/FEsOM/woa01_ts.nc'
+    # Read the World Ocean Atlas temperature and salinity climatology from a text
+    # file (FESOM input format), and save in a NetCDF file for easier use.
 
-print 'Reading text file'
-file = open(in_file, 'r')
+    # File paths
+    in_file = '/short/y99/kaa561/FESOM/annual_woa01_ts.out'
+    out_file = '/short/y99/kaa561/FEsOM/woa01_ts.nc'
 
-# First line contains the number of longitude, latitude, and depth values
-sizes = (file.readline()).split()
-num_lon = int(sizes[0])
-num_lat = int(sizes[1])
-num_depth = int(sizes[2])
+    print 'Reading text file'
+    file = open(in_file, 'r')
 
-# Read longitude values
-lon = []
-while True:
-    # Each row has multiple values, so split based on white space
-    lon_vals = (file.readline()).split()
-    for elm in lon_vals:
-        lon.append(float(elm))
-    if len(lon) == num_lon:
-        # Finished with longitude
-        break
-# Repeat for latitude
-lat = []
-while True:
-    lat_vals = (file.readline()).split()
-    for elm in lat_vals:
-        lat.append(float(elm))
-    if len(lat) == num_lat:
-        break
-# Repeat for depth
-depth = []
-while True:
-    depth_vals = (file.readline()).split()
-    for elm in depth_vals:
-        # Convert from negative to positive depth
-        depth.append(-1*float(elm))
-    if len(depth) == num_depth:
-        break
-# Read temperature values (save in one long 1D array for now)
-temp = []
-while True:
-    temp_vals = (file.readline()).split()
-    for elm in temp_vals:
-        temp.append(float(elm))
-    if len(temp) == num_lon*num_lat*num_depth:
-        break
-# Repeat for salinity
-salt = []
-while True:
-    salt_vals = (file.readline()).split()
-    for elm in salt_vals:
-        salt.append(float(elm))
-    if len(salt) == num_lon*num_lat*num_depth:
-        break
-file.close()
+    # First line contains the number of longitude, latitude, and depth values
+    sizes = (file.readline()).split()
+    num_lon = int(sizes[0])
+    num_lat = int(sizes[1])
+    num_depth = int(sizes[2])
 
-# Copy contents of the long 1D temp and salt arrays into 3D arrays
-# (depth x latitude x longitude)
-print 'Reshaping temperature and salinity arrays'
-temp_3d = zeros((num_depth, num_lat, num_lon))
-salt_3d = zeros((num_depth, num_lat, num_lon))
-posn = 0
-for i in range(num_lon):
-    for j in range(num_lat):
-        for k in range(num_depth):
-            temp_3d[k,j,i] = temp[posn]
-            salt_3d[k,j,i] = salt[posn]
-            posn = posn+1
+    # Read longitude values
+    lon = []
+    while True:
+        # Each row has multiple values, so split based on white space
+        lon_vals = (file.readline()).split()
+        for elm in lon_vals:
+            lon.append(float(elm))
+        if len(lon) == num_lon:
+            # Finished with longitude
+            break
+    # Repeat for latitude
+    lat = []
+    while True:
+        lat_vals = (file.readline()).split()
+        for elm in lat_vals:
+            lat.append(float(elm))
+        if len(lat) == num_lat:
+            break
+    # Repeat for depth
+    depth = []
+    while True:
+        depth_vals = (file.readline()).split()
+        for elm in depth_vals:
+            # Convert from negative to positive depth
+            depth.append(-1*float(elm))
+        if len(depth) == num_depth:
+            break
+    # Read temperature values (save in one long 1D array for now)
+    temp = []
+    while True:
+        temp_vals = (file.readline()).split()
+        for elm in temp_vals:
+            temp.append(float(elm))
+        if len(temp) == num_lon*num_lat*num_depth:
+            break
+    # Repeat for salinity
+    salt = []
+    while True:
+        salt_vals = (file.readline()).split()
+        for elm in salt_vals:
+            salt.append(float(elm))
+        if len(salt) == num_lon*num_lat*num_depth:
+            break
+    file.close()
 
-# Output to NetCDF file
-print 'Writing NetCDF file'
-id = Dataset(out_file, 'w')
-id.createDimension('longitude', num_lon)
-id.createDimension('latitude', num_lat)
-id.createDimension('depth', num_depth)
-id.createVariable('longitude', 'f8', ('longitude'))
-id.variables['longitude'].units = 'degrees'
-id.variables['longitude'][:] = lon
-id.createVariable('latitude', 'f8', ('latitude'))
-id.variables['latitude'].units = 'degrees'
-id.variables['latitude'][:] = lat
-id.createVariable('depth', 'f8', ('depth'))
-id.variables['depth'].units = 'metres'
-id.variables['depth'][:] = depth
-id.createVariable('temp', 'f8', ('depth', 'latitude', 'longitude'))
-id.variables['temp'].units = 'C'
-id.variables['temp'][:,:,:] = temp_3d
-id.createVariable('salt', 'f8', ('depth', 'latitude', 'longitude'))
-id.variables['salt'].units = 'psu'
-id.variables['salt'][:,:,:] = salt_3d
-id.close()
+    # Copy contents of the long 1D temp and salt arrays into 3D arrays
+    # (depth x latitude x longitude)
+    print 'Reshaping temperature and salinity arrays'
+    temp_3d = zeros((num_depth, num_lat, num_lon))
+    salt_3d = zeros((num_depth, num_lat, num_lon))
+    posn = 0
+    for i in range(num_lon):
+        for j in range(num_lat):
+            for k in range(num_depth):
+                temp_3d[k,j,i] = temp[posn]
+                salt_3d[k,j,i] = salt[posn]
+                posn = posn+1
+
+    # Output to NetCDF file
+    print 'Writing NetCDF file'
+    id = Dataset(out_file, 'w')
+    id.createDimension('longitude', num_lon)
+    id.createDimension('latitude', num_lat)
+    id.createDimension('depth', num_depth)
+    id.createVariable('longitude', 'f8', ('longitude'))
+    id.variables['longitude'].units = 'degrees'
+    id.variables['longitude'][:] = lon
+    id.createVariable('latitude', 'f8', ('latitude'))
+    id.variables['latitude'].units = 'degrees'
+    id.variables['latitude'][:] = lat
+    id.createVariable('depth', 'f8', ('depth'))
+    id.variables['depth'].units = 'metres'
+    id.variables['depth'][:] = depth
+    id.createVariable('temp', 'f8', ('depth', 'latitude', 'longitude'))
+    id.variables['temp'].units = 'C'
+    id.variables['temp'][:,:,:] = temp_3d
+    id.createVariable('salt', 'f8', ('depth', 'latitude', 'longitude'))
+    id.variables['salt'].units = 'psu'
+    id.variables['salt'][:,:,:] = salt_3d
+    id.close()