CSTMS UNSW2008

Adding Sediment to the Latte Example

On Friday, 3 April 2009, Chris Sherwood presented a tutorial for new ROMS users interested in sediments. The tutorial was held at the Sydney Institute of Marine Sciences as part of the ROMS/TOMS Asia-Pacific Workshop. This tutorial continues using the LATTE example presented during the Monday tutorial.

UNDER CONSTRUCTION - I'll post a note on the forum when this provides a consistent, step-by-step approach. Until then, user beware.

Set up and test the LATTE case

These instructions differ a little from the Monday tutorial because we assume you are working on your own computer.

• If you have not already done so, check out a version of the code. We'll assume you keep it in a directory called ~/src. If you want to get the exact version used in this tutorial, then cd to the parent directory for your source code and then:

mkdir src
svn checkout -r 342 --username bruce https://www.myroms.org/svn/src/trunk src

replacing bruce with your ROMS user name. Omitting the -r 342 will get you the most current (HEAD) version of the code.

• Make sure you can compile and run UPWELLING.

• Make a directory structure that looks like this:

../Latte
../Latte/original
../Latte/in
../Latte/out
../Latte/run01

We will add a few more directories later. Files for this case are at http://coast-enviro.er.usgs.gov/models/Latte/. Use your browser to copy the contents of original and in.

Build and run the Latte example

• Copy of the original .h and .in files into the ../Latte/Run01 directory.
• Try to build and run the Latte case.
  • If you can use the build.bash script, set the correct entries for environment variables:
    ROMS_APPLICATION=LATTE_C causes build.bash to look for the file latte_c.h in order to set the CPP options
    MY_PROJECT_DIR=${HOME}/Projects/latte_c will instruct build.bash where to look for the latte_c.h file.
  • No one has figured out how to use build.bash with the Intel (ifort) compiler under Cygwin on a Windows machine. If this is your set up, then edit the makefile in the src directory so that it finds and builds the application.
    ROMS_APPLICATION ?= LATTE_C
    MY_HEADER_DIR ?= /cygdrive/d/crs/proj/CSTMS_Tutorial/Latte/run01
  • If you get a complaint about a missing stations file, you can either #UNDEF STATIONS in the .h file (and rebuild), or copy a valid stations file to the local (Run01) directory, or change the filename mentioned in the .in file to point to a valid stations file. Unless you change the path of the output file, it will be end up in stations.nc in the run directory, while other output goes to ../out.
  • If you get an error that T_cline or S_cline don't match the input file, follow the instructions here Changing netCDF variables.

Change the advection scheme (Run01)

• You need to use a "positive-definite" advection scheme with sediment. The small over- and undershoots that occur with many advection algorithms will cause mass-conservation problems with sediment. (What, for example, is deposition of negative sediment concentrations? Erosion?). TS_MPDATA is a Multidimensional Positive Definite Advection Transport Algorithm method based on the work of Piotr Smolarkiewicz. See the [[1]wiki discussion.]
• Modify the .h file
  #define TS_MPDATA
  and
  #undef TS_U3HADVECTION
  #undef TS_SVADVECTION
  (or remove) any other horizontal advection switches. This will be the only change, so recompile and make sure it works again after this small change...if something goes awry, it is good not to have accumulated lots of changes.

Turn on sediment

What is different about sediment?

External differences

The external differences are required changes in user input and output. They include:

• Switches in .h file. A summary of sediment-related cppdefs is listed in sediment_cppdefs.h
• Changes in .in file
• Need sediment.in file
• Need to modify (a copy of) ana_sediment.F
• New output in his.nc file

Many applications also need

• More info in init.nc file (initial bed)
• Wave input
• River sediment input
• Boundary conditions

Internal differences

Components of the sediment model include:

• TS_MPDATA – Positive-definite advection scheme
• Wave model (as input via a forcing file, or directly coupled)
• Wave-current combined bottom stresses
• Erosion / deposition / bed model (sand)
• Settling
• Bedload transport and flux divergence
• Morphological evolution
• Wetting and drying
• Cohesive and mixed sediments

Simplest possible case (Run 02)

It should be possible to turn sediment on by making the following simple changes. However, this did not work for the Latte case, and you will have take the additional steps. For now, skip to Run 03 and when this is working correctly, we'll describe it here.

In the latte_c.h file:

#define TS_MPDATA
#define sediment
#define ana_sediment

In the roms_lattec.in file:

NNS = 1  ! Number of non-cohesive (sand) sediment tracers
...
SPARNAM = ./sediment_lattec.in

The best way to get a template for the sediment input file is to copy one from the source code to your local directory.

cp ~/yoursourcedir/ROMS/External/sediment_sed_test1.in ./run02/sediment_lattec.in

then edit to suit.

Make an init file (Run03)

• Copy the .h and .in files from Run01 to a directory called ../Latte/Run03
• Edit the .h file to activate sediment and the SSW version of the bottom boundary layer calculations. First, the following declarations are needed. Some of these were already defined...the additional ones represent fluxes that are required but not provided by a sediment bed, and will be set to zero in their respective ana_xx.F files.
  #define ANA_BPFLUX /* needed when sediment and bulk fluxes are on */
  #define ANA_SPFLUX /* needed when sediment and bulk fluxes are on */
  #define ANA_BSFLUX /* bottom temp flux...defaults to zero */
  #define ANA_BTFLUX /* bottom heat flux...defaults to zero */
  #define ANA_SSFLUX /* analytical surface salt flux defaults to zero */
  Now add the declarations for the wave-current bottom stress calculations, and for sediment. We have to turn off the quadratic drag specification, because we are replacing it with SSW_BBL.
  /* stuff needed for sediment */
  #define SSW_BBL /* Sherwood et al. BBL closure */
  #ifdef SSW_BBL
  # define SSW_CALC_ZNOT /* Computing bottom roughness internally */
  # undef SSW_LOGINT /* Logarithmic interpolation of (Ur,Vr) */
  # define SSW_CALC_UB /* Computing bottom orbital velocity internally */
  # undef SSW_FORM_DRAG_COR /* Activate form drag coefficient */
  # undef SSW_ZOBIO /* Biogenic bedform roughness from ripples */
  # undef SSW_ZOBL /* Bedload roughness for ripples */
  # undef SSW_ZORIP /* Bedform roughness from ripples */
  #endif
  Finally, to reduce the output, we can also:
  #undef AVERAGES /* reduce output */
  #undef STATIONS /* reduce output */
  You can download the resulting .h file from https://coast-enviro.er.usgs.gov/Latte/run03.
• Copy a sediment.in file from an example application and make some changes.
  cp ~/src/ROMS/External/sediment_sed_test1.in ./Run03/sediment_lattec.in
  We will only have one bed layer; it should be thick enough so that it does not erode completely in most of the model domain but, if it is too thick, changes in composition during the model run will be less apparent. We are going to specify only one non-cohesive sediment class, so we need only one value for each parameter.
  NEWLAYER_THICK == 0.01d0
  
  ...
  
  !------------------------------------------------------------------------------
  ! Non-cohesive Sediment Parameters, [1:NNS,1:Ngrids] values expected.
  !------------------------------------------------------------------------------
  
  ! Median sediment grain diameter (mm).
  
  SAND_SD50 == .150d0
  
  !Sediment concentration (kg/m3).
  
  SAND_CSED == 0.0d0
  
  ! Sediment grain density (kg/m3).
  
  SAND_SRHO == 2650.0d0
  
  ! Particle settling velocity (mm/s).
  
  SAND_WSED == .20d0
  
  ! Surface erosion rate (kg/m2/s).
  
  SAND_ERATE == 5.0d-4
  
  ! Critical shear for erosion and deposition (N/m2).
  
  SAND_TAU_CE == 0.1d0
  SAND_TAU_CD == 0.1d0
  
  ! Porosity (nondimensional: 0.0-1.0): Vwater/(Vwater+Vsed).
  
  SAND_POROS == 0.5d0
  
  If we had several sediment classes, we would need more (NNS) numbers arrayed across after each variable.
• Make several changes to the .in file, starting with the title.
  TITLE = Latte Tutorial, Run03
  Specify the number of sediment tracers...in this case, one non-cohesive size.
  NCS = 0  ! Number of cohesive (mud) sediment tracers
  NNS = 1  ! Number of non-cohesive (sand) sediment tracers
  Turn on additional output to appear in the .his output file.
  ...
  Hout(idUbrs) == T  ! bottom U-current stress
  Hout(idVbrs) == T  ! bottom V-current stress
  Hout(idUbws) == T  ! bottom U-wave stress
  Hout(idVbws) == T  ! bottom V-wave stress
  Hout(idUbcs) == T  ! bottom max wave-current U-stress
  Hout(idVbcs) == T  ! bottom max wave-current V-stress
  
  Hout(idUbot) == T  ! bed wave orbital U-velocity
  Hout(idVbot) == T  ! bed wave orbital V-velocity
  Hout(idUbur) == T  ! bottom U-velocity above bed
  Hout(idVbvr) == T  ! bottom V-velocity above bed
  ...
  Hout(idBott) == T T T T T T T T T F F F F F F F
  ...
• Add more information to the initialization file (discussed in the next section), so we will point to that in the .in file
  ININAME == ../in/lattec_sed1_ini_94.nc
• We want to add some waves, so a wave forcing file (see below) is added to the list and NFFILES is increased.
  NFFILES == 10  ! number of forcing files
  
  FRCNAME == ../in/frc_latte_wrf_Lwrf.nc \
  ../in/frc_latte_wrf_Pair.nc \
  ../in/frc_latte_wrf_Qair.nc \
  ../in/frc_latte_wrf_Swrf.nc \
  ../in/frc_latte_wrf_Tair.nc \
  ../in/frc_latte_wrf_Uwind.nc \
  ../in/frc_latte_wrf_Vwind.nc \
  ../in/roms_lattec_river.nc \
  ../in/frc_tides_lattec.nc \
  ../in/frc_lattec_small_waves.nc
  
  * Also, rename the output files
   ! Output NetCDF file names, [1:Ngrids].
  
  GSTNAME == ocean_gst.nc
  RSTNAME == ../out/lattec_rst_run03.nc
  HISNAME == ../out/lattec_his_run03.nc
  

Needs waves!

In most coastal sediment-transport applications, waves are critical. Waves are often responsible for mobilizing sediment when tidal or wind-generated currents are too weak...but currents are needed to mix sediment out of the wave boundary layer and transport it. Sometimes a time series of wave input is available, and can be passed to ROMS directly as a forcing file. To make use of the wave input, one of the three available wave-current bottom boundary layer routines must be enabled. These are SSW_BBL, MB_BBL, and SG_BBL (see https://www.myroms.org/wiki/index.php/Sediment_Model#Bottom_Stress).

Wave input can be as simple as a single time series of wave height, period, and direction that applies to the whole model domain (handy for test cases and simplified models), or as complicated as calculated near-bottom orbital velocities and associated periods and directions that vary with time and space. In the first instance, the near bottom wave-orbital velocities are calculated from H_s and T using linear wave theory, using #define SSW_CALC_UB

Note that the linear wave calculations made using #define SSW_CALC_UB are seldom exactly correct (see, for example, Soulsby, 1997). A better representation of near-bottom wave-orbital velocities can be made from measured wave spectra, or from Hs and T with a parameterized spectra. Both methods are discussed, and Matlab procedure are provided by Wiberg and Sherwood (2008). These orbital velocities are e

Make and use a simple wave forcing file (Run04)

In this section, we will use buoy data to make a simple wave forcing file.

SWAN - A Digression

Getting SWAN set up correctly is easier if you do it separately. For many applications, one-way coupling (SWAN influences ROMS, but ROMS does not feed back to SWAN) captures the essential physics. So, before we couple ROMS and SWAN, we'll run SWAN indepedently. This effort won't be wasted because the coupling to ROMS is not seamless...you still have to specify the SWAN INPUT file, bathymetry, and some of the forcing for SWAN as separate files.

Directions in SWAN can be tricky. For our purposes, we will use these commands in the INPUT file:

SET NAUTICAL
COORDINATES SPHERICAL
CGRID CURVILINEAR 145 81 EXC 9.999000e+003 &
CIRCLE 36 0.0418 1.0 24
READGRID COORDINATES 1 './lattec_swan.grd' 4 0 0 FREE

NAUTICAL implies that directions are specified in degrees and (for wind and waves) indicate the direction from which the wind/waves are approaching, For example, a uniform southwesterly wind (blowing toward the northeast) at 5 m/s for a stationary run would be specified as

WIND 5 225

SWAN output in .mat files should be interpreted this way as well.

In the Latte example, we have good wind forcing files that work with ROMS, and we want to use the same files to force waves in SWAN. We will use the SWAN commands

INPGRID WIND CURVILINEAR 0 0 145 81 EXC 9.999000e+003 &
NONSTATIONARY 20060405.000000 1 HR 20060608.230000
READINP WIND 1 'swan_lattec_wind.dat' 4 0 FREE

We can translate the frc_latte_wrf_Uwind.nc and frc_latte_wrf_Vwind.nc files to SWAN input using create_swan_wind.m. The resulting file swan_lattec_wind.dat has sequences of u and v velocities written for each time step. The u- and v- components are aligned with the x- and y-axes of the SWAN grid. The order in which they are written is determined by the 'layout' parameter (in this case, 4), which means the

file looks like this.

u(1,1,1) u(1,1,2) ... u(1,1,nj)
u(1,2,1) u(1,2,2) ... u(1,2,nj)
...
u(1,ni,1) u(1,ni,2)... u(1,ni,nj)
v(1,1,1) v(1,1,2) ... v(1,1,nj)
v(1,2,1) v(1,2,2) ... v(1,2,nj)
...
v(1,ni,1) v(1,ni,2)... v(1,ni,nj)
u(1,1,1) u(1,1,2) ... u(1,1,nj)
u(1,2,1) u(1,2,2) ... u(1,2,nj)
...
u(1,ni,1) u(1,ni,2)... u(1,ni,nj)
v(2,1,1) v(2,1,2) ... v(2,1,nj)
v(2,2,1) v(2,2,2) ... v(2,2,nj)
...
v(2,ni,1) v(2,ni,2)... v(2,ni,nj)
...etc

where the first indice is time step (1...nt), the second is x-direction (1...ni) and the third is y-direction (1...nj).