From ef17f74fc60c34c9aea9aa1a0a14b527da33733d Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Fri, 17 Jun 2016 17:04:59 +1000 Subject: [PATCH 1/5] Doxygen macro expansion and variable grouping We expand only the macro ALLOCABLE_ to the Fortran keyword 'allocatable'. This cleans up some type definitions, where doxygen thought that allocable_ was a variable name and got confused about what the type members actually were. We also turn on variable grouping with DISTRIBUTE_GROUP_DOC If there's a set of variables that all do the same thing, such as diagnostic identifiers, they can be grouped and the same docstring will be used for all members. For example: !>@{ !! Diagnostic identifiers integer :: id_u = -1, id_v = -1, id_h = -1 !>@} --- .doxygen | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.doxygen b/.doxygen index 43bb7d8eb6..81fa549885 100644 --- a/.doxygen +++ b/.doxygen @@ -341,7 +341,7 @@ 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 # 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 @@ -1940,7 +1940,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 +1948,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 +1980,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 From cb5ee49e87499931ad3e5eacf0e9152a2d246dbb Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Fri, 17 Jun 2016 17:20:34 +1000 Subject: [PATCH 2/5] Doxygenise MOM_vert_friction The existing comments have been tweaked so that they are Doxygen-aware for the subroutines and types defined in this module. I've also made an attempt to document the tridiagonal solver and the vertvisc routine --- .../vertical/MOM_vert_friction.F90 | 599 +++++++++--------- 1 file changed, 301 insertions(+), 298 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 0a2cbf0d17..7fdd3fb0a4 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 @@ -99,137 +28,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 @@ -311,6 +249,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)) @@ -325,6 +286,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 @@ -449,31 +413,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, @@ -565,37 +522,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). @@ -972,39 +921,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. @@ -1195,22 +1153,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. @@ -1393,34 +1355,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 @@ -1582,13 +1533,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 @@ -1620,6 +1573,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,3 +1584,52 @@ subroutine vertvisc_end(CS) end subroutine vertvisc_end end module MOM_vert_friction +!> \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: +!! +!! 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). From a20ae1154452b6618809ccbe971c0d2a5d73b2e6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 17 Jun 2016 10:52:13 -0400 Subject: [PATCH 3/5] Doxygen: Updated doxygen settings file to 1.8.12 - Regenerated the .doxygen file to have up to date documentation of settings, consistent with Doxygen version 1.8.12 . - I had to add *.F90 to the FILE_PATTERNS which seems to not be default anymore in newer doxygen versions. --- .doxygen | 106 ++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 89 insertions(+), 17 deletions(-) diff --git a/.doxygen b/.doxygen index 81fa549885..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 @@ -343,6 +352,13 @@ IDL_PROPERTY_SUPPORT = YES 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 # type (e.g. under the Public Functions section). Set it to NO to prevent @@ -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 #--------------------------------------------------------------------------- @@ -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. From 44a0091dc43893bed7e0fad5a33d75427fe6ab21 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 17 Jun 2016 10:57:11 -0400 Subject: [PATCH 4/5] Doxygen: Changed case of some \namespace arguments - Doxygen considers Fortran to be case insensitive and lowers the case of all objects. When documenting a module with the \namespace command we must use lower case for the module name even though the code uses capitalization. - Also added a \verbatim around a schematic of the staggered grid. --- src/ALE/MOM_regridding.F90 | 2 +- src/parameterizations/vertical/MOM_vert_friction.F90 | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) 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/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 7fdd3fb0a4..ffe469a694 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1583,8 +1583,7 @@ subroutine vertvisc_end(CS) deallocate(CS) end subroutine vertvisc_end -end module MOM_vert_friction -!> \namespace MOM_vert_friction +!> \namespace mom_vert_friction !! \author Robert Hallberg !! \date April 1994 - October 2006 !! @@ -1623,7 +1622,7 @@ end module MOM_vert_friction !! 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 @@ -1631,5 +1630,7 @@ end module MOM_vert_friction !! 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 From 04f7a9d29bbdcc1b1459305d9f86653db048fd62 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 17 Jun 2016 12:07:42 -0400 Subject: [PATCH 5/5] Added passive tracers to ISOMIP Created a module that sets up two passive tracers, one that is injected in the sponge layer and one that is injected based on the ice shelf melt water field. These tracers are activated by using USE_ISOMIP_TRACER = True. To add a tracer based on the ice shelf melt water I first included the ice shelf melt field (iceshelf_melt) in the structure "forcing", which contains pointers to the boundary forcing (defined in MOM_forcing_type.F90) so this field could be accessed from ISOMIP_tracer.F90. iceshelf_melt is passed to the forcing structure in MOM_ice_shelf.F90. Finally, the tracer is set in ISOMIP_tracer_column_physics (ISOMIP_tracer.F90). The first tracer injected in the sponge layer works fine, but I could not test the melt tracer yet since I first need to figure out how to properly activate the boundary layer scheme. The model compiles sucessfully. --- src/core/MOM_forcing_type.F90 | 26 +- src/ice_shelf/MOM_ice_shelf.F90 | 7 +- src/tracer/ISOMIP_tracer.F90 | 425 +++++++++++++++++++++++++ src/tracer/MOM_tracer_flow_control.F90 | 20 ++ 4 files changed, 469 insertions(+), 9 deletions(-) create mode 100644 src/tracer/ISOMIP_tracer.F90 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/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 3a71eff968..3b1591d031 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -52,6 +52,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 @@ -83,6 +86,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. @@ -91,6 +95,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() @@ -143,6 +148,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.) @@ -178,6 +186,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) @@ -244,6 +255,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) @@ -363,6 +377,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) @@ -558,6 +575,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) & @@ -581,6 +600,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)