Difference between revisions of "FJORD TIDAL CASE"

From WikiROMS
Jump to navigationJump to search
Line 56: Line 56:


===Set up USER SETTINGS===
===Set up USER SETTINGS===
Before you run EASYGRID to create your ROMS grid, you need to edit EASYGRID's user settings, which is the first part of script (after the help lines). is to edit t
Before you run EASYGRID to create your ROMS grid, you need to edit EASYGRID's user settings, which is the first part of script (after the help lines). Beside each entry line there are detailed instructions. Below is a copy of EASYGRID's user settings / instructions. The default settings are for this Fjord test case, hence you can just run EASYGRID (after downloading the bathymetry and coastline files above) to generate the grid and initial-conditions files required for this tutorial.


  <span style="color:green; font-size:85%; line-height:90%;">%*******************************************************************************************************************************
  <span style="color:green; font-size:85%; line-height:90%;">%*******************************************************************************************************************************
Line 74: Line 74:
  %
  %
  % Geographical and Grid parameters --------
  % Geographical and Grid parameters --------
       <span class="black">lat =  44.8125;</span>         % Latitude  (degrees) of the bottom-left corner of the grid.
       <span class="black">lat =  44.8125;</span>       % Latitude  (degrees) of the bottom-left corner of the grid.
       <span class="black">lon = -62.8855;</span>         % Longitude (degrees) of the bottom-left corner of the grid.  
       <span class="black">lon = -62.8855;</span>       % Longitude (degrees) of the bottom-left corner of the grid.  
  %
  %
         <span class="black">X = 9100;</span>            % Width of domain (meters)
         <span class="black">X = 9100;</span>            % Width of domain (meters)
Line 110: Line 110:
  % OUTPUT: File naming and NetCDF descriptors ------------------------------
  % OUTPUT: File naming and NetCDF descriptors ------------------------------
  % -------------------------------------------------------------------------
  % -------------------------------------------------------------------------
  <span class="black">Grid_filename =</span> <span class="purple">'Fjord'</span>; % Filename for grid and initial conditions files (don't include extension).  
  <span class="black">Grid_filename =</span> <span class="purple">'Fjord'</span>;     % Filename for grid and initial conditions files (don't include extension).  
                            % "_grd.nc" is added to grid file and "_ini.nc" is added to initial conditions file
                              % "_grd.nc" is added to grid file and "_ini.nc" is added to initial conditions file
  <span class="black">Descrip_grd  =</span> <span class="purple">'Test grid'</span>;              %Description for grid .nc file
  <span class="black">Descrip_grd  =</span> <span class="purple">'Test grid'</span>;              %Description for grid .nc file
  <span class="black">Descrip_ini  =</span> <span class="purple">'Test initial conditions'</span>; %Description for initial conditions .nc file
  <span class="black">Descrip_ini  =</span> <span class="purple">'Test initial conditions'</span>; %Description for initial conditions .nc file

Revision as of 16:36, 16 April 2008

Fjord Tidal Test Case
Warning This WikiROMS article is currently under HEAVY construction. This message will be erased when active construction is finished (in 1 or 2 days).
April 14, 2008

This tutorial will go over some of the basic steps to set up a ROMS realistic application (yet, this is a very simple one). This tutorial is also a demonstration of how to use EASYGRID, a "quick-and-dirty" matlab script to make ROMS grids and initialization files.


Warning PREREQUISITES: This tutorial assumes that you have (1) downloaded ROMS, (2) installed it in your computer (or cluster), (3) tested it by compiling and running one of the included test cases, and (4) installed MEXNC, SNCTOOLS and the ROMS Matlab tool-kit. If you haven't done all the above, check the "Getting Started" and "Tutorials" WikiROMS sections.



Geographical Preamble: This application is for Ship Harbour, an estuarine fjord in Nova Scotia, Canada. Click here to see the location in Google. Tides are semidiurnal and tidal range is 1.4 m on average and 2 m on spring tides. For now I will only include tidal forcing, however, there is a river at the uppermost end of the estuary, which discharges freshwater at an annual average rate of 18 m3 s-1. I plan to write another tutorial on how to add a river, but for now is only tides.



Grid Generation

The first step to set up a realistic application is to set up a realistic grid. There are several software packages to generate ROMS grids. Here I will use EASYGRID.


Download EASYGRID

  • Download EASYGRID here (I will add link a bit later...)
  • Add the location where you placed EASYGRID to Matlab's search path.
  • Dependencies: EASYGRID uses MEXNC, SNCTOOLS and a few scripts (like spheric_dist.m and smth_bath.m) from the ROMS Matlab tool-kit


Get bathymetry

You will need a bathymetry .mat file containing 3 vectors (xbathy, ybathy and zbathy), where xbathy = Longitude, ybathy = Latitude and zbathy = depth (positive, in meters). Vectors xbathy and ybathy are in decimal degrees. For more information about bathymetry for ROMS, click here.


There easiest way to get bathymetry data (that I know) is to use Rich Signell's read_srtm30plus.m Matlab script. Instant bathymetry in 1 line of Matlab code! Click HERE for instruction to get and use read_srtm30plus.m.


The most difficult way to get bathymetry data (that I know) is to scan a chart and then digitize it using your mouse and a digitizing software like Didger (sorry, I am not aware of any good digitizing free software).


For this tutorial, you should download this bathymetry file (I will add link soon). I created it by merging two bathymetry files... one high-resolution digitized from a chart and another coarse-resolution using read_srtm30plus.m


Get coastline

In EASYGRID, the coastline is only used for plotting and therefore not indispensable. However, a later step of mask editing (using editmask.m) requires a coastline file. For both, EASYGRID and editmask.m, the coastline should be in a .mat file and should contain 2 vectors named "lat" and "lon" (both in units of decimal degrees).


An easy way to get coastline data is at NOAA's Coastline Extractor]. Remember to select Matlab as format option (not default). Again, SEAGRID's wikiROMS page has excellent tools and descriptions on how to use the Coastline Extractor...


For this tutorial, you should download this coastline file (I will add link soon). I created from the digitized chart I use to meka the bathymetry file.


Set up USER SETTINGS

Before you run EASYGRID to create your ROMS grid, you need to edit EASYGRID's user settings, which is the first part of script (after the help lines). Beside each entry line there are detailed instructions. Below is a copy of EASYGRID's user settings / instructions. The default settings are for this Fjord test case, hence you can just run EASYGRID (after downloading the bathymetry and coastline files above) to generate the grid and initial-conditions files required for this tutorial. 

%*******************************************************************************************************************************
% USER SETTINGS ****************************************************************************************************************
%*******************************************************************************************************************************
%
% Switches -------------------------
     save_grid = 1;   % Save GRID in a NetCDF file ( yes = 1; no = 0 )
     save_init = 1;   % Create (and save) INITIAL CONDITIONS (from grid) ( yes = 1; no = 0 )
 screen_output = 1;   % Displays (on the screen) many parameters that need to be copy-pasted in the ocean.in file ( yes = 1; no = 0 )
     plot_grid = 1;   % Plot grid ( yes = 1; no = 0 )
   smooth_grid = 1;   % Smooth bathymetry ( yes = 1; no = 0 ) using H. Arango's smth_bath.m , which applies a Shapiro filter to the bathymetry
%
% -------------------------------------------------------------------------
% GRID SETTINGS -----------------------------------------------------------
% -------------------------------------------------------------------------
%
% Geographical and Grid parameters --------
     lat =  44.8125;        % Latitude  (degrees) of the bottom-left corner of the grid.
     lon = -62.8855;        % Longitude (degrees) of the bottom-left corner of the grid. 
%
       X = 9100;            % Width of domain (meters)
       Y = 2550;            % Length of domain (meters)
rotangle = -43;             % Angle (degrees) to rotate the grid conterclock-wise
   resol = 30;              % Cell width and height (i.e. Resolution)in meters. Grid cells are forced to be (almost) square.
       N = 10;              % Number of vertical levels
%   
%     
% Bathymetry -------------- % Bathymetry for ROMS is measured positive downwards (zeros are not allowed) see: https://www.myroms.org/wiki/index.php/Grid_Generation#Bathymetry
                            % To allow variations in surface elevation (eg. tides) while keeping all depths positive,
                            % an arbitrary offset (see minh below) is added to the depth vector.
%    
      hh = nan;             % Analytical Depth (meters) used to create a uniform-depth grid. If using a bathymetry file, leave hh = nan;
    minh = 4;               % Arbitrary depth offset in meters (see above). minh should be a little more than the maximum expected tidal variation.
   Bathy = 'Test_bathy.mat';% Bathymetry file. If using the analytical depth above (i.e. hh ~= nan), then Bathy will not be used.
                            % The bathymetry file should be a .mat file containing 3 vectors (xbathy, ybathy and zbathy). where xbathy = Longitude, 
                            % ybathy = Latitude and zbathy = depth (positive, in meters). xbathy and ybathy are in decimal degrees See: http://en.wikipedia.org/wiki/Decimal_degrees
       %Bathymetry smoothing routine...  See "Switches section" (above) to turn this ON or OFF
       if smooth_grid == 1;
           order = 2;       % Order of Shapiro filter (2,4,8)... default: 2
            rlim = 0.35;    % Maximum r-factor allowed (0.35)... default: 0.35
           npass = 50;      % Maximum number of passes.......... default: 50
       end
%---------------------------------
%                      
%
% Coastline ----------------------
   Coast = load('Test_coast.mat'); % If there isn't a coastline file... comment-out this line: e.g. %Coast = load('test_coast.mat');
                                   % The coastline is only used for plotting. The coastline .mat file should contain 2 vectors named "lat" and "lon"
%                                    
%      
% -------------------------------------------------------------------------
% OUTPUT: File naming and NetCDF descriptors ------------------------------
% -------------------------------------------------------------------------
Grid_filename = 'Fjord';      % Filename for grid and initial conditions files (don't include extension). 
                              % "_grd.nc" is added to grid file and "_ini.nc" is added to initial conditions file
Descrip_grd   = 'Test grid';               %Description for grid .nc file
Descrip_ini   = 'Test initial conditions'; %Description for initial conditions .nc file
Author        = 'John Smith';
Computer      = 'My Computer';
%
% -------------------------------------------------------------------------
% INITIAL CONDITIONS ------------------------------------------------------
% -------------------------------------------------------------------------
if save_init == 1; % See "Switches section" (above) to turn this ON or OFF
    % Initial conditions will be constant throught the grid domain
    %--------------------------------------------------------------------------
    NH4          = 0.1;     % Ammonium concentration (millimole_NH4 meter-3)
    NO3          = 10;      % Nitrate concentration (millimole_N03 meter-3)
    chlorophyll1 = 0.3;     % Chlorophyll concentration in small phytoplankyon (milligrams_chlorophyll meter-3)
    chlorophyll2 = 0.3;     % Chlorophyll concentration in large phytoplankyon (milligrams_chlorophyll meter-3)
    detritus1    = 0.03;    % Small detritus concentration (millimole_nitrogen meter-3)
    detritus2    = 0.03;    % Large detritus concentration (millimole_nitrogen meter-3)
    detritusC1   = 1;       % Small detritus carbon concentration (millimole_carbon meter-3)
    detritusC2   = 0.2;     % Large detritus carbon concentration (millimole_carbon meter-3)
    phyto1       = 0.02;    % Small phytoplankton concentration (millimole_nitrogen meter-3)
    phyto2       = 0.02;    % Large phytoplankton concentration (millimole_nitrogen meter-3)
    phytoC1      = 0.2;     % Small phytoplankton carbon concentration (millimole_carbon meter-3)
    phytoC2      = 0.1;     % Small phytoplankton carbon concentration (millimole_carbon meter-3)
    salt         = 30;      % Salinity (PSU)
    temp         = 9;       % Potential temperature (Celsius)
    u            = 0;       % u-momentum component (meter second-1)
    ubar         = 0;       % Vertically integrated u-momentum component (meter second-1)
    v            = 0;       % v-momentum component (meter second-1)
    vbar         = 0;       % Vertically integrated v-momentum component (meter second-1)
    zeta         = 0;       % Free-surface (meters)
    zooplankton  = 0.01;    % Zooplankton concentration "millimole_nitrogen meter-3"
    zooplanktonC = 0.5;     % Zooplankton carbon concentration "millimole_carbon meter-3"
    %--------------------------------------------------------------------------
end
%
%
% 
%*******************************************************************************************************************************
% END OF USER SETTINGS *********************************************************************************************************
%*******************************************************************************************************************************

Initialization

Dummy text... No real information here yet

Compiling ROMS

Before we compile ROMS, we need to create a header (.h) file and to modify some analytical Fortran files.

Create header (.h) file

Create a copy of the basin.h test-case header file, and rename the copy fjord.h


fjord.h file: Use basin.h as a template. Erase red and Add green (black remains the same). 

/*
** svn $Id: basin.h 139 2008-01-10 00:17:29Z arango $
*******************************************************************************
** Copyright (c) 2002-2008 The ROMS/TOMS Group **
** Licensed under a MIT/X style license **
** See License_ROMS.txt **
*******************************************************************************
**
** Options for Big Bad Basin.
** Options for Tidal Fjord.
**
** Application flag: BASIN
** Application flag: FJORD
** Input script: ocean_basin.in
** Input script: ocean_fjord.in
*/

#define UV_ADV
#define UV_COR
#define UV_QDRAG
#define UV_VIS4
#define MIX_S_UV
#define DJ_GRADPS
#define TS_U3HADVECTION
#define TS_C4VADVECTION
#define SOLVE3D
#define SPLINES
#define EASTERN_WALL
#define WESTERN_WALL
#define SOUTHERN_WALL
#define NORTHERN_WALL
#define BODYFORCE
#define ANA_GRID
#define ANA_INITIAL
#define ANA_SMFLUX
#define ANA_STFLUX
#define MASKING
#define EAST_FSCHAPMAN
#define EAST_M2FLATHER
#define EAST_M3RADIATION
#define EAST_TRADIATION
#define ANA_FSOBC
#define ANA_M2OBC

Modify ana_fsobc.h

fjord.h file: Use basin.h as a template. Erase red and Add green (black remains the same). 

#elif defined TEST_CHAN
     IF (WESTERN_EDGE) THEN
       cff=0.0_r8
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_west(j)=cff
       END DO
     END IF
     IF (EASTERN_EDGE) THEN
       cff=-0.4040_r8*MIN(time(ng)/150000.0_r8,1.0_r8)
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_east(j)=cff
       END DO
     END IF
#elif defined WEDDELL
     IF (WESTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       val=0.53_r8+(0.53_r8-0.48_r8)/REAL(Iend+1,r8)
       phase=(277.0_r8+(277.0_r8-240.0_r8)/REAL(Iend+1,r8))*deg2rad
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_west(j)=fac*val*COS(omega-phase)
       END DO
     END IF
     IF (EASTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       val=0.53_r8+(0.53_r8-0.48_r8)
       phase=(277.0_r8+(277.0_r8-240.0_r8))*deg2rad
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_east(j)=fac*val*COS(omega-phase)
       END DO
     END IF
#elif defined FJORD
     IF (WESTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       val=0.53_r8+(0.53_r8-0.48_r8)/REAL(Iend+1,r8)
       phase=(277.0_r8+(277.0_r8-240.0_r8)/REAL(Iend+1,r8))*deg2rad
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_west(j)=fac*val*COS(omega-phase)
       END DO
     END IF
     IF (EASTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       val=0.53_r8+(0.53_r8-0.48_r8)
       phase=(277.0_r8+(277.0_r8-240.0_r8))*deg2rad
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_east(j)=fac*val*COS(omega-phase)
       END DO
     END IF      
#else
     IF (EASTERN_EDGE) THEN
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_east(j)=0.0_r8
       END DO
     END IF
     IF (WESTERN_EDGE) THEN
       DO j=JstrR,JendR
         BOUNDARY(ng)%zeta_west(j)=0.0_r8
       END DO
     END IF
     IF (SOUTHERN_EDGE) THEN
       DO i=IstrR,IendR
         BOUNDARY(ng)%zeta_south(i)=0.0_r8
       END DO
     END IF
     IF (NORTHERN_EDGE) THEN
       DO i=IstrR,IendR
         BOUNDARY(ng)%zeta_north(i)=0.0_r8
       END DO
     END IF
#endif
     RETURN
     END SUBROUTINE ana_fsobc_tile

Modify ana_m2obc.h

Text 

#elif defined WEDDELL
     IF (WESTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       minor=0.0143_r8+(0.0143_r8+0.010_r8)/REAL(Iend+1,r8)
       major=0.1144_r8+(0.1144_r8-0.013_r8)/REAL(Iend+1,r8)
       phase=(318.0_r8+(318.0_r8-355.0_r8)/REAL(Iend+1,r8))*deg2rad
       angle=(125.0_r8+(125.0_r8- 25.0_r8)/REAL(Iend+1,r8))*deg2rad
       DO j=JstrR,JendR
         val=0.5_r8*(angler(Istr-1,j)+angler(Istr,j))
         BOUNDARY(ng)%ubar_west(j)=fac*(major*COS(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         SIN(omega-phase))
       END DO
       DO j=Jstr,JendR
         val=0.5_r8*(angler(Istr-1,j-1)+angler(Istr-1,j))
         BOUNDARY(ng)%vbar_west(j)=fac*(major*SIN(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         COS(omega-phase))
       END DO
     END IF
     IF (EASTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       minor=0.0143_r8+(0.0143_r8+0.010_r8)
       major=0.1144_r8+(0.1144_r8-0.013_r8)
       phase=(318.0_r8+(318.0_r8-355.0_r8))*deg2rad
       angle=(125.0_r8+(125.0_r8- 25.0_r8))*deg2rad
       DO j=JstrR,JendR
         val=0.5_r8*(angler(Iend,j)+angler(Iend+1,j))
         BOUNDARY(ng)%ubar_east(j)=fac*(major*COS(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         SIN(omega-phase))
       END DO
       DO j=Jstr,JendR
         val=0.5_r8*(angler(Iend+1,j-1)+angler(Iend+1,j))
         BOUNDARY(ng)%vbar_east(j)=fac*(major*SIN(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         COS(omega-phase))
       END DO
     END IF
#elif defined FJORD
     IF (WESTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       minor=0.0143_r8+(0.0143_r8+0.010_r8)/REAL(Iend+1,r8)
       major=0.1144_r8+(0.1144_r8-0.013_r8)/REAL(Iend+1,r8)
       phase=(318.0_r8+(318.0_r8-355.0_r8)/REAL(Iend+1,r8))*deg2rad
       angle=(125.0_r8+(125.0_r8- 25.0_r8)/REAL(Iend+1,r8))*deg2rad
       DO j=JstrR,JendR
         val=0.5_r8*(angler(Istr-1,j)+angler(Istr,j))
         BOUNDARY(ng)%ubar_west(j)=fac*(major*COS(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         SIN(omega-phase))
       END DO
       DO j=Jstr,JendR
         val=0.5_r8*(angler(Istr-1,j-1)+angler(Istr-1,j))
         BOUNDARY(ng)%vbar_west(j)=fac*(major*SIN(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         COS(omega-phase))
       END DO
     END IF
     IF (EASTERN_EDGE) THEN
       fac=TANH((tdays(ng)-dstart)/1.0_r8)
       omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8)  !  M2 Tide period
       minor=0.0143_r8+(0.0143_r8+0.010_r8)
       major=0.1144_r8+(0.1144_r8-0.013_r8)
       phase=(318.0_r8+(318.0_r8-355.0_r8))*deg2rad
       angle=(125.0_r8+(125.0_r8- 25.0_r8))*deg2rad
       DO j=JstrR,JendR
         val=0.5_r8*(angler(Iend,j)+angler(Iend+1,j))
         BOUNDARY(ng)%ubar_east(j)=fac*(major*COS(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         SIN(omega-phase))
       END DO
       DO j=Jstr,JendR
         val=0.5_r8*(angler(Iend+1,j-1)+angler(Iend+1,j))
         BOUNDARY(ng)%vbar_east(j)=fac*(major*SIN(angle-val)*          &
    &                                         COS(omega-phase)-        &
    &                                   minor*SIN(angle-val)*          &
    &                                         COS(omega-phase))
       END DO
     END IF      
#else
     IF (EASTERN_EDGE) THEN
       DO j=JstrR,JendR
         BOUNDARY(ng)%ubar_east(j)=0.0_r8
       END DO
       DO j=Jstr,JendR
         BOUNDARY(ng)%vbar_east(j)=0.0_r8
       END DO
     END IF
     IF (WESTERN_EDGE) THEN
       DO j=JstrR,JendR
         BOUNDARY(ng)%ubar_west(j)=0.0_r8
       END DO
       DO j=Jstr,JendR
         BOUNDARY(ng)%vbar_west(j)=0.0_r8
       END DO
     END IF
     IF (SOUTHERN_EDGE) THEN
       DO i=Istr,IendR
         BOUNDARY(ng)%ubar_south(i)=0.0_r8
       END DO
       DO i=IstrR,IendR
         BOUNDARY(ng)%vbar_south(i)=0.0_r8
       END DO
     END IF
     IF (NORTHERN_EDGE) THEN
       DO i=Istr,IendR
         BOUNDARY(ng)%ubar_north(i)=0.0_r8
       END DO
       DO i=IstrR,IendR
         BOUNDARY(ng)%vbar_north(i)=0.0_r8
       END DO
     END IF
#endif
     RETURN
     END SUBROUTINE ana_m2obc_tile

Tidal Forcing

Dummy text... No real information here yet

Retrieved from "https://www.myroms.org/wiki/index.php?title=FJORD_TIDAL_CASE&oldid=2936"