<p><b>mhoffman@lanl.gov</b> 2013-03-08 16:01:58 -0700 (Fri, 08 Mar 2013)</p><p>BRANCH COMMIT - ice-ocean coupler<br>
<br>
Initial work on a design document and a script for doing land ice-ocean coupling. <br>
</p><hr noshade><pre><font color="gray">Added: branches/landice_ocean_coupler/coupler-time.png
===================================================================
(Binary files differ)
Index: branches/landice_ocean_coupler/coupler-time.png
===================================================================
--- branches/landice_ocean_coupler/coupler-time.png 2013-03-08 22:58:18 UTC (rev 2576)
+++ branches/landice_ocean_coupler/coupler-time.png 2013-03-08 23:01:58 UTC (rev 2577)
Property changes on: branches/landice_ocean_coupler/coupler-time.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property
Added: branches/landice_ocean_coupler/coupler.png
===================================================================
(Binary files differ)
Index: branches/landice_ocean_coupler/coupler.png
===================================================================
--- branches/landice_ocean_coupler/coupler.png 2013-03-08 22:58:18 UTC (rev 2576)
+++ branches/landice_ocean_coupler/coupler.png 2013-03-08 23:01:58 UTC (rev 2577)
Property changes on: branches/landice_ocean_coupler/coupler.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property
Added: branches/landice_ocean_coupler/ice-ocean-coupling-DesignDocument.tex
===================================================================
--- branches/landice_ocean_coupler/ice-ocean-coupling-DesignDocument.tex (rev 0)
+++ branches/landice_ocean_coupler/ice-ocean-coupling-DesignDocument.tex 2013-03-08 23:01:58 UTC (rev 2577)
@@ -0,0 +1,199 @@
+ \documentclass[11pt]{report}
+
+\usepackage{epsf,amsmath,amsfonts}
+\usepackage{graphicx}
+
+\usepackage{fullpage}
+\usepackage{mathtools}
+\usepackage{listings}
+\usepackage{float}
+
+</font>
<font color="gray">ewcommand{\vect}[1]{\mathbf{#1}}
+\restylefloat{figure}
+
+\begin{document}
+
+\title{ MPAS Ice-Ocean Offline Coupling: \\
+Requirements and Design}
+\author{M. Hoffman, D. Jacobsen}
+
+\maketitle
+\tableofcontents
+
+%-----------------------------------------------------------------------
+
+\chapter{Summary}
+
+We present the design for a coupler between the MPAS Ocean and Land Ice models that runs the two models independently for short intervals, modifying the geometry and forcing for each model.
+
+
+%-----------------------------------------------------------------------
+
+\chapter{Requirements}
+
+\section{Requirement: Use existing models with minimal modifications.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+\section{Requirement: The coupler calculates boundary-layer physics}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+Ice sheet model receives updated basal mass balance and basal heat flux.
+
+Ocean model receives updated ocean surface pressure and mass, heat, and salinity fluxes.
+
+\section{Requirement: The two models can run with different time steps.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+%-----------------------------------------------------------------------
+
+\chapter{Algorithmic Formulations}
+
+\section{Design Solution: Use existing models with minimal modifications.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+Using MPAS' restart capability, write a script in Python or Bash that calls each model for set periods of time, recalculates forcings for the two models, and then restarts each.
+
+\section{Design Solution: The coupler calculates boundary-layer physics}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+
+
+\begin{figure}
+ \center{\includegraphics[width=8cm]{./coupler.png}}
+ \caption{Diagram of information the coupler needs to calculate (black circles) and pass to each model (arrows).}
+ \label{coupler}
+\end{figure}
+
+\subsection{Long term}
+ Apply boundary-layer physics from ISOMIP experiments (Hunter 2006, Losch 2008, Holland and Jenkins 1999).
+
+\subsection{Short term}
+Simplified calculation of forcings:
+
+Ice Sheet Forcings
+\begin{itemize}
+\item
+BMB = linear function of ocean temperature
+\item
+BHF not initially needed since MPAS land ice currently does not perform column temperature diffusion, which is the calculation that this would affect.
+\end{itemize}
+
+Ocean Forcings
+\begin{itemize}
+\item
+ocean surface pressure
+\begin{equation}
+ \label{pressure}
+ P_{os} = \rho_i g H
+\end{equation}
+where $P_{os}$ is the ocean surface pressure (Pa), $\rho_i$ is ice density (kg m$^{-3}$), and $H$ is ice thickness (m).
+\item
+ocean surface mass flux = BMB
+
+\item
+ocean surface heat flux
+
+Assume ice at 0$^{\circ}$ C.
+
+\item
+ocean surface salinity flux
+
+Assume melt is fresh. What about refreezing?
+
+\end{itemize}
+
+
+\section{Design Solution: The two models can run with different time steps.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+For some applications we may want to run both models at the same time step and for others we will want to allow the ice sheet model to have a longer time step for efficiency. The ice sheet model typically runs with time steps of $\sim$0.1-1.0 years. The ocean model typically runs with time steps of $\sim$1 hr.
+
+At each coupling, the coupler will calculate $P_{os}$ and the various fluxes using the current state of both models. (An alternative would be to make the fluxes based on some backward-looking average of the state in the ocean model.) These fields will then be held constant within each model until the next coupling. Because the fluxes are specified as rates, no additional modification needs to be done. This also ensures that mass and energy will be conserved since both models are using the same fluxes over the same total time. We impose the requirement that the ocean time step be evenly divisible into the land ice time step.
+
+
+\begin{figure}[hbt]
+ \center{\includegraphics[width=12cm]{./coupler-time.png}}
+ \caption{Diagram of typical time steps for the two models.}
+ \label{coupler-time}
+\end{figure}
+
+%-----------------------------------------------------------------------
+
+\chapter{Design and Implementation}
+
+
+\section{Implementation: Use existing models with minimal modifications.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+We will use Python to write the coupler because it is cross-platform, allows both system calls (running the model) and calculations (netCDF I/O, matrix math), and is familiar to the developers (MJH, DJ).
+
+Each model will need to be setup (compiled, input file created, namelist file created) separately. Initially we will require that both models use the same grid, though we could explore relaxing that requirement later. The coupler does not need to know the models' time steps (unless the flux calculations are based on a time-average from the ocean model). The order in which the two models are run does not matter, and we could explore running them in parallel.
+
+\begin{lstlisting}[numbers=left, numberstyle=\tiny, basicstyle=\small, stringstyle=\ttfamily, mathescape, escapeinside={@}{@}]
+Initialize - identify executables and input/output files for each model
+Read ice sheet input file for thickness and basal temperature initial state
+Read ocean model input file for surface temperature and salinity initial state
+Calculate @$P_{os}$@
+Calculate boundary layer fluxes
+Write BMB, BHF to ice sheet model input file on first call
+Write @$P_{os}$@, fluxes to ocean model input file
+Run ocean model
+Run ice sheet model
+for t in range(coupleIntervals): # How do we want to keep time?
+ Copy ice sheet and ocean output files
+ Read ice sheet output file for thickness and basal temperature
+ Read ocean model output file for surface temperature and salinity
+ Calculate @$P_{os}$@
+ Calculate boundary layer fluxes
+ Write BMB, BHF to ice sheet model restart file
+ Write @$P_{os}$@ to ocean model restart file
+ Run ocean model
+ Run ice sheet model
+\end{lstlisting}
+
+
+\section{Implementation: The coupler calculates boundary-layer physics}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+It will be cleaner to put these calculations in Python functions are a separate python module so the simplified calculations can be swapped out for the more sophisticated boundary-layer physics calculations.
+
+\section{Implementation: The two models can run with different time steps.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+This is implemented by specifying different time steps in the namelist files for each model.
+
+
+%-----------------------------------------------------------------------
+
+\chapter{Testing}
+
+\section{Testing: Use existing models with minimal modifications.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+\section{Testing: The coupler calculates boundary-layer physics}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+Run ISOMIP through the coupler.
+
+\section{Testing: The two models can run with different time steps.}
+Date last modified: 2013/03/08 \\
+Contributors: MJH \\
+
+
+\begin{itemize}
+\item
+\end{itemize}
+%-----------------------------------------------------------------------
+
+\end{document}
Added: branches/landice_ocean_coupler/landice-ocean-coupler.py
===================================================================
--- branches/landice_ocean_coupler/landice-ocean-coupler.py (rev 0)
+++ branches/landice_ocean_coupler/landice-ocean-coupler.py 2013-03-08 23:01:58 UTC (rev 2577)
@@ -0,0 +1,173 @@
+#!/usr/bin/python
+# Coupler for the MPAS Ocean and Land Ice models
+# Matt Hoffman, March 8, 2013
+
+import numpy
+import os, sys, subprocess
+#import matplotlib.pyplot as plt
+try:
+ from Scientific.IO.NetCDF import NetCDFFile
+ netCDF_module = 'Scientific.IO.NetCDF'
+except ImportError:
+ try:
+ from netCDF4 import Dataset as NetCDFFile
+ netCDF_module = 'netCDF4'
+ except ImportError:
+ print 'Unable to import any of the following python modules:'
+ print ' Scientific.IO.NetCDF </font>
<font color="blue"> netcdf4 '
+ print 'One of them must be installed.'
+ raise ImportError('No netCDF module found')
+
+
+# In the future we may want to add command line arguments using the optParse module.
+
+# The script assumes you have setup a directory structure:
+# ocean model is setup in ./ocean
+# land ice model is setup in ./land_ice
+
+# =======================================
+# Setup needed files
+ocnInput = 'grid.nc'
+ocnOutput = 'output.nc'
+
+liceInput = 'land_ice_grid.nc'
+liceOutput = 'output.nc'
+
+# Setup needed executables
+runLICE = 'mpirun -np 4 land_ice_model.exe'
+runOCN = 'mpirun -np 4 ocean_model.exe'
+
+# Define number of coupling intervals (assumes you have set dt for each model accordingly, and ocean dt is evenly divisble in land ice dt)
+# (We may want to do more sophisticated time-keeping of the coupling if there is a need.)
+numCoupleIntervals = 3
+
+# Some parameters
+rhoi = 900.0 # land ice model's ice density (kg m^-3)
+grav = 9.8101 # gravitational acceleration
+
+# =======================================
+
+
+
+
+for t in range(coupleIntervals):
+ print 'Starting coupling number ', t
+
+ # Get the appropriate source and destination files for the input needed by the coupler
+ if t==0:
+ # Use I.C. from input files at initial time
+ liceSourceFile = NetCDFFile('./land_ice/'+liceInput,'r')
+ ocnSourceFile = NetCDFFile('./ocean/'+ocnInput,'r')
+
+ # At initial time, we want to write to the input files
+ liceDestFile = NetCDFFile('./land_ice/'+liceInput,'r+')
+ ocnDestFile = NetCDFFile('./ocean/'+ocnInput,'r+')
+
+ else:
+ # Use values from output files at other times.
+ liceSourceFile = NetCDFFile('./land_ice/'+liceOutput,'r')
+ ocnSourceFile = NetCDFFile('./ocean/'+ocnOutput,'r')
+
+ # At other times we want to write to the restart files
+ # Need to choose the appropriate restart file - there may be smarter ways to do this, but an easy way is to take the most recent
+ # os.listdir does not have a sort ability, so use a shell command (unpythonic but works)
+ p = os.popen('ls -1t ./land_ice/restart.*.nc | head -n 1')
+ liceRestart = p.read().rstrip() # get stdout
+ liceDestFile = NetCDFFile(liceRestart,'r+')
+ p = os.popen('ls -1t ./ocean/restart.*.nc | head -n 1')
+ ocnRestart = p.read().rstrip() # get stdout
+ ocnDestFile = NetCDFFile(ocnRestart,'r+')
+
+
+ # Get the source fields needed by the coupler
+ # (Get the -1 time level for all variables - this will be the last whether this is from input or output
+ liceThickness = liceSourceFile.variables['thickness'][-1,:]
+ liceBasalTemp = liceSourceFile.variables['temperature'][-1,:,-1]
+ ocnSurfTemp = ocnSourceFile.variables['temperature'][-1,:,0]
+ ocnSurfSalin = ocnSourceFile.variables['salinity'][-1,:,0]
+
+ # Get the destination fields needed by the coupler
+ liceBMB = liceDestFile.variables['marineBasalMassBalanceTimeSeries'][:,0] # There should only be on time level
+ #liceBHF = liceDestFile.variables['basalHeatFluxTimeSeries'][:,0] # NEED TO RESTRICT THIS TO CELLS WITH ICE
+ ocnSurfPressure = ocnDestFile.variables['ssh'][0,:] # TODO WHAT IS THE PROPER VARIABLE HERE?
+ # TODO WHAT ARE THE PROPER VARIABLE FOR MASS, HEAT, AND SALINITY FLUXES?
+
+ # Calculate the ocean surface pressure - this could be a function call
+ # TODO WHAT TO DO ABOUT NON-ICE SHELF OCEAN CELLS?
+ ocnSurfPressure = rhoi * grav * liceThickness # TODO ARE THESE THE PROPER UNITS?
+
+ # Calculate boundary layer fluxes - this could be an internal or external function call
+ #calculateBdyLyrFluxes(liceBasalTemp, ocnSurfTemp, ocnSurfSalin, liceBHF, ocnMassFlux, ocnHeatFlux, ocnSalinFlux) #TODO
+
+ # Write new quantities to the destination files
+ # TODO Does this need to be done explicitly? I always forget how netCDF i/o works in python precisely.
+
+ # =======================================
+ # Run the models # TODO Do we want to do this in parallel?
+ try:
+ os.chdir('ocean')
+ # Set the namelist file properly for a restart or no
+ print 'Setting up ocean namelist file'
+ if t==0:
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_do_restart.*/ config_do_restart = .false./" namelist.input" ' , shell=True, executable='/bin/bash') # no restart
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_write_output_on_startup.*/ config_write_output_on_startup = .true./" namelist.input" ' , shell=True, executable='/bin/bash') # write the initial time
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_start_time.*/ config_start_time = '0000-01-01_00:00:00'/" namelist.input" ' , shell=True, executable='/bin/bash') # set the start time
+ else:
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_do_restart.*/ config_do_restart = .true./" namelist.input" ' , shell=True, executable='/bin/bash') # do a restart
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_write_output_on_startup.*/ config_write_output_on_startup = .false./" namelist.input" ' , shell=True, executable='/bin/bash') # don't write the initial time
+ # Get the old stop time so we can set it to be the new start time.
+ p = os.popen('ls -1t output*nc | head -n 1')
+ oldoutputfile=p.read().rstrip()
+ p = os.popen('ncks -v xtime -H -s "%c" ' + oldoutputfile + ' | tail -r -c 65')
+ oldstoptime=p.read().rstrip()
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_start_time.*/ config_start_time = "' + oldstoptime '"/" namelist.input" ' , shell=True, executable='/bin/bash') # set the start time to the old stop time
+
+ print 'Starting ocean model.'
+ subprocess.check_call(runOCN, shell=True, executable='/bin/bash')
+ print 'Ocean model completed.'
+ # Copy output file so we don't lose it when the model restarts - could use Python's shutil module. TODO Also may want to change format of the string conversion of t
+ subprocess.check_call('cp ' + ocnOutput + ' output.' + str(t) + '.nc', shell=True, executable='/bin/bash')
+ os.chdir('..')
+ except:
+ sys.exit('Ocean model failed!')
+
+ try:
+ os.chdir('land_ice')
+ # Set the namelist file properly for a restart or no and set the start and stop times
+ print 'Setting up land ice namelist file'
+ if t==0:
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_do_restart.*/ config_do_restart = .false./" namelist.input" ' , shell=True, executable='/bin/bash') # no restart
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_write_output_on_startup.*/ config_write_output_on_startup = .true./" namelist.input" ' , shell=True, executable='/bin/bash') # write the initial time
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_start_time.*/ config_start_time = '0000-01-01_00:00:00'/" namelist.input" ' , shell=True, executable='/bin/bash') # set the start time
+ else:
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_do_restart.*/ config_do_restart = .true./" namelist.input" ' , shell=True, executable='/bin/bash') # do a restart
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_write_output_on_startup.*/ config_write_output_on_startup = .false./" namelist.input" ' , shell=True, executable='/bin/bash') # don't write the initial time
+ # Get the old stop time so we can set it to be the new start time.
+ p = os.popen('ls -1t output*nc | head -n 1')
+ oldoutputfile=p.read().rstrip()
+ p = os.popen('ncks -v xtime -H -s "%c" ' + oldoutputfile + ' | tail -r -c 65')
+ oldstoptime=p.read().rstrip()
+ subprocess.check_call('sed -i.SEDBACKUP "s/^.*config_start_time.*/ config_start_time = "' + oldstoptime '"/" namelist.input" ' , shell=True, executable='/bin/bash') # set the start time to the old stop time
+
+ print 'Starting land ice model.'
+ subprocess.check_call(runLICE, shell=True, executable='/bin/bash')
+ print 'Land ice model completed.'
+ # Copy output file so we don't lose it when the model restarts - could use Python's shutil module. TODO Also may want to change format of the string conversion of t
+ subprocess.check_call('cp ' + liceOutput + ' output.' + str(t) + '.nc', shell=True, executable='/bin/bash')
+ os.chdir('..')
+ except:
+ sys.exit('Land Ice model failed!')
+
+# =======================================
+
+print 'Coupled model run complete!'
+
+# After the coupled run is completed, we may want to go in and use NCO to patch together all of the output files for each model.
+
+
+
+
+
+
+
+
</font>
</pre>