diff --git a/.doxygen b/.doxygen
index 43bb7d8eb6..fceb7fea9b 100644
--- a/.doxygen
+++ b/.doxygen
@@ -1,4 +1,4 @@
-# Doxyfile 1.8.9.1
+# Doxyfile 1.8.12
# This file describes the settings to be used by the documentation system
# doxygen (www.doxygen.org) for a project.
@@ -293,6 +293,15 @@ EXTENSION_MAPPING =
MARKDOWN_SUPPORT = YES
+# When the TOC_INCLUDE_HEADINGS tag is set to a non-zero value, all headings up
+# to that level are automatically included in the table of contents, even if
+# they do not have an id attribute.
+# Note: This feature currently applies only to Markdown headings.
+# Minimum value: 0, maximum value: 99, default value: 0.
+# This tag requires that the tag MARKDOWN_SUPPORT is set to YES.
+
+TOC_INCLUDE_HEADINGS = 0
+
# When enabled doxygen tries to link words that correspond to documented
# classes, or namespaces to their corresponding documentation. Such a link can
# be prevented in individual cases by putting a % sign in front of the word or
@@ -341,7 +350,14 @@ IDL_PROPERTY_SUPPORT = YES
# all members of a group must be documented explicitly.
# The default value is: NO.
-DISTRIBUTE_GROUP_DOC = NO
+DISTRIBUTE_GROUP_DOC = YES
+
+# If one adds a struct or class to a group and this option is enabled, then also
+# any nested class or struct is added to the same group. By default this option
+# is disabled and one has to add nested compounds explicitly via \ingroup.
+# The default value is: NO.
+
+GROUP_NESTED_COMPOUNDS = NO
# Set the SUBGROUPING tag to YES to allow class member groups of the same type
# (for instance a group of public functions) to be put as a subgroup of that
@@ -732,6 +748,12 @@ WARN_IF_DOC_ERROR = YES
WARN_NO_PARAMDOC = NO
+# If the WARN_AS_ERROR tag is set to YES then doxygen will immediately stop when
+# a warning is encountered.
+# The default value is: NO.
+
+WARN_AS_ERROR = NO
+
# The WARN_FORMAT tag determines the format of the warning messages that doxygen
# can produce. The string should contain the $file, $line, and $text tags, which
# will be replaced by the file and line number from which the warning originated
@@ -771,14 +793,38 @@ INPUT_ENCODING = UTF-8
# If the value of the INPUT tag contains directories, you can use the
# FILE_PATTERNS tag to specify one or more wildcard patterns (like *.cpp and
-# *.h) to filter out the source-files in the directories. If left blank the
-# following patterns are tested:*.c, *.cc, *.cxx, *.cpp, *.c++, *.java, *.ii,
-# *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h, *.hh, *.hxx, *.hpp,
-# *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc, *.m, *.markdown,
-# *.md, *.mm, *.dox, *.py, *.f90, *.f, *.for, *.tcl, *.vhd, *.vhdl, *.ucf,
-# *.qsf, *.as and *.js.
-
-FILE_PATTERNS =
+# *.h) to filter out the source-files in the directories.
+#
+# Note that for custom extensions or not directly supported extensions you also
+# need to set EXTENSION_MAPPING for the extension otherwise the files are not
+# read by doxygen.
+#
+# If left blank the following patterns are tested:*.c, *.cc, *.cxx, *.cpp,
+# *.c++, *.java, *.ii, *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h,
+# *.hh, *.hxx, *.hpp, *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc,
+# *.m, *.markdown, *.md, *.mm, *.dox, *.py, *.pyw, *.f90, *.f, *.for, *.tcl,
+# *.vhd, *.vhdl, *.ucf and *.qsf.
+
+FILE_PATTERNS = *.c \
+ *.cc \
+ *.cxx \
+ *.cpp \
+ *.c++ \
+ *.h \
+ *.hh \
+ *.hxx \
+ *.hpp \
+ *.h++ \
+ *.inc \
+ *.m \
+ *.markdown \
+ *.md \
+ *.mm \
+ *.dox \
+ *.f90 \
+ *.f \
+ *.for \
+ *.F90
# The RECURSIVE tag can be used to specify whether or not subdirectories should
# be searched for input files as well.
@@ -833,7 +879,7 @@ EXAMPLE_PATH =
# *.h) to filter out the source-files in the directories. If left blank all
# files are included.
-EXAMPLE_PATTERNS =
+EXAMPLE_PATTERNS = *
# If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be
# searched for input files to be used with the \include or \dontinclude commands
@@ -862,6 +908,10 @@ IMAGE_PATH =
# Note that the filter must not add or remove lines; it is applied before the
# code is scanned, but not when the output code is generated. If lines are added
# or removed, the anchors will not be placed correctly.
+#
+# Note that for custom extensions or not directly supported extensions you also
+# need to set EXTENSION_MAPPING for the extension otherwise the files are not
+# properly processed by doxygen.
INPUT_FILTER =
@@ -871,6 +921,10 @@ INPUT_FILTER =
# (like *.cpp=my_cpp_filter). See INPUT_FILTER for further information on how
# filters are used. If the FILTER_PATTERNS tag is empty or if none of the
# patterns match the file name, INPUT_FILTER is applied.
+#
+# Note that for custom extensions or not directly supported extensions you also
+# need to set EXTENSION_MAPPING for the extension otherwise the files are not
+# properly processed by doxygen.
FILTER_PATTERNS =
@@ -1129,6 +1183,7 @@ HTML_COLORSTYLE_GAMMA = 80
# If the HTML_TIMESTAMP tag is set to YES then the footer of each generated HTML
# page will contain the date and time when the page was generated. Setting this
+# to YES can help to show when doxygen was last run and thus if the
# to NO can help when comparing the output of multiple runs.
# The default value is: YES.
# This tag requires that the tag GENERATE_HTML is set to YES.
@@ -1604,9 +1659,12 @@ COMPACT_LATEX = NO
PAPER_TYPE = a4
# The EXTRA_PACKAGES tag can be used to specify one or more LaTeX package names
-# that should be included in the LaTeX output. To get the times font for
-# instance you can specify
-# EXTRA_PACKAGES=times
+# that should be included in the LaTeX output. The package can be specified just
+# by its name or with the correct syntax as to be used with the LaTeX
+# \usepackage command. To get the times font for instance you can specify :
+# EXTRA_PACKAGES=times or EXTRA_PACKAGES={times}
+# To use the option intlimits with the amsmath package you can specify:
+# EXTRA_PACKAGES=[intlimits]{amsmath}
# If left blank no extra packages will be included.
# This tag requires that the tag GENERATE_LATEX is set to YES.
@@ -1709,6 +1767,14 @@ LATEX_SOURCE_CODE = NO
LATEX_BIB_STYLE = plain
+# If the LATEX_TIMESTAMP tag is set to YES then the footer of each generated
+# page will contain the date and time when the page was generated. Setting this
+# to NO can help when comparing the output of multiple runs.
+# The default value is: NO.
+# This tag requires that the tag GENERATE_LATEX is set to YES.
+
+LATEX_TIMESTAMP = NO
+
#---------------------------------------------------------------------------
# Configuration options related to the RTF output
#---------------------------------------------------------------------------
@@ -1940,7 +2006,7 @@ ENABLE_PREPROCESSING = YES
# The default value is: NO.
# This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
-MACRO_EXPANSION = NO
+MACRO_EXPANSION = YES
# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES then
# the macro expansion is limited to the macros specified with the PREDEFINED and
@@ -1948,7 +2014,7 @@ MACRO_EXPANSION = NO
# The default value is: NO.
# This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
-EXPAND_ONLY_PREDEF = NO
+EXPAND_ONLY_PREDEF = YES
# If the SEARCH_INCLUDES tag is set to YES, the include files in the
# INCLUDE_PATH will be searched if a #include is found.
@@ -1980,7 +2046,7 @@ INCLUDE_FILE_PATTERNS =
# recursively expanded use the := operator instead of the = operator.
# This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
-PREDEFINED =
+PREDEFINED = ALLOCABLE_=allocatable
# If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
# tag can be used to specify a list of macro names that should be expanded. The
@@ -2207,7 +2273,8 @@ INCLUDED_BY_GRAPH = YES
#
# Note that enabling this option will significantly increase the time of a run.
# So in most cases it will be better to enable call graphs for selected
-# functions only using the \callgraph command.
+# functions only using the \callgraph command. Disabling a call graph can be
+# accomplished by means of the command \hidecallgraph.
# The default value is: NO.
# This tag requires that the tag HAVE_DOT is set to YES.
@@ -2218,7 +2285,8 @@ CALL_GRAPH = YES
#
# Note that enabling this option will significantly increase the time of a run.
# So in most cases it will be better to enable caller graphs for selected
-# functions only using the \callergraph command.
+# functions only using the \callergraph command. Disabling a caller graph can be
+# accomplished by means of the command \hidecallergraph.
# The default value is: NO.
# This tag requires that the tag HAVE_DOT is set to YES.
@@ -2241,11 +2309,15 @@ GRAPHICAL_HIERARCHY = YES
DIRECTORY_GRAPH = YES
# The DOT_IMAGE_FORMAT tag can be used to set the image format of the images
-# generated by dot.
+# generated by dot. For an explanation of the image formats see the section
+# output formats in the documentation of the dot tool (Graphviz (see:
+# http://www.graphviz.org/)).
# Note: If you choose svg you need to set HTML_FILE_EXTENSION to xhtml in order
# to make the SVG files visible in IE 9+ (other browsers do not have this
# requirement).
-# Possible values are: png, jpg, gif and svg.
+# Possible values are: png, jpg, gif, svg, png:gd, png:gd:gd, png:cairo,
+# png:cairo:gd, png:cairo:cairo, png:cairo:gdiplus, png:gdiplus and
+# png:gdiplus:gdiplus.
# The default value is: png.
# This tag requires that the tag HAVE_DOT is set to YES.
diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90
index b59cebd277..6649985423 100644
--- a/src/ALE/MOM_regridding.F90
+++ b/src/ALE/MOM_regridding.F90
@@ -2770,7 +2770,7 @@ subroutine regridding_memory_deallocation( CS )
end subroutine regridding_memory_deallocation
-!> \namespace MOM_regridding
+!> \namespace mom_regridding
!!
!! A vertical grid is defined solely by the cell thicknesses, \f$h\f$.
!! Most calculations in this module start with the coordinate at the bottom
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index e96888250a..e2f9998b6a 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -121,14 +121,15 @@ module MOM_forcing_type
! land ice-shelf related inputs
real, pointer, dimension(:,:) :: &
- ustar_shelf => NULL(), & !< friction velocity under ice-shelves (m/s)
- !! as computed by the ocean at the previous time step.
- frac_shelf_h => NULL(), & !< Fractional ice shelf coverage of h-, u-, and v-
- frac_shelf_u => NULL(), & !< cells, nondimensional from 0 to 1. These are only
- frac_shelf_v => NULL(), & !< associated if ice shelves are enabled, and are
- !! exactly 0 away from shelves or on land.
- rigidity_ice_u => NULL(),& !< Depth-integrated lateral viscosity of ice
- rigidity_ice_v => NULL() !< shelves or sea ice at u- or v-points (m3/s)
+ ustar_shelf => NULL(), & !< friction velocity under ice-shelves (m/s)
+ !! as computed by the ocean at the previous time step.
+ frac_shelf_h => NULL(), & !< Fractional ice shelf coverage of h-, u-, and v-
+ frac_shelf_u => NULL(), & !< cells, nondimensional from 0 to 1. These are only
+ frac_shelf_v => NULL(), & !< associated if ice shelves are enabled, and are
+ !! exactly 0 away from shelves or on land.
+ iceshelf_melt => NULL(), & !< ice shelf melt rate (positive) or freezing (negative) ( m/year )
+ rigidity_ice_u => NULL(),& !< Depth-integrated lateral viscosity of ice
+ rigidity_ice_v => NULL() !< shelves or sea ice at u- or v-points (m3/s)
! Scalars set by surface forcing modules
real :: vPrecGlobalAdj !< adjustment to restoring vprec to zero out global net ( kg/(m^2 s) )
@@ -1615,6 +1616,13 @@ subroutine forcing_accumulate(flux_tmp, fluxes, dt, G, wt2)
fluxes%ustar_shelf(i,j) = flux_tmp%ustar_shelf(i,j)
enddo ; enddo
endif
+
+ if (associated(fluxes%iceshelf_melt) .and. associated(flux_tmp%iceshelf_melt)) then
+ do i=isd,ied ; do j=jsd,jed
+ fluxes%iceshelf_melt(i,j) = flux_tmp%iceshelf_melt(i,j)
+ enddo ; enddo
+ endif
+
if (associated(fluxes%frac_shelf_h) .and. associated(flux_tmp%frac_shelf_h)) then
do i=isd,ied ; do j=jsd,jed
fluxes%frac_shelf_h(i,j) = flux_tmp%frac_shelf_h(i,j)
@@ -2198,6 +2206,7 @@ subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, p
call myAlloc(fluxes%frac_shelf_u,IsdB,IedB,jsd,jed, shelf)
call myAlloc(fluxes%frac_shelf_v,isd,ied,JsdB,JedB, shelf)
call myAlloc(fluxes%ustar_shelf,isd,ied,jsd,jed, shelf)
+ call myAlloc(fluxes%iceshelf_melt,isd,ied,jsd,jed, shelf)
call myAlloc(fluxes%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf)
call myAlloc(fluxes%rigidity_ice_v,isd,ied,JsdB,JedB, shelf)
@@ -2263,6 +2272,7 @@ subroutine deallocate_forcing_type(fluxes)
if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal)
if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal)
if (associated(fluxes%ustar_shelf)) deallocate(fluxes%ustar_shelf)
+ if (associated(fluxes%iceshelf_melt)) deallocate(fluxes%iceshelf_melt)
if (associated(fluxes%frac_shelf_h)) deallocate(fluxes%frac_shelf_h)
if (associated(fluxes%frac_shelf_u)) deallocate(fluxes%frac_shelf_u)
if (associated(fluxes%frac_shelf_v)) deallocate(fluxes%frac_shelf_v)
diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90
index d1e0644f57..ef2b752cd9 100644
--- a/src/ice_shelf/MOM_ice_shelf.F90
+++ b/src/ice_shelf/MOM_ice_shelf.F90
@@ -678,10 +678,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS)
else !not shelf
CS%t_flux(i,j) = 0.0
endif
-
enddo ! i-loop
enddo ! j-loop
+ ! melt in m/year
+ fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice)
+
if (CS%DEBUG) then
call hchksum (CS%h_shelf, "melt rate", G, haloshift=0)
endif
@@ -1382,6 +1384,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti
allocate( fluxes%frac_shelf_u(Isdq:Iedq,jsd:jed) ) ; fluxes%frac_shelf_u(:,:) = 0.0
allocate( fluxes%frac_shelf_v(isd:ied,Jsdq:Jedq) ) ; fluxes%frac_shelf_v(:,:) = 0.0
allocate( fluxes%ustar_shelf(isd:ied,jsd:jed) ) ; fluxes%ustar_shelf(:,:) = 0.0
+ allocate( fluxes%iceshelf_melt(isd:ied,jsd:jed) ) ; fluxes%iceshelf_melt(:,:) = 0.0
if (.not.associated(fluxes%p_surf)) then
allocate( fluxes%p_surf(isd:ied,jsd:jed) ) ; fluxes%p_surf(:,:) = 0.0
endif
@@ -1445,6 +1448,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti
if (.not. solo_ice_sheet) then
vd = var_desc("ustar_shelf","m s-1","Friction velocity under ice shelves",z_grid='1')
call register_restart_field(fluxes%ustar_shelf, vd, .true., CS%restart_CSp)
+ vd = var_desc("iceshelf_melt","m year-1","Ice Shelf Melt Rate",z_grid='1')
+ call register_restart_field(fluxes%iceshelf_melt, vd, .true., CS%restart_CSp)
endif
CS%restart_output_dir = dirs%restart_output_dir
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 0ec1bc5dbe..7f32001ace 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1,78 +1,7 @@
+!> Implements vertical viscosity (vertvisc)
+
module MOM_vert_friction
-!***********************************************************************
-!* GNU General Public License *
-!* This file is a part of MOM. *
-!* *
-!* MOM is free software; you can redistribute it and/or modify it and *
-!* are expected to follow the terms of the GNU General Public License *
-!* as published by the Free Software Foundation; either version 2 of *
-!* the License, or (at your option) any later version. *
-!* *
-!* MOM is distributed in the hope that it will be useful, but WITHOUT *
-!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
-!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public *
-!* License for more details. *
-!* *
-!* For the full text of the GNU General Public License, *
-!* write to: Free Software Foundation, Inc., *
-!* 675 Mass Ave, Cambridge, MA 02139, USA. *
-!* or see: http://www.gnu.org/licenses/gpl.html *
-!***********************************************************************
-
-!********+*********+*********+*********+*********+*********+*********+**
-!* *
-!* By Robert Hallberg, April 1994 - October 2006 *
-!* *
-!* This file contains the subroutine that implements vertical *
-!* viscosity (vertvisc). *
-!* *
-!* The vertical diffusion of momentum is fully implicit. This is *
-!* necessary to allow for vanishingly small layers. The coupling *
-!* is based on the distance between the centers of adjacent layers, *
-!* except where a layer is close to the bottom compared with a *
-!* bottom boundary layer thickness when a bottom drag law is used. *
-!* A stress top b.c. and a no slip bottom b.c. are used. There *
-!* is no limit on the time step for vertvisc. *
-!* *
-!* Near the bottom, the horizontal thickness interpolation scheme *
-!* changes to an upwind biased estimate to control the effect of *
-!* spurious Montgomery potential gradients at the bottom where *
-!* nearly massless layers layers ride over the topography. Within a *
-!* few boundary layer depths of the bottom, the harmonic mean *
-!* thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity *
-!* is from the thinner side and the arithmetic mean thickness *
-!* (i.e. (h+ + h-)/2) is used if the velocity is from the thicker *
-!* side. Both of these thickness estimates are second order *
-!* accurate. Above this the arithmetic mean thickness is used. *
-!* *
-!* In addition, vertvisc truncates any velocity component that *
-!* exceeds maxvel to truncvel. This basically keeps instabilities *
-!* spatially localized. The number of times the velocity is *
-!* truncated is reported each time the energies are saved, and if *
-!* exceeds CS%Maxtrunc the model will stop itself and change the time *
-!* to a large value. This has proven very useful in (1) diagnosing *
-!* model failures and (2) letting the model settle down to a *
-!* meaningful integration from a poorly specified initial condition. *
-!* *
-!* The same code is used for the two velocity components, by *
-!* indirectly referencing the velocities and defining a handful of *
-!* direction-specific defined variables. *
-!* *
-!* Macros written all in capital letters are defined in MOM_memory.h. *
-!* *
-!* A small fragment of the grid is shown below: *
-!* *
-!* j+1 x ^ x ^ x At x: q *
-!* j+1 > o > o > At ^: v, frhatv, tauy *
-!* j x ^ x ^ x At >: u, frhatu, taux *
-!* j > o > o > At o: h *
-!* j-1 x ^ x ^ x *
-!* i-1 i i+1 At x & ^: *
-!* i i+1 At > & o: *
-!* *
-!* The boundaries always run through q grid points (x). *
-!* *
-!********+*********+*********+*********+*********+*********+*********+**
+! This file is part of MOM6. See LICENSE.md for the license.
use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : diag_ctrl
@@ -100,137 +29,146 @@ module MOM_vert_friction
public updateCFLtruncationValue
type, public :: vertvisc_CS ; private
- real :: Hmix ! The mixed layer thickness in m.
- real :: Hmix_stress ! The mixed layer thickness over which the wind
- ! stress is applied with direct_stress, in m.
- real :: Kvml ! The mixed layer vertical viscosity in m2 s-1.
- real :: Kv ! The interior vertical viscosity in m2 s-1.
- real :: Hbbl ! The static bottom boundary layer thickness, in m.
- real :: Kvbbl ! The vertical viscosity in the bottom boundary
- ! layer, in m2 s-1.
-
- real :: maxvel ! Velocity components greater than maxvel,
- ! in m s-1, are truncated.
- logical :: CFL_based_trunc ! If true, base truncations on CFL numbers, not
- ! absolute velocities.
- real :: CFL_trunc ! Velocity components will be truncated when they
- ! are large enough that the corresponding CFL number
- ! exceeds this value, nondim.
- real :: CFL_report ! The value of the CFL number that will cause the
- ! accelerations to be reported, nondim. CFL_report
- ! will often equal CFL_trunc.
- real :: truncRampTime ! The time-scale over which to ramp up the value of
- ! CFL_trunc from zero to CFK_trunc0
- real :: CFL_truncS ! The start value of CFL_trunc
- real :: CFL_truncE ! The end/target value of CFL_trunc
- logical :: CFLrampingIsActivated = .false. ! True of the ramping has been initialized
- type(time_type) :: rampStartTime ! The time that the ramping of CFL_trunc starts
+ real :: Hmix !< The mixed layer thickness in m.
+ real :: Hmix_stress !< The mixed layer thickness over which the wind
+ !! stress is applied with direct_stress, in m.
+ real :: Kvml !< The mixed layer vertical viscosity in m2 s-1.
+ real :: Kv !< The interior vertical viscosity in m2 s-1.
+ real :: Hbbl !< The static bottom boundary layer thickness, in m.
+ real :: Kvbbl !< The vertical viscosity in the bottom boundary
+ !! layer, in m2 s-1.
+
+ real :: maxvel !< Velocity components greater than maxvel,
+ !! in m s-1, are truncated.
+ logical :: CFL_based_trunc !< If true, base truncations on CFL numbers, not
+ !! absolute velocities.
+ real :: CFL_trunc !< Velocity components will be truncated when they
+ !! are large enough that the corresponding CFL number
+ !! exceeds this value, nondim.
+ real :: CFL_report !< The value of the CFL number that will cause the
+ !! accelerations to be reported, nondim. CFL_report
+ !! will often equal CFL_trunc.
+ real :: truncRampTime !< The time-scale over which to ramp up the value of
+ !! CFL_trunc from CFL_truncS to CFL_truncE
+ real :: CFL_truncS !< The start value of CFL_trunc
+ real :: CFL_truncE !< The end/target value of CFL_trunc
+ logical :: CFLrampingIsActivated = .false. !< True if the ramping has been initialized
+ type(time_type) :: rampStartTime !< The time at which the ramping of CFL_trunc starts
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &
- a_u ! The u-drag coefficient across an interface, in m s-1.
+ a_u !< The u-drag coefficient across an interface, in m s-1.
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
- h_u ! The effective layer thickness at u-points, m or kg m-2.
+ h_u !< The effective layer thickness at u-points, m or kg m-2.
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &
- a_v ! The v-drag coefficient across an interface, in m s-1.
+ a_v !< The v-drag coefficient across an interface, in m s-1.
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
- h_v ! The effective layer thickness at v-points, m or kg m-2.
+ h_v !< The effective layer thickness at v-points, m or kg m-2.
+ !>@{
+ !! The surface coupling coefficient under ice shelves
+ !! in m s-1. Retained to determine stress under shelves.
real, pointer, dimension(:,:) :: &
- a1_shelf_u => NULL(), & ! The surface coupling coefficients under ice
- a1_shelf_v => NULL() ! shelves, in m s-1. These are retained to determine
- ! the stress under shelves.
-
- logical :: split ! If true, use the split time stepping scheme.
- logical :: bottomdraglaw ! If true, the bottom stress is calculated with a
- ! drag law c_drag*|u|*u. The velocity magnitude
- ! may be an assumed value or it may be based on the
- ! actual velocity in the bottommost HBBL, depending
- ! on whether linear_drag is true.
- logical :: Channel_drag ! If true, the drag is exerted directly on each
- ! layer according to what fraction of the bottom
- ! they overlie.
- logical :: harmonic_visc ! If true, the harmonic mean thicknesses are used
- ! to calculate the viscous coupling between layers
- ! except near the bottom. Otherwise the arithmetic
- ! mean thickness is used except near the bottom.
- real :: harm_BL_val ! A scale to determine when water is in the boundary
- ! layers based solely on harmonic mean thicknesses
- ! for the purpose of determining the extent to which
- ! the thicknesses used in the viscosities are upwinded.
- logical :: direct_stress ! If true, the wind stress is distributed over the
- ! topmost Hmix_stress of fluid and KVML may be very small.
- logical :: dynamic_viscous_ML ! If true, use the results from a dynamic
- ! calculation, perhaps based on a bulk Richardson
- ! number criterion, to determine the mixed layer
- ! thickness for viscosity.
- logical :: debug ! If true, write verbose checksums for debugging purposes.
- integer :: nkml ! The number of layers in the mixed layer.
- integer, pointer :: ntrunc ! The number of times the velocity has been
- ! truncated since the last call to write_energy.
- character(len=200) :: u_trunc_file ! The complete path to files in which a
- character(len=200) :: v_trunc_file ! column's worth of accelerations are
- ! written when velocity truncations occur.
- type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
- ! timing of diagnostic output.
+ a1_shelf_u => NULL(), &
+ a1_shelf_v => NULL()
+ !>@}
+
+ logical :: split !< If true, use the split time stepping scheme.
+ logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
+ !! drag law c_drag*|u|*u. The velocity magnitude
+ !! may be an assumed value or it may be based on the
+ !! actual velocity in the bottommost HBBL, depending
+ !! on whether linear_drag is true.
+ logical :: Channel_drag !< If true, the drag is exerted directly on each
+ !! layer according to what fraction of the bottom
+ !! they overlie.
+ logical :: harmonic_visc !< If true, the harmonic mean thicknesses are used
+ !! to calculate the viscous coupling between layers
+ !! except near the bottom. Otherwise the arithmetic
+ !! mean thickness is used except near the bottom.
+ real :: harm_BL_val !< A scale to determine when water is in the boundary
+ !! layers based solely on harmonic mean thicknesses
+ !! for the purpose of determining the extent to which
+ !! the thicknesses used in the viscosities are upwinded.
+ logical :: direct_stress !< If true, the wind stress is distributed over the
+ !! topmost Hmix_stress of fluid and KVML may be very small.
+ logical :: dynamic_viscous_ML !< If true, use the results from a dynamic
+ !! calculation, perhaps based on a bulk Richardson
+ !! number criterion, to determine the mixed layer
+ !! thickness for viscosity.
+ logical :: debug !< If true, write verbose checksums for debugging purposes.
+ integer :: nkml !< The number of layers in the mixed layer.
+ integer, pointer :: ntrunc !< The number of times the velocity has been
+ !! truncated since the last call to write_energy.
+ !>@{
+ !! The complete path to files in which a column's worth of
+ !! accelerations are written when velocity truncations occur.
+ character(len=200) :: u_trunc_file
+ character(len=200) :: v_trunc_file
+ !>@}
+
+ type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
+ !! timing of diagnostic output.
+
+ !>@{
+ !! Diagnostic identifiers
integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1
integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1
integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1
+ !>@}
type(PointAccel_CS), pointer :: PointAccel_CSp => NULL()
end type vertvisc_CS
contains
+!> Perform a fully implicit vertical diffusion
+!! of momentum. Stress top and bottom boundary conditions are used.
+!!
+!! This is solving the tridiagonal system
+!! \f[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1}
+!! = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \f]
+!! where \f$a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}\f$
+!! is the interfacial coupling thickness per time step,
+!! encompassing background viscosity as well as contributions from
+!! enhanced mixed and bottom layer viscosities.
+!! $r_k$ is a Rayleight drag term due to channel drag.
+!! There is an additional stress term on the right-hand side
+!! if DIRECT_STRESS is true, applied to the surface layer.
+
subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, &
taux_bot, tauy_bot)
-! This subroutine does a fully implicit vertical diffusion
-! of momentum. Stress top and bottom b.c.s are used.
- type(ocean_grid_type), intent(in) :: G
- type(verticalGrid_type), intent(in) :: GV
- real, intent(inout), dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u
- real, intent(inout), dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v
- real, intent(in), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h
- type(forcing), intent(in) :: fluxes
- type(vertvisc_type), intent(inout) :: visc
- real, intent(in) :: dt
- type(ocean_OBC_type), pointer :: OBC
- type(accel_diag_ptrs), intent(inout) :: ADp
- type(cont_diag_ptrs), intent(inout) :: CDp
- type(vertvisc_CS), pointer :: CS
- real, dimension(SZIB_(G),SZJ_(G)), optional, intent(out) :: taux_bot
- real, dimension(SZI_(G),SZJB_(G)), optional, intent(out) :: tauy_bot
-
-! Arguments: u - Zonal velocity, in m s-1. Intent in/out.
-! (in/out) v - Meridional velocity, in m s-1.
-! (in) h - Layer thickness, in m.
-! (in) fluxes - A structure containing pointers to any possible
-! forcing fields. Unused fields have NULL ptrs.
-! (in) visc - The vertical viscosity type, containing information about
-! viscosities and bottom drag-related quantities.
-! (in) dt - Time increment in s.
-! (in) OBC - This open boundary condition type specifies whether, where,
-! and what open boundary conditions are used.
-! (in) ADp - A structure pointing to the various accelerations in
-! the momentum equations, to enable the later calculation
-! of derived diagnostics, like energy budgets.
-! (in) CDp - A structure with pointers to various terms in the continuity
-! equations.
-! (in) G - The ocean's grid structure.
-! (in) CS - The control structure returned by a previous call to
-! vertvisc_init.
-! (out,opt) taux_bot - The zonal bottom stress from ocean to rock, in Pa.
-! (out,opt) tauy_bot - The meridional bottom stress from ocean to rock, in Pa.
-!
-! Fields from fluxes used in this subroutine:
-! taux: Zonal wind stress in Pa.
-! tauy: Meridional wind stress in Pa.
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ real, intent(inout), &
+ dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u !< Zonal velocity in m s-1
+ real, intent(inout), &
+ dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v !< Meridional velocity in m s-1
+ real, intent(in), &
+ dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< Layer thickness in H
+ type(forcing), intent(in) :: fluxes !< Structure containing forcing fields
+ type(vertvisc_type), intent(inout) :: visc !< Viscosities and bottom drag
+ real, intent(in) :: dt !< Time increment in s
+ type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
+ type(accel_diag_ptrs), intent(inout) :: ADp !< Accelerations in the momentum
+ !! equations for diagnostics
+ type(cont_diag_ptrs), intent(inout) :: CDp !< Continuity equation terms
+ type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
+ !> Zonal bottom stress from ocean to rock in Pa
+ real, optional, intent(out), dimension(SZIB_(G),SZJ_(G)) :: taux_bot
+ !> Meridional bottom stress from ocean to rock in Pa
+ real, optional, intent(out), dimension(SZI_(G),SZJB_(G)) :: tauy_bot
+
+ ! Fields from fluxes used in this subroutine:
+ ! taux: Zonal wind stress in Pa.
+ ! tauy: Meridional wind stress in Pa.
+
+ ! Local variables
real :: b1(SZIB_(G)) ! b1 and c1 are variables used by the
real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver. c1 is nondimensional,
! while b1 has units of inverse thickness.
real :: d1(SZIB_(G)) ! d1=1-c1 is used by the tridiagonal solver, ND.
- real :: Ray(SZIB_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity times the
- ! time step, in m.
- real :: b_denom_1 ! The first term in the denominator of b1, in m or kg m-2.
+ real :: Ray(SZIB_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity in m s-1
+ real :: b_denom_1 ! The first term in the denominator of b1, in H.
real :: Hmix ! The mixed layer thickness over which stress
! is applied with direct_stress, translated into
@@ -312,6 +250,29 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, &
Ray(I,k) = visc%Ray_u(I,j,k)
enddo ; enddo ; endif
+ ! perform forward elimination on the tridiagonal system
+ !
+ ! denote the diagonal of the system as b_i, the subdiagonal as a_i
+ ! and the superdiagonal as c_i. The right-hand side terms are d_i.
+ !
+ ! ignoring the rayleigh drag contribution,
+ ! we have a_i = -dt_m_to_H * a_u(i)
+ ! b_i = h_u(i) + dt_m_to_H * (a_u(i) + a_u(i+1))
+ ! c_i = -dt_m_to_H * a_u(i+1)
+ !
+ ! for forward elimination, we want to:
+ ! calculate c'_i = - c_i / (b_i + a_i c'_(i-1))
+ ! and d'_i = (d_i - a_i d'_(i-1)) / (b_i + a_i c'_(i-1))
+ ! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1
+ ! (see Thomas' tridiagonal matrix algorithm)
+ !
+ ! b1 is the denominator term 1 / (b_i + a_i c'_(i-1))
+ ! b_denom_1 is (b_i + a_i + c_i) - a_i(1 - c'_(i-1))
+ ! = (b_i + c_i + c'_(i-1))
+ ! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(i-1)
+ ! c1(i) is -c'_(i - 1)
+ ! and the right-hand-side is destructively updated to be d'_i
+ !
do I=Isq,Ieq ; if (do_i(I)) then
b_denom_1 = CS%h_u(I,j,1) + dt_m_to_H * (Ray(I,1) + CS%a_u(I,j,1))
b1(I) = 1.0 / (b_denom_1 + dt_m_to_H*CS%a_u(I,j,2))
@@ -326,6 +287,9 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, &
u(I,j,k) = (CS%h_u(I,j,k) * u(I,j,k) + &
dt_m_to_H * CS%a_u(I,j,K) * u(I,j,k-1)) * b1(I)
endif ; enddo ; enddo
+
+ ! back substitute to solve for the new velocities
+ ! u_i = d'_i - c'_i x_(i+1)
do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
u(I,j,k) = u(I,j,k) + c1(I,k+1) * u(I,j,k+1)
endif ; enddo ; enddo ! i and k loops
@@ -450,31 +414,24 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, &
end subroutine vertvisc
+!> Calculate the fraction of momentum originally in a layer that remains
+!! after a time-step of viscosity, and the fraction of a time-step's
+!! worth of barotropic acceleration that a layer experiences after
+!! viscosity is applied.
subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, CS)
-! This subroutine does a fully implicit vertical diffusion
-! of momentum. Stress top and bottom b.c.s are used.
- type(ocean_grid_type), intent(in) :: G
- type(verticalGrid_type), intent(in) :: GV
- type(vertvisc_type), intent(in) :: visc
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
+ !> Fraction of a time-step's worth of a barotopic acceleration that
+ !! a layer experiences after viscosity is applied in the zonal direction
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: visc_rem_u
+ !> Fraction of a time-step's worth of a barotopic acceleration that
+ !! a layer experiences after viscosity is applied in the meridional direction
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: visc_rem_v
- real, intent(in) :: dt
- type(vertvisc_CS), pointer :: CS
-! Arguments: visc - The vertical viscosity type, containing information about
-! viscosities and bottom drag-related quantities, intent in.
-! (out) visc_rem_u - Both the fraction of the momentum originally in a
-! (out) visc_rem_v - layer that remains after a time-step of viscosity,
-! and the fraction of a time-step's worth of a
-! barotropic acceleration that a layer experiences
-! after viscosity is applied, in the zonal (_u) and
-! meridional (_v) directions. Nondimensional between
-! 0 (at the bottom) and 1 (far above the bottom).
-! (in) dt - Time increment in s.
-! (in) G - The ocean's grid structure.
-! (in) GV - The ocean's vertical grid structure.
-! (in) CS - The control structure returned by a previous call to
-! vertvisc_init.
-!
+ real, intent(in) :: dt !< Time increment in s
+ type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
+
+ ! Local variables
real :: b1(SZIB_(G)) ! b1 and c1 are variables used by the
real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver. c1 is nondimensional,
@@ -566,37 +523,29 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, CS)
end subroutine vertvisc_remnant
+!> Calculate the coupling coefficients (CS%a_u and CS%a_v)
+!! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the
+!! applying the implicit vertical viscosity via vertvisc().
subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS)
-! This subroutine calculates the coupling coefficients (CS%a_u and CS%a_v)
-! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the
-! applying the implicit vertical viscosity via vertvisc.
- type(ocean_grid_type), intent(in) :: G
- type(verticalGrid_type), intent(in) :: GV
- real, intent(in), dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u
- real, intent(in), dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v
- real, intent(in), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h
- type(forcing), intent(in) :: fluxes
- type(vertvisc_type), intent(in) :: visc
- real, intent(in) :: dt
- type(vertvisc_CS), pointer :: CS
-
-! Arguments: u - Zonal velocity, in m s-1. Intent in.
-! (in) v - Meridional velocity, in m s-1.
-! (in) h - Layer thickness, in H.
-! (in) fluxes - A structure containing pointers to any possible
-! forcing fields. Unused fields have NULL ptrs.
-! (in) visc - The vertical viscosity type, containing information about
-! viscosities and bottom drag-related quantities.
-! (in) dt - Time increment in s.
-! (in) G - The ocean's grid structure.
-! (in) GV - The ocean's vertical grid structure.
-! (in/out) CS - The control structure returned by a previous call to
-! vertvisc_init.
-!
-! Field from fluxes used in this subroutine:
-! ustar: the friction velocity in m s-1, used here as the mixing
-! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
-!
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ real, intent(in), &
+ dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u !< Zonal velocity in m s-1
+ real, intent(in), &
+ dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v !< Meridional velocity in m s-1
+ real, intent(in), &
+ dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< Layer thickness in H
+ type(forcing), intent(in) :: fluxes !< Structure containing forcing fields
+ type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
+ real, intent(in) :: dt !< Time increment in s
+ type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
+
+ ! Field from fluxes used in this subroutine:
+ ! ustar: the friction velocity in m s-1, used here as the mixing
+ ! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
+
+ ! Local variables
+
real, dimension(SZIB_(G),SZK_(G)) :: &
h_harm, & ! Harmonic mean of the thicknesses around a velocity grid point,
! given by 2*(h+ * h-)/(h+ + h-), in m or kg m-2 (H for short).
@@ -973,39 +922,48 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS)
end subroutine vertvisc_coef
+!> Calculate the 'coupling coefficient' (a[k]) at the
+!! interfaces. If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the
+!! adjacent layer thicknesses are used to calculate a[k] near the bottom.
subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
dt, j, G, GV, CS, visc, fluxes, work_on_u, shelf)
-! This subroutine calculates the 'coupling coefficient' (a[k]) at the
-! interfaces. If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the
-! adjacent layer thicknesses are used to calculate a[k] near the bottom.
- type(ocean_grid_type), intent(in) :: G
- type(verticalGrid_type), intent(in) :: GV
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ !> Coupling coefficient across interfaces, in m s-1
real, dimension(SZIB_(G),SZK_(GV)+1), intent(out) :: a
+ !> Thickness at velocity points, in H
real, dimension(SZIB_(G),SZK_(GV)), intent(in) :: hvel
+ !> If true, determine coupling coefficient for a column
logical, dimension(SZIB_(G)), intent(in) :: do_i
+ !> Harmonic mean of thicknesses around a velocity grid point, in H
real, dimension(SZIB_(G),SZK_(GV)), intent(in) :: h_harm
- real, dimension(SZIB_(G)), intent(in) :: bbl_thick, kv_bbl
+ !> Bottom boundary layer thickness, in H
+ real, dimension(SZIB_(G)), intent(in) :: bbl_thick
+ !> Bottom boundary layer viscosity, in m2 s-1
+ real, dimension(SZIB_(G)), intent(in) :: kv_bbl
+ !> Estimate of interface heights above the bottom,
+ !! normalised by the bottom boundary layer thickness
real, dimension(SZIB_(G),SZK_(GV)+1), intent(in) :: z_i
+ !> Mixed layer depth, in H
real, dimension(SZIB_(G)), intent(out) :: h_ml
+ !> j-index to find coupling coefficient for
integer, intent(in) :: j
+ !> Time increment, in s
real, intent(in) :: dt
+ !> Vertical viscosity control structure
type(vertvisc_CS), pointer :: CS
+ !> Structure containing viscosities and bottom drag
type(vertvisc_type), intent(in) :: visc
+ !> Structure containing forcing fields
type(forcing), intent(in) :: fluxes
+ !> If true, u-points are being calculated, otherwise v-points
logical, intent(in) :: work_on_u
+ !> If present and true, use a surface boundary condition
+ !! appropriate for an ice shelf.
logical, optional, intent(in) :: shelf
-! Arguments: a - The coupling coefficent across interfaces, in m/s. Intent out.
-! (in) hvel - The thickness at velocity points, in H.
-! (in) do_i - If true, determine the a for a column.
-! ...
-! (in) dt - The amount of time covered by this call, in s.
-! (in) G - The ocean's grid structure.
-! (in) GV - The ocean's vertical grid structure.
-! ...
-! (in) work_on_u - If true, u-points are being worked on, otherwise this
-! call is for v-points.
-! (in) shelf - If present and true, use a surface boundary condition
-! appropriate for an ice shelf.
+
+ ! Local variables
+
real, dimension(SZIB_(G)) :: &
u_star, & ! ustar at a velocity point, in m s-1.
absf, & ! The average of the neighboring absolute values of f, in s-1.
@@ -1196,22 +1154,26 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m
end subroutine find_coupling_coef
-
+!> Velocity components which exceed a threshold for physically
+!! reasonable values are truncated. Optionally, any column with excessive
+!! velocities may be sent to a diagnostic reporting subroutine.
subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, GV, CS)
-! Within this subroutine, velocity components which exceed a threshold for
-! physically reasonable values are truncated. Optionally, any column with
-! excessive velocities may be sent to a diagnostic reporting subroutine.
- type(ocean_grid_type), intent(in) :: G
- type(verticalGrid_type), intent(in) :: GV
- real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u
- real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v
- real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h
- type(accel_diag_ptrs), intent(in) :: ADp
- type(cont_diag_ptrs), intent(in) :: CDp
- type(forcing), intent(in) :: fluxes
- type(vertvisc_type), intent(in) :: visc
- real, intent(in) :: dt
- type(vertvisc_CS), pointer :: CS
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ real, intent(inout), &
+ dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u !< Zonal velocity in m s-1
+ real, intent(inout), &
+ dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v !< Meridional velocity in m s-1
+ real, intent(in), &
+ dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< Layer thickness in H
+ type(accel_diag_ptrs), intent(in) :: ADp !< Acceleration diagnostic pointers
+ type(cont_diag_ptrs), intent(in) :: CDp !< Continuity diagnostic pointers
+ type(forcing), intent(in) :: fluxes !< Forcing fields
+ type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
+ real, intent(in) :: dt !< Time increment in s
+ type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
+
+ ! Local variables
real :: maxvel ! Velocities components greater than maxvel
real :: truncvel ! are truncated to truncvel, both in m s-1.
@@ -1394,34 +1356,23 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, GV, CS)
end subroutine vertvisc_limit_vel
+!> Initialise the vertical friction module
subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, &
ntrunc, CS)
+ !> "MOM Internal State", a set of pointers to the fields and accelerations
+ !! that make up the ocean's physical state
type(ocean_internal_state), target, intent(in) :: MIS
- type(time_type), target, intent(in) :: Time
- type(ocean_grid_type), intent(in) :: G
- type(verticalGrid_type), intent(in) :: GV
- type(param_file_type), intent(in) :: param_file
- type(diag_ctrl), target, intent(inout) :: diag
- type(accel_diag_ptrs), intent(inout) :: ADp
- type(directories), intent(in) :: dirs
- integer, target, intent(inout) :: ntrunc
- type(vertvisc_CS), pointer :: CS
-! Arguments: MIS - For "MOM Internal State" a set of pointers to the fields and
-! accelerations that make up the ocean's physical state.
-! (in) Time - The current model time.
-! (in) G - The ocean's grid structure.
-! (in) GV - The ocean's vertical grid structure.
-! (in) param_file - A structure indicating the open file to parse for
-! model parameter values.
-! (in) diag - A structure that is used to regulate diagnostic output.
-! (inout) ADp - A structure pointing to the various accelerations in
-! the momentum equations, to enable the later calculation
-! of derived diagnostics, like energy budgets.
-! (in) dirs - A structure containing several relevant directory paths.
-! (in/out) ntrunc - The integer that stores the number of times the velocity
-! has been truncated since the last call to write_energy.
-! (in/out) CS - A pointer that is set to point to the control structure
-! for this module
+ type(time_type), target, intent(in) :: Time !< Current model time
+ type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
+ type(param_file_type), intent(in) :: param_file !< File to parse for parameters
+ type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic control structure
+ type(accel_diag_ptrs), intent(inout) :: ADp !< Acceleration diagnostic pointers
+ type(directories), intent(in) :: dirs !< Relevant directory paths
+ integer, target, intent(inout) :: ntrunc !< Number of velocity truncations
+ type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
+
+ ! Local variables
real :: hmix_str_dflt
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
@@ -1583,13 +1534,15 @@ subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, &
end subroutine vertvisc_init
+!> Update the CFL truncation value as a function of time.
+!! If called with the optional argument activate=.true., record the
+!! value of Time as the beginning of the ramp period.
subroutine updateCFLtruncationValue(Time, CS, activate)
- ! This routine updates the CFL truncation value as a function of time
- ! If called with the optional argument activate=.true., records the
- ! value of Time as the beginning of the ramp period.
- type(time_type), target, intent(in) :: Time
- type(vertvisc_CS), pointer :: CS
+ type(time_type), target, intent(in) :: Time !< Current model time
+ type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
+ !> Whether to record the value of Time as the beginning of the ramp period
logical, optional, intent(in) :: activate
+
! Local variables
real :: deltaTime, wghtA
character(len=12) :: msg
@@ -1621,6 +1574,7 @@ subroutine updateCFLtruncationValue(Time, CS, activate)
" limit to "//trim(msg))
end subroutine updateCFLtruncationValue
+!> Clean up and deallocate the vertical friction module
subroutine vertvisc_end(CS)
type(vertvisc_CS), pointer :: CS
DEALLOC_(CS%a_u) ; DEALLOC_(CS%h_u)
@@ -1630,4 +1584,54 @@ subroutine vertvisc_end(CS)
deallocate(CS)
end subroutine vertvisc_end
+!> \namespace mom_vert_friction
+!! \author Robert Hallberg
+!! \date April 1994 - October 2006
+!!
+!! The vertical diffusion of momentum is fully implicit. This is
+!! necessary to allow for vanishingly small layers. The coupling
+!! is based on the distance between the centers of adjacent layers,
+!! except where a layer is close to the bottom compared with a
+!! bottom boundary layer thickness when a bottom drag law is used.
+!! A stress top b.c. and a no slip bottom b.c. are used. There
+!! is no limit on the time step for vertvisc.
+!!
+!! Near the bottom, the horizontal thickness interpolation scheme
+!! changes to an upwind biased estimate to control the effect of
+!! spurious Montgomery potential gradients at the bottom where
+!! nearly massless layers layers ride over the topography. Within a
+!! few boundary layer depths of the bottom, the harmonic mean
+!! thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity
+!! is from the thinner side and the arithmetic mean thickness
+!! (i.e. (h+ + h-)/2) is used if the velocity is from the thicker
+!! side. Both of these thickness estimates are second order
+!! accurate. Above this the arithmetic mean thickness is used.
+!!
+!! In addition, vertvisc truncates any velocity component that
+!! exceeds maxvel to truncvel. This basically keeps instabilities
+!! spatially localized. The number of times the velocity is
+!! truncated is reported each time the energies are saved, and if
+!! exceeds CS%Maxtrunc the model will stop itself and change the time
+!! to a large value. This has proven very useful in (1) diagnosing
+!! model failures and (2) letting the model settle down to a
+!! meaningful integration from a poorly specified initial condition.
+!!
+!! The same code is used for the two velocity components, by
+!! indirectly referencing the velocities and defining a handful of
+!! direction-specific defined variables.
+!!
+!! Macros written all in capital letters are defined in MOM_memory.h.
+!!
+!! A small fragment of the grid is shown below:
+!! \verbatim
+!! j+1 x ^ x ^ x At x: q
+!! j+1 > o > o > At ^: v, frhatv, tauy
+!! j x ^ x ^ x At >: u, frhatu, taux
+!! j > o > o > At o: h
+!! j-1 x ^ x ^ x
+!! i-1 i i+1 At x & ^:
+!! i i+1 At > & o:
+!! \endverbatim
+!!
+!! The boundaries always run through q grid points (x).
end module MOM_vert_friction
diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90
new file mode 100644
index 0000000000..072f5b4d87
--- /dev/null
+++ b/src/tracer/ISOMIP_tracer.F90
@@ -0,0 +1,425 @@
+!> This module contains the routines used to set up and use a set of (one for now)
+!! dynamically passive tracers. For now, just one passive tracer is injected in
+!! the sponge layer.
+!! Set up and use passive tracers requires the following:
+!! (1) register_ISOMIP_tracer
+!! (2)
+
+!********+*********+*********+*********+*********+*********+*********+**
+!* *
+!* By Robert Hallberg, 2002 *
+!* Adapted to the ISOMIP test case by Gustavo Marques, May 2016 * !* *
+!********+*********+*********+*********+*********+*********+*********+**
+
+module ISOMIP_tracer
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
+use MOM_diag_mediator, only : diag_ctrl
+use MOM_diag_to_Z, only : register_Z_tracer, diag_to_Z_CS
+use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
+use MOM_forcing_type, only : forcing
+use MOM_hor_index, only : hor_index_type
+use MOM_grid, only : ocean_grid_type
+use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc
+use MOM_restart, only : register_restart_field, MOM_restart_CS
+use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
+use MOM_time_manager, only : time_type, get_time
+use MOM_tracer_registry, only : register_tracer, tracer_registry_type
+use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
+use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_variables, only : surface, ocean_OBC_type
+use MOM_verticalGrid, only : verticalGrid_type
+
+use coupler_util, only : set_coupler_values, ind_csurf
+use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux
+
+implicit none ; private
+
+#include
+
+!< Publicly available functions
+public register_ISOMIP_tracer, initialize_ISOMIP_tracer
+public ISOMIP_tracer_column_physics, ISOMIP_tracer_surface_state, ISOMIP_tracer_end
+
+!< ntr is the number of tracers in this module.
+integer, parameter :: ntr = 2
+
+type p3d
+ real, dimension(:,:,:), pointer :: p => NULL()
+end type p3d
+
+!> tracer control structure
+type, public :: ISOMIP_tracer_CS ; private
+ logical :: coupled_tracers = .false. !< These tracers are not offered to the
+ !< coupler.
+ character(len = 200) :: tracer_IC_file !< The full path to the IC file, or " "
+ !< to initialize internally.
+ type(time_type), pointer :: Time !< A pointer to the ocean model's clock.
+ type(tracer_registry_type), pointer :: tr_Reg => NULL()
+ real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this
+ !< subroutine, in g m-3?
+ real, pointer :: tr_aux(:,:,:,:) => NULL() !< The masked tracer concentration
+ !< for output, in g m-3.
+ type(p3d), dimension(NTR) :: &
+ tr_adx, &!< Tracer zonal advective fluxes in g m-3 m3 s-1.
+ tr_ady, &!< Tracer meridional advective fluxes in g m-3 m3 s-1.
+ tr_dfx, &!< Tracer zonal diffusive fluxes in g m-3 m3 s-1.
+ tr_dfy !< Tracer meridional diffusive fluxes in g m-3 m3 s-1.
+ real :: land_val(NTR) = -1.0 !< The value of tr used where land is masked out.
+ logical :: mask_tracers !< If true, tracers are masked out in massless layers.
+ logical :: use_sponge
+
+ integer, dimension(NTR) :: ind_tr !< Indices returned by aof_set_coupler_flux
+ !< if it is used and the surface tracer concentrations are to be
+ !< provided to the coupler.
+
+ type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
+ !< timing of diagnostic output.
+ integer, dimension(NTR) :: id_tracer = -1, id_tr_adx = -1, id_tr_ady = -1
+ integer, dimension(NTR) :: id_tr_dfx = -1, id_tr_dfy = -1
+
+ type(vardesc) :: tr_desc(NTR)
+end type ISOMIP_tracer_CS
+
+contains
+
+
+!> This subroutine is used to register tracer fields
+function register_ISOMIP_tracer(HI, GV, param_file, CS, tr_Reg, &
+ restart_CS)
+ type(hor_index_type), intent(in) :: HI ! NULL()
+ logical :: register_ISOMIP_tracer
+ integer :: isd, ied, jsd, jed, nz, m
+ isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
+
+ if (associated(CS)) then
+ call MOM_error(WARNING, "ISOMIP_register_tracer called with an "// &
+ "associated control structure.")
+ return
+ endif
+ allocate(CS)
+
+ ! Read all relevant parameters and write them to the model log.
+ call log_version(param_file, mod, version, "")
+ call get_param(param_file, mod, "ISOMIP_TRACER_IC_FILE", CS%tracer_IC_file, &
+ "The name of a file from which to read the initial \n"//&
+ "conditions for the ISOMIP tracers, or blank to initialize \n"//&
+ "them internally.", default=" ")
+ if (len_trim(CS%tracer_IC_file) >= 1) then
+ call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
+ inputdir = slasher(inputdir)
+ CS%tracer_IC_file = trim(inputdir)//trim(CS%tracer_IC_file)
+ call log_param(param_file, mod, "INPUTDIR/ISOMIP_TRACER_IC_FILE", &
+ CS%tracer_IC_file)
+ endif
+ call get_param(param_file, mod, "SPONGE", CS%use_sponge, &
+ "If true, sponges may be applied anywhere in the domain. \n"//&
+ "The exact location and properties of those sponges are \n"//&
+ "specified from MOM_initialization.F90.", default=.false.)
+
+ allocate(CS%tr(isd:ied,jsd:jed,nz,NTR)) ; CS%tr(:,:,:,:) = 0.0
+ if (CS%mask_tracers) then
+ allocate(CS%tr_aux(isd:ied,jsd:jed,nz,NTR)) ; CS%tr_aux(:,:,:,:) = 0.0
+ endif
+
+ do m=1,NTR
+ if (m < 10) then ; write(name,'("tr_D",I1.1)') m
+ else ; write(name,'("tr_D",I2.2)') m ; endif
+ write(longname,'("Concentration of ISOMIP Tracer ",I2.2)') m
+ CS%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mod)
+
+ ! This is needed to force the compiler not to do a copy in the registration
+ ! calls. Curses on the designers and implementers of Fortran90.
+ tr_ptr => CS%tr(:,:,:,m)
+ ! Register the tracer for the restart file.
+ call register_restart_field(tr_ptr, CS%tr_desc(m), .true., restart_CS)
+ ! Register the tracer for horizontal advection & diffusion.
+ call register_tracer(tr_ptr, CS%tr_desc(m), param_file, HI, GV, tr_Reg, &
+ tr_desc_ptr=CS%tr_desc(m))
+
+ ! Set coupled_tracers to be true (hard-coded above) to provide the surface
+ ! values to the coupler (if any). This is meta-code and its arguments will
+ ! currently (deliberately) give fatal errors if it is used.
+ if (CS%coupled_tracers) &
+ CS%ind_tr(m) = aof_set_coupler_flux(trim(name)//'_flux', &
+ flux_type=' ', implementation=' ', caller="register_ISOMIP_tracer")
+ enddo
+
+ CS%tr_Reg => tr_Reg
+ register_ISOMIP_tracer = .true.
+end function register_ISOMIP_tracer
+
+!> Initializes the NTR tracer fields in tr(:,:,:,:)
+! and it sets up the tracer output.
+subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS, ALE_sponge_CSp, &
+ diag_to_Z_CSp)
+ type(ocean_grid_type), intent(in) :: G !< Grid structure.
+ type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
+ logical, intent(in) :: restart !< .true. if the fields have already been read from a restart file.
+ type(time_type), target, intent(in) :: day !< Time of the start of the run.
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2.
+ type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, where, and what open boundary conditions are used. This is not being used for now.
+ type(ISOMIP_tracer_CS), pointer :: CS !< The control structure returned by a previous call to ISOMIP_register_tracer.
+ type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< A pointer to the control structure for the sponges, if they are in use. Otherwise this may be unassociated.
+ type(diag_to_Z_CS), pointer :: diag_to_Z_CSp !< A pointer to the control structure for diagnostics in depth space.
+
+ real, allocatable :: temp(:,:,:)
+ real, pointer, dimension(:,:,:) :: &
+ OBC_tr1_u => NULL(), & ! These arrays should be allocated and set to
+ OBC_tr1_v => NULL() ! specify the values of tracer 1 that should come
+ ! in through u- and v- points through the open
+ ! boundary conditions, in the same units as tr.
+ character(len=16) :: name ! A variable's name in a NetCDF file.
+ character(len=72) :: longname ! The long name of that variable.
+ character(len=48) :: units ! The dimensions of the variable.
+ character(len=48) :: flux_units ! The units for tracer fluxes, usually
+ ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
+ real, pointer :: tr_ptr(:,:,:) => NULL()
+ real :: PI ! 3.1415926... calculated as 4*atan(1)
+ real :: tr_y ! Initial zonally uniform tracer concentrations.
+ real :: dist2 ! The distance squared from a line, in m2.
+ real :: h_neglect ! A thickness that is so small it is usually lost
+ ! in roundoff and can be neglected, in m.
+ real :: e(SZK_(G)+1), e_top, e_bot, d_tr
+ integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
+ integer :: IsdB, IedB, JsdB, JedB
+
+ if (.not.associated(CS)) return
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+ isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
+ IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
+ h_neglect = GV%H_subroundoff
+
+ CS%Time => day
+
+ if (.not.restart) then
+ if (len_trim(CS%tracer_IC_file) >= 1) then
+ ! Read the tracer concentrations from a netcdf file.
+ if (.not.file_exists(CS%tracer_IC_file, G%Domain)) &
+ call MOM_error(FATAL, "ISOMIP_initialize_tracer: Unable to open "// &
+ CS%tracer_IC_file)
+ do m=1,NTR
+ call query_vardesc(CS%tr_desc(m), name, caller="initialize_ISOMIP_tracer")
+ call read_data(CS%tracer_IC_file, trim(name), &
+ CS%tr(:,:,:,m), domain=G%Domain%mpp_domain)
+ enddo
+ else
+ do m=1,NTR
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ CS%tr(i,j,k,m) = 1.0e-20 ! This could just as well be 0.
+ enddo ; enddo ; enddo
+ enddo
+ endif
+ endif ! restart
+
+ if ( CS%use_sponge ) then
+! If sponges are used, this example damps tracers in sponges in the
+! northern half of the domain to 1 and tracers in the southern half
+! to 0. For any tracers that are not damped in the sponge, the call
+! to set_up_sponge_field can simply be omitted.
+ if (.not.associated(ALE_sponge_CSp)) &
+ call MOM_error(FATAL, "ISOMIP_initialize_tracer: "// &
+ "The pointer to ALEsponge_CSp must be associated if SPONGE is defined.")
+
+ allocate(temp(G%isd:G%ied,G%jsd:G%jed,nz))
+
+ do j=js,je ; do i=is,ie
+ if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then
+ temp(i,j,:) = 1.0
+ else
+ temp(i,j,:) = 0.0
+ endif
+ enddo ; enddo
+
+! do m=1,NTR
+ do m=1,1
+ ! This is needed to force the compiler not to do a copy in the sponge
+ ! calls. Curses on the designers and implementers of Fortran90.
+ tr_ptr => CS%tr(:,:,:,m)
+ call set_up_ALE_sponge_field(temp, G, tr_ptr, ALE_sponge_CSp)
+ enddo
+ deallocate(temp)
+ endif
+
+ ! This needs to be changed if the units of tracer are changed above.
+ if (GV%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
+ else ; flux_units = "kg s-1" ; endif
+
+ do m=1,NTR
+ ! Register the tracer for the restart file.
+ call query_vardesc(CS%tr_desc(m), name, units=units, longname=longname, &
+ caller="initialize_ISOMIP_tracer")
+ CS%id_tracer(m) = register_diag_field("ocean_model", trim(name), CS%diag%axesTL, &
+ day, trim(longname) , trim(units))
+ CS%id_tr_adx(m) = register_diag_field("ocean_model", trim(name)//"_adx", &
+ CS%diag%axesCuL, day, trim(longname)//" advective zonal flux" , &
+ trim(flux_units))
+ CS%id_tr_ady(m) = register_diag_field("ocean_model", trim(name)//"_ady", &
+ CS%diag%axesCvL, day, trim(longname)//" advective meridional flux" , &
+ trim(flux_units))
+ CS%id_tr_dfx(m) = register_diag_field("ocean_model", trim(name)//"_dfx", &
+ CS%diag%axesCuL, day, trim(longname)//" diffusive zonal flux" , &
+ trim(flux_units))
+ CS%id_tr_dfy(m) = register_diag_field("ocean_model", trim(name)//"_dfy", &
+ CS%diag%axesCvL, day, trim(longname)//" diffusive zonal flux" , &
+ trim(flux_units))
+ if (CS%id_tr_adx(m) > 0) call safe_alloc_ptr(CS%tr_adx(m)%p,IsdB,IedB,jsd,jed,nz)
+ if (CS%id_tr_ady(m) > 0) call safe_alloc_ptr(CS%tr_ady(m)%p,isd,ied,JsdB,JedB,nz)
+ if (CS%id_tr_dfx(m) > 0) call safe_alloc_ptr(CS%tr_dfx(m)%p,IsdB,IedB,jsd,jed,nz)
+ if (CS%id_tr_dfy(m) > 0) call safe_alloc_ptr(CS%tr_dfy(m)%p,isd,ied,JsdB,JedB,nz)
+
+! Register the tracer for horizontal advection & diffusion.
+ if ((CS%id_tr_adx(m) > 0) .or. (CS%id_tr_ady(m) > 0) .or. &
+ (CS%id_tr_dfx(m) > 0) .or. (CS%id_tr_dfy(m) > 0)) &
+ call add_tracer_diagnostics(name, CS%tr_Reg, CS%tr_adx(m)%p, &
+ CS%tr_ady(m)%p, CS%tr_dfx(m)%p, CS%tr_dfy(m)%p)
+
+ call register_Z_tracer(CS%tr(:,:,:,m), trim(name), longname, units, &
+ day, G, diag_to_Z_CSp)
+ enddo
+
+end subroutine initialize_ISOMIP_tracer
+
+!> This subroutine applies diapycnal diffusion and any other column
+! tracer physics or chemistry to the tracers from this file.
+! This is a simple example of a set of advected passive tracers.
+subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+ type(ocean_grid_type), intent(in) :: G
+ type(verticalGrid_type), intent(in) :: GV
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
+ type(forcing), intent(in) :: fluxes
+ real, intent(in) :: dt
+ type(ISOMIP_tracer_CS), pointer :: CS
+
+! Arguments: h_old - Layer thickness before entrainment, in m or kg m-2.
+! (in) h_new - Layer thickness after entrainment, in m or kg m-2.
+! (in) ea - an array to which the amount of fluid entrained
+! from the layer above during this call will be
+! added, in m or kg m-2.
+! (in) eb - an array to which the amount of fluid entrained
+! from the layer below during this call will be
+! added, in m or kg m-2.
+! (in) fluxes - A structure containing pointers to any possible
+! forcing fields. Unused fields have NULL ptrs.
+! (in) dt - The amount of time covered by this call, in s.
+! (in) G - The ocean's grid structure.
+! (in) GV - The ocean's vertical grid structure.
+! (in) CS - The control structure returned by a previous call to
+! ISOMIP_register_tracer.
+!
+! The arguments to this subroutine are redundant in that
+! h_new[k] = h_old[k] + ea[k] - eb[k-1] + eb[k] - ea[k+1]
+
+ real :: b1(SZI_(G)) ! b1 and c1 are variables used by the
+ real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
+ integer :: i, j, k, is, ie, js, je, nz, m
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+
+ if (.not.associated(CS)) return
+
+
+ ! dye melt water (m=2)
+ ! converting melt (m/year) to (m/s)
+ do m=2,NTR
+ do j=js,je ; do i=is,ie
+ if (fluxes%iceshelf_melt(i,j) > 0.0) then
+! CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) &
+! + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1)
+ CS%tr(i,j,1,m) = 1.0
+ endif
+ enddo ; enddo
+ enddo
+
+ do m=1,NTR
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+
+ if (CS%mask_tracers) then
+ do m = 1,NTR ; if (CS%id_tracer(m) > 0) then
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ if (h_new(i,j,k) < 1.1*GV%Angstrom) then
+ CS%tr_aux(i,j,k,m) = CS%land_val(m)
+ else
+ CS%tr_aux(i,j,k,m) = CS%tr(i,j,k,m)
+ endif
+ enddo ; enddo ; enddo
+ endif ; enddo
+ endif
+
+ do m=1,NTR
+ if (CS%mask_tracers) then
+ if (CS%id_tracer(m)>0) &
+ call post_data(CS%id_tracer(m),CS%tr_aux(:,:,:,m),CS%diag)
+ else
+ if (CS%id_tracer(m)>0) &
+ call post_data(CS%id_tracer(m),CS%tr(:,:,:,m),CS%diag)
+ endif
+ if (CS%id_tr_adx(m)>0) &
+ call post_data(CS%id_tr_adx(m),CS%tr_adx(m)%p(:,:,:),CS%diag)
+ if (CS%id_tr_ady(m)>0) &
+ call post_data(CS%id_tr_ady(m),CS%tr_ady(m)%p(:,:,:),CS%diag)
+ if (CS%id_tr_dfx(m)>0) &
+ call post_data(CS%id_tr_dfx(m),CS%tr_dfx(m)%p(:,:,:),CS%diag)
+ if (CS%id_tr_dfy(m)>0) &
+ call post_data(CS%id_tr_dfy(m),CS%tr_dfy(m)%p(:,:,:),CS%diag)
+ enddo
+
+end subroutine ISOMIP_tracer_column_physics
+
+!> This particular tracer package does not report anything back to the coupler.
+! The code that is here is just a rough guide for packages that would.
+subroutine ISOMIP_tracer_surface_state(state, h, G, CS)
+ type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
+ type(surface), intent(inout) :: state !< A structure containing fields that describe the surface state of the ocean.
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2.
+ type(ISOMIP_tracer_CS), pointer :: CS !< The control structure returned by a previous call to ISOMIP_register_tracer.
+ integer :: i, j, m, is, ie, js, je, nz
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+
+ if (.not.associated(CS)) return
+
+ if (CS%coupled_tracers) then
+ do m=1,ntr
+ ! This call loads the surface vlues into the appropriate array in the
+ ! coupler-type structure.
+ call set_coupler_values(CS%tr(:,:,1,1), state%tr_fields, CS%ind_tr(m), &
+ ind_csurf, is, ie, js, je)
+ enddo
+ endif
+
+end subroutine ISOMIP_tracer_surface_state
+
+subroutine ISOMIP_tracer_end(CS)
+ type(ISOMIP_tracer_CS), pointer :: CS
+ integer :: m
+
+ if (associated(CS)) then
+ if (associated(CS%tr)) deallocate(CS%tr)
+ if (associated(CS%tr_aux)) deallocate(CS%tr_aux)
+ do m=1,NTR
+ if (associated(CS%tr_adx(m)%p)) deallocate(CS%tr_adx(m)%p)
+ if (associated(CS%tr_ady(m)%p)) deallocate(CS%tr_ady(m)%p)
+ if (associated(CS%tr_dfx(m)%p)) deallocate(CS%tr_dfx(m)%p)
+ if (associated(CS%tr_dfy(m)%p)) deallocate(CS%tr_dfy(m)%p)
+ enddo
+
+ deallocate(CS)
+ endif
+end subroutine ISOMIP_tracer_end
+
+end module ISOMIP_tracer
diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90
index 49a40296c7..bc0127faf5 100644
--- a/src/tracer/MOM_tracer_flow_control.F90
+++ b/src/tracer/MOM_tracer_flow_control.F90
@@ -53,6 +53,9 @@ module MOM_tracer_flow_control
use DOME_tracer, only : register_DOME_tracer, initialize_DOME_tracer
use DOME_tracer, only : DOME_tracer_column_physics, DOME_tracer_surface_state
use DOME_tracer, only : DOME_tracer_end, DOME_tracer_CS
+use ISOMIP_tracer, only : register_ISOMIP_tracer, initialize_ISOMIP_tracer
+use ISOMIP_tracer, only : ISOMIP_tracer_column_physics, ISOMIP_tracer_surface_state
+use ISOMIP_tracer, only : ISOMIP_tracer_end, ISOMIP_tracer_CS
use ideal_age_example, only : register_ideal_age_tracer, initialize_ideal_age_tracer
use ideal_age_example, only : ideal_age_tracer_column_physics, ideal_age_tracer_surface_state
use ideal_age_example, only : ideal_age_stock, ideal_age_example_end, ideal_age_tracer_CS
@@ -84,6 +87,7 @@ module MOM_tracer_flow_control
type, public :: tracer_flow_control_CS ; private
logical :: use_USER_tracer_example = .false.
logical :: use_DOME_tracer = .false.
+ logical :: use_ISOMIP_tracer = .false.
logical :: use_ideal_age = .false.
logical :: use_regional_dyes = .false.
logical :: use_oil = .false.
@@ -92,6 +96,7 @@ module MOM_tracer_flow_control
logical :: use_MOM_generic_tracer = .false.
type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL()
type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL()
+ type(ISOMIP_tracer_CS), pointer :: ISOMIP_tracer_CSp => NULL()
type(ideal_age_tracer_CS), pointer :: ideal_age_tracer_CSp => NULL()
type(dye_tracer_CS), pointer :: dye_tracer_CSp => NULL()
type(oil_tracer_CS), pointer :: oil_tracer_CSp => NULL()
@@ -144,6 +149,9 @@ subroutine call_tracer_register(HI, GV, param_file, CS, tr_Reg, restart_CS)
call get_param(param_file, mod, "USE_DOME_TRACER", CS%use_DOME_tracer, &
"If true, use the DOME_tracer tracer package.", &
default=.false.)
+ call get_param(param_file, mod, "USE_ISOMIP_TRACER", CS%use_ISOMIP_tracer, &
+ "If true, use the ISOMIP_tracer tracer package.", &
+ default=.false.)
call get_param(param_file, mod, "USE_IDEAL_AGE_TRACER", CS%use_ideal_age, &
"If true, use the ideal_age_example tracer package.", &
default=.false.)
@@ -179,6 +187,9 @@ subroutine call_tracer_register(HI, GV, param_file, CS, tr_Reg, restart_CS)
if (CS%use_DOME_tracer) CS%use_DOME_tracer = &
register_DOME_tracer(HI, GV, param_file, CS%DOME_tracer_CSp, &
tr_Reg, restart_CS)
+ if (CS%use_ISOMIP_tracer) CS%use_ISOMIP_tracer = &
+ register_ISOMIP_tracer(HI, GV, param_file, CS%ISOMIP_tracer_CSp, &
+ tr_Reg, restart_CS)
if (CS%use_ideal_age) CS%use_ideal_age = &
register_ideal_age_tracer(HI, GV, param_file, CS%ideal_age_tracer_CSp, &
tr_Reg, restart_CS)
@@ -245,6 +256,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, diag, OB
if (CS%use_DOME_tracer) &
call initialize_DOME_tracer(restart, day, G, GV, h, diag, OBC, CS%DOME_tracer_CSp, &
sponge_CSp, diag_to_Z_CSp)
+ if (CS%use_ISOMIP_tracer) &
+ call initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS%ISOMIP_tracer_CSp, &
+ ALE_sponge_CSp, diag_to_Z_CSp)
if (CS%use_ideal_age) &
call initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS%ideal_age_tracer_CSp, &
sponge_CSp, diag_to_Z_CSp)
@@ -364,6 +378,9 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, o
if (CS%use_DOME_tracer) &
call DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
G, GV, CS%DOME_tracer_CSp)
+ if (CS%use_ISOMIP_tracer) &
+ call ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%ISOMIP_tracer_CSp)
if (CS%use_ideal_age) &
call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
G, GV, CS%ideal_age_tracer_CSp)
@@ -559,6 +576,8 @@ subroutine call_tracer_surface_state(state, h, G, CS)
call USER_tracer_surface_state(state, h, G, CS%USER_tracer_example_CSp)
if (CS%use_DOME_tracer) &
call DOME_tracer_surface_state(state, h, G, CS%DOME_tracer_CSp)
+ if (CS%use_ISOMIP_tracer) &
+ call ISOMIP_tracer_surface_state(state, h, G, CS%ISOMIP_tracer_CSp)
if (CS%use_ideal_age) &
call ideal_age_tracer_surface_state(state, h, G, CS%ideal_age_tracer_CSp)
if (CS%use_regional_dyes) &
@@ -582,6 +601,7 @@ subroutine tracer_flow_control_end(CS)
if (CS%use_USER_tracer_example) &
call USER_tracer_example_end(CS%USER_tracer_example_CSp)
if (CS%use_DOME_tracer) call DOME_tracer_end(CS%DOME_tracer_CSp)
+ if (CS%use_ISOMIP_tracer) call ISOMIP_tracer_end(CS%ISOMIP_tracer_CSp)
if (CS%use_ideal_age) call ideal_age_example_end(CS%ideal_age_tracer_CSp)
if (CS%use_regional_dyes) call regional_dyes_end(CS%dye_tracer_CSp)
if (CS%use_oil) call oil_tracer_end(CS%oil_tracer_CSp)