Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions src/core_ocean/Registry.xml
Original file line number Diff line number Diff line change
Expand Up @@ -680,6 +680,10 @@
description="Value of period of monochromatic tide."
possible_values="Any positive real number."
/>
<nml_option name="config_tidal_forcing_monochromatic_phaseLag" type="real" default_value="0.0" units=""
description="Value of phase of monochromatic tide."
possible_values="Any real number between."
/>
<nml_option name="config_tidal_forcing_monochromatic_baseline" type="real" default_value="0.0" units="days"
description="Value of baseline monochromatic tide, e.g., sea level rise."
possible_values="Any positive real number."
Expand Down
48 changes: 48 additions & 0 deletions src/core_ocean/mode_init/Registry_tidal_boundary.xml
Original file line number Diff line number Diff line change
Expand Up @@ -59,4 +59,52 @@
description="Bottom drag of right side of the boundary."
possible_values="Any real number"
/>
<nml_option name="config_use_idealized_transect" type="logical" default_value="false" units=""
description="Logical to determine if idealized tidal flat profile is defined."
possible_values="True or False"
/>
<nml_option name="config_idealized_transect_Lshore" type="real" default_value="0.6" units=""
description="Ratio of shore length in the idealized coastal profile."
possible_values="Any positive real value between 0 and 1."
/>
<nml_option name="config_idealized_transect_Sshore" type="real" default_value="0.001" units=""
description="Shore slope."
possible_values="Any positive real value, should be a small value."
/>
<nml_option name="config_idealized_transect_Lcoast" type="real" default_value="0.3" units=""
description="Ratio of coast length in the idealized coastal profile"
possible_values="Any positive real value between 0 and 1."
/>
<nml_option name="config_idealized_transect_Scoast" type="real" default_value="0.001" units=""
description="Coast slope."
possible_values="Any positive real value."
/>
<nml_option name="config_idealized_transect_Lmarsh" type="real" default_value="0.1" units=""
description="Ratio of marsh length in the idealized coastal profile."
possible_values="Lmarsh=1-Lshore-Lcoast; any positive real value between 0 and 1."
/>
<nml_option name="config_idealized_transect_Smarsh" type="real" default_value="0.0" units=""
description="Marsh slope"
possible_values="Any real value, can be negative"
/>
<nml_option name="config_idealized_transect_roughness" type="real" default_value="0.025" units=""
description="Bottom roughness (Cd or Manning roughness) at non-vegetated region"
possible_values="Any real value, can be negative"
/>
<nml_option name="config_idealized_transect_roughness_marsh" type="real" default_value="0.075" units=""
description="Bottom roughness (Cd or Manning roughness) at vegetated region"
possible_values="Any real value, can be negative"
/>
<nml_option name="config_idealized_vegetation_diameter" type="real" default_value="0.05" units="m"
description="Constant vegetation diameter for idealized transect case"
possible_values="Any non-negative real value"
/>
<nml_option name="config_idealized_vegetation_height" type="real" default_value="0.2" units="m"
description="Constant vegetation height for idealized transect case"
possible_values="Any non-negative real value "
/>
<nml_option name="config_idealized_vegetation_density" type="real" default_value="1000" units="m^{-2}"
description="Constant vegetation density for indealized transect case"
possible_values="Any non-negative real value"
/>
</nml_record>
77 changes: 74 additions & 3 deletions src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F
Original file line number Diff line number Diff line change
Expand Up @@ -93,23 +93,25 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{
type (mpas_pool_type), pointer :: verticalMeshPool
type (mpas_pool_type), pointer :: tracersPool

integer :: iCell, k
integer :: iCell, k, N

! Define dimensions
integer, pointer :: nCellsSolve, nEdgesSolve, nVertLevels, nVertLevelsP1
integer, pointer :: index_temperature, index_salinity, index_tracer1

! Define arrays
integer, dimension(:), pointer :: maxLevelCell
integer, dimension(:), pointer :: vegetationMask
real (kind=RKIND), dimension(:), pointer :: vegetationHeight,vegetationDensity, vegetationDiameter
real (kind=RKIND), dimension(:), pointer :: yCell, refBottomDepth, bottomDepth, vertCoordMovementWeights, dcEdge
real (kind=RKIND), dimension(:), pointer :: tidalInputMask
real (kind=RKIND), dimension(:), pointer :: bottomDrag
real (kind=RKIND), dimension(:), pointer :: ssh
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness, zMid
real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers, debugTracers

real (kind=RKIND), dimension(:), pointer :: interfaceLocations
real (kind=RKIND), parameter :: eps=1.0e-12
real (kind=RKIND) :: cff1,cff2,cff3, dep_mark1, dep_mark2 ! intermediate variables

iErr = 0

Expand Down Expand Up @@ -225,6 +227,69 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{
end do
end if

if (config_use_vegetation_drag) then
call mpas_pool_get_array(forcingPool, 'vegetationMask', vegetationMask)
vegetationMask = 0
call mpas_pool_get_array(forcingPool, 'vegetationDiameter', vegetationDiameter)
vegetationDiameter = config_idealized_vegetation_diameter
call mpas_pool_get_array(forcingPool, 'vegetationHeight', vegetationHeight)
vegetationHeight = config_idealized_vegetation_height
call mpas_pool_get_array(forcingPool, 'vegetationDensity', vegetationDensity)
vegetationDensity = config_idealized_vegetation_density
endif

if (config_use_idealized_transect) then
call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag)
config_idealized_transect_Lmarsh = 1.0_RKIND - config_idealized_transect_Lcoast &
- config_idealized_transect_Lshore
if (config_idealized_transect_Lmarsh .lt. 0.0_RKIND) then
call mpas_log_write("Lshore+Lcoast cannot be bigger than 1.0")
iErr = 1
endif
cff1 = yMax*config_idealized_transect_Lcoast
cff2 = yMax*config_idealized_transect_Lmarsh
cff3 = yMax*config_idealized_transect_Lshore
! by defining the slopes and left_bottom_depth, the pre-defined right_bottom_depth won't work.
! Need redefine it.
config_tidal_boundary_right_bottom_depth = config_tidal_boundary_left_bottom_depth + &
cff1*config_idealized_transect_Scoast + &
cff2*config_idealized_transect_Smarsh + &
cff3*config_idealized_transect_Sshore
do iCell = 1, nCellsSolve
if (yCell(iCell) .lt. cff1) then
bottomDepth(iCell) = config_tidal_boundary_left_bottom_depth &
+ (yCell(iCell)-yMin)*config_idealized_transect_Scoast
dep_mark1 = bottomDepth(iCell)
bottomDrag(iCell) = config_idealized_transect_roughness
elseif (yCell(iCell) .lt. (cff1+cff2)) then
bottomDepth(iCell) = dep_mark1 + (yCell(iCell)-cff1)*config_idealized_transect_Smarsh
dep_mark2 = bottomDepth(iCell)
bottomDrag(iCell) = config_idealized_transect_roughness_marsh
if (config_use_vegetation_drag) vegetationMask(iCell) = 1
else
if (cff2 .eq. 0.0_RKIND) dep_mark2 = dep_mark1
bottomDepth(iCell) = dep_mark2 + (yCell(iCell)-cff1-cff2)*config_idealized_transect_Sshore
bottomDrag(iCell) = config_idealized_transect_roughness
endif
end do !! do iCell
else
! if not config_idealized_transect, assign a constant slope
do iCell = 1, nCellsSolve
bottomDepth(iCell) = config_tidal_boundary_left_bottom_depth &
+ (yCell(iCell) - yMin) / (yMax - yMin) * &
(config_tidal_boundary_right_bottom_depth - config_tidal_boundary_left_bottom_depth)
end do
end if !! if config_idealized_transect

if (config_tidal_boundary_right_bottom_depth < config_tidal_boundary_left_bottom_depth) then
call mpas_log_write('Right boundary must be deeper than left boundary!', MPAS_LOG_CRIT)
end if

! Set refBottomDepth, bottomDepth, and maxLevelCell
do k = 1, nVertLevels
refBottomDepth(k) = config_tidal_boundary_right_bottom_depth * interfaceLocations(k+1)
end do

if (config_use_wetting_drying .and. config_tidal_start_dry .and. &
trim(config_tidal_boundary_layer_type) == 'zstar') then
do iCell = 1, nCellsSolve
Expand Down Expand Up @@ -301,6 +366,7 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{
if (config_use_wetting_drying .and. config_tidal_start_dry) then
do iCell = 1, nCellsSolve
ssh(iCell) = -bottomDepth(iCell) + config_drying_min_cell_height*maxLevelCell(iCell)
ssh(iCell) = MAX(ssh(iCell),-config_tidal_forcing_monochromatic_baseline)
! also computes zMid
do k = 1, maxLevelCell(iCell)
layerThickness(k,iCell) = (ssh(iCell) + bottomDepth(iCell))/maxLevelCell(iCell)
Expand All @@ -327,14 +393,19 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{
! Set tidal boundary mask
do iCell = 1, nCellsSolve
tidalInputMask(iCell) = 0.0_RKIND
if (yCell(iCell) > (25.0e3 - dcEdgeMinGlobal/2.0_RKIND)) then
if (yCell(iCell) > (yMax - dcEdgeMinGlobal/2.0_RKIND)) then
tidalInputMask(iCell) = 1.0_RKIND
! spread it over multiple cells
!if (yCell(iCell) > (25.0e3 - 3*dcEdgeMinGlobal)) then
! tidalInputMask(iCell) = exp(-((yCell(iCell)-25.0e3)/dcEdgeMinGlobal)**2.0)
end if
end do

! check that there is some tidalInputMask
if (.not. sum(tidalInputMask) > 0) then
call mpas_log_write('Input mask for tidal case is not set!', MPAS_LOG_CRIT)
end if

! Set salinity
if ( associated(activeTracers) ) then
do iCell = 1, nCellsSolve
Expand Down
3 changes: 2 additions & 1 deletion src/core_ocean/shared/mpas_ocn_tidal_forcing.F
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,8 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, diagnost
! compute the tidalHeight
if (trim(config_tidal_forcing_model) == 'monochromatic') then
tidalHeight = config_tidal_forcing_monochromatic_amp * &
SIN(2.0_RKIND*pii/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim) - &
SIN(2.0_RKIND*pii/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim - &
pii*config_tidal_forcing_monochromatic_phaseLag/180.0_RKIND) - &
config_tidal_forcing_monochromatic_baseline
!else if (trim(config_tidal_forcing_type) == 'data') then
! ! data option
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
<template>
<namelist>
<option name="config_Rayleigh_friction">.false.</option>
<option name="config_Rayleigh_damping_depth_variable">.false.</option>
<option name="config_use_implicit_bottom_drag">.false.</option>
<option name="config_use_implicit_bottom_drag_variable">.false.</option>
<option name="config_use_implicit_bottom_drag_variable_mannings">.true.</option>
<option name="config_thickness_flux_type">'centered'</option>
<option name="config_dt">'0000_00:00:05'</option>
<option name="config_run_duration">'0000_48:00:01'</option>
<option name="config_tidal_forcing_monochromatic_amp">1.5</option>
<option name="config_tidal_forcing_monochromatic_period">1.0</option>
<option name="config_tidal_forcing_monochromatic_phaseLag">0.0</option>
<option name="config_tidal_forcing_monochromatic_baseline">3.0</option>
</namelist>
</template>
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
<template>
<streams>
<stream name="forcing">
<attribute name="type">input</attribute>
<attribute name="input_interval">initial_only</attribute>
<attribute name="filename_template">forcing.nc</attribute>
<add_contents>
<member type="var" name="vegetationMask"/>
<member type="var" name="vegetationDiameter"/>
<member type="var" name="vegetationHeight"/>
<member type="var" name="vegetationDensity"/>
</add_contents>
</stream>
<stream name="output">
<attribute name="type">output</attribute>
<attribute name="clobber_mode">truncate</attribute>
<attribute name="filename_template">output.nc</attribute>
<attribute name="output_interval">0000-00-00_00:12:00</attribute>
<add_contents>
<member type="var" name="vegetationManning"/>
</add_contents>
</stream>
</streams>
</template>
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
<template>
<namelist>
<option name="config_use_variable_drag">.true.</option>
<option name="config_use_idealized_transect">.true.</option>
<option name="config_idealized_transect_Lcoast">0.3</option>
<option name="config_idealized_transect_Scoast">0.00125</option>
<option name="config_idealized_transect_Lshore">0.2</option>
<option name="config_idealized_transect_Sshore">0.0025</option>
<option name="config_idealized_transect_Lmarsh">0.5</option>
<option name="config_idealized_transect_Smarsh">0.0</option>
<option name="config_idealized_transect_roughness">0.025</option>
<option name="config_idealized_transect_roughness_marsh">0.025</option>
<option name="config_tidal_boundary_left_bottom_depth">0.0</option>
<option name="config_tidal_boundary_left_value">0.0</option>
<option name="config_tidal_boundary_right_value">8.0e3</option>
<option name="config_tidal_forcing_monochromatic_baseline">3.0</option>
<option name="config_alter_ICs_for_pbcs">.true.</option>
<option name="config_thickness_flux_type">'centered'</option>
</namelist>
</template>
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
<template>
<streams>
<stream name="forcing_data_init">
<attribute name="type">output</attribute>
<attribute name="filename_template">forcing.nc</attribute>
<attribute name="clobber_mode">truncate</attribute>
<attribute name="output_interval">0000_00:00:01</attribute>
<add_contents>
<member name="vegetationMask" type="var"/>
<member name="vegetationDiameter" type="var"/>
<member name="vegetationHeight" type="var"/>
<member name="vegetationDensity" type="var"/>
</add_contents>
</stream>
</streams>
</template>
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
#!/usr/bin/env python
# coding: utf-8
# Authors: Zhendong Cao, Phillip Wolfram
# Date: May 2020

# import libs
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import os
#import cv2
# plot tidal amplitude at open boundary

fname1 = 'NoVegetation'
fname2 = 'ConstantVegManning'
fname3 = 'VegManningEquation'

# mpas-o
ds1 = xr.open_dataset(fname1+'.nc')
ds2 = xr.open_dataset(fname2+'.nc')
ds3 = xr.open_dataset(fname3+'.nc')

x = np.linspace(0, 1.0, len(ds1.xtime))*48.0
ympas1 = ds1.ssh.where(ds1.tidalInputMask).mean('nCells').values
ympas2 = ds2.ssh.where(ds2.tidalInputMask).mean('nCells').values
ympas3 = ds3.ssh.where(ds3.tidalInputMask).mean('nCells').values

plt.plot(x, ympas1, marker='o', label=fname1)
plt.plot(x, ympas2, marker='^', label=fname2)
plt.plot(x, ympas3, marker='*', label=fname3)

# analytical
x = np.linspace(0,48.0,100)
y = 1.5*np.sin(x*np.pi/12.0) - 3.0
plt.plot(x, y, 'k-', lw=3, label='analytical')

plt.legend(frameon=False, loc='upper right', fontsize=8)
plt.ylabel('Tidal amplitude (m)')
plt.xlabel('Time (hrs)')
plt.suptitle('Tidal amplitude forcing for test cases and analytical')
plt.savefig('tidalcomparison.png',dpi=300)
plt.close()

# extract bottomDepth
bottomDepth = ds1.bottomDepth.values.reshape((116, 4))[:,1]
x_along = np.linspace(0,8,116)
# plot water elevation at selected time slices
times = np.linspace(0.0,2.0,41)
ii = 0
figname_prefix = 'Waterlevel'
for atime in times:
plt.figure(figsize=(12,6))
plt.xlim(0, 8)
plt.ylim(-4, 0)
plottime = int((float(atime)/0.2)*24.0)
timestr = ds1.isel(Time=plottime).xtime.values
ymean = ds1.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS)
x = ymean.yCell.values/1000.0
y = ymean.ssh.values
plt.plot(x, y,'-k',label=fname1)

ymean = ds2.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS)
x = ymean.yCell.values/1000.0
y = ymean.ssh.values
plt.plot(x, y,'-b',label=fname2)

ymean = ds3.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS)
x = ymean.yCell.values/1000.0
y = ymean.ssh.values
plt.plot(x, y,'-g',label=fname3)

plt.fill_between(x=x_along,y1=-10*np.ones(x_along.shape), y2=-bottomDepth, color='grey',zorder=9)
plt.xlabel('Cross shore distance (km)')
plt.ylabel('Bathymetry (m)')
atime='%.2f'%atime
plt.title('time='+atime+'days')
plt.legend(frameon=False, loc='upper right')
figname = '{}{:02d}.png'.format(figname_prefix,ii)
plt.savefig(figname, dpi=100)
plt.close()
print('Time '+atime+ ' done!')
ii += 1
# make a movie
os.system("ffmpeg -r 1 -i Waterlevel%02d.png -vcodec mpeg4 -y movie.mp4")
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
<?xml version="1.0"?>
<config case="analysis">
<add_link path_base="script_test_dir" source="comparison.py" dest="comparison.py"/>
<add_link source="../forward1/output.nc" dest="NoVegetation.nc"/>
<add_link source="../forward2/output.nc" dest="ConstantVegManning.nc"/>
<add_link source="../forward3/output.nc" dest="VegManningEquation.nc"/>

<run_script name="run.py">
<step executable="./comparison.py">
</step>
<step executable="paraview_vtk_field_extractor.py">
<argument flag="-f">output1.nc</argument>
<argument flag="-o">vtk_output1</argument>
<argument flag="-v">allOnCells</argument>
<argument flag="-d">maxEdges=0</argument>
<argument flag="">nVertLevels=0:10</argument>
<argument flag="--combine"></argument>
</step>
<step executable="paraview_vtk_field_extractor.py">
<argument flag="-f">output2.nc</argument>
<argument flag="-o">vtk_output2</argument>
<argument flag="-v">allOnCells</argument>
<argument flag="-d">maxEdges=0</argument>
<argument flag="">nVertLevels=0:10</argument>
<argument flag="--combine"></argument>
</step>
</run_script>
</config>
Loading