From 0440ed45fd1460080910818685441b72def88a3b Mon Sep 17 00:00:00 2001 From: Yi-Cheng Teng - NOAA GFDL <143743249+yichengt900@users.noreply.github.com> Date: Fri, 19 Jul 2024 10:02:24 -0400 Subject: [PATCH 01/16] Perform unit conversion for internal heat only when it is applicable --- src/tracer/MOM_generic_tracer.F90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 5591e74d97..e998e3029c 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -584,13 +584,24 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga) else - call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + ! tv%internal_heat is a null pointer unless DO_GEOTHERMAL = True, + ! so we have to check and only do the scaling if it is associated. + if(associated(tv%internal_heat)) then + call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & optics%nbands, optics%max_wavelength_band, & sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & internal_heat=G%US%RZ_to_kg_m2*US%C_to_degC*tv%internal_heat(:,:), & frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga) + else + call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, & + sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & + opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga) + endif endif ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes From a8d43f718075bed7e082accb70d5de27042cf78f Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 23 Jul 2024 11:11:05 -0800 Subject: [PATCH 02/16] Adding Bodner reference --- docs/parameterizations_lateral.rst | 3 ++- docs/zotero.bib | 14 ++++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 3a3266a2bb..279a5749fe 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -32,7 +32,8 @@ Mixed layer restratification by sub-mesoscale eddies ---------------------------------------------------- Mixed layer restratification from :cite:`fox-kemper2008` and -:cite:`fox-kemper2008-2` is implemented in MOM_mixed_layer_restrat. +:cite:`fox-kemper2008-2` is implemented in MOM_mixed_layer_restrat, +which now also contains the mixed layer restratication comes from :cite: Bodner2023. Lateral diffusion ----------------- diff --git a/docs/zotero.bib b/docs/zotero.bib index c0f1ddccbb..997b1d24cb 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2760,3 +2760,17 @@ @article{Adcroft2019 journal = {J. Adv. Mod. Earth Sys.} } +@article{Bodner2023, + title={Modifying the Mixed Layer Eddy Parameterization to Include Frontogenesis Arrest by Boundary Layer Turbulence}, + volume={53}, + ISSN={1520-0485}, + url={http://dx.doi.org/10.1175/JPO-D-21-0297.1}, + DOI={10.1175/jpo-d-21-0297.1}, + number={1}, + journal={Journal of Physical Oceanography}, + publisher={American Meteorological Society}, + author={Bodner, Abigail S. and Fox-Kemper, Baylor and Johnson, Leah and Van Roekel, Luke P. and McWilliams, James C. and Sullivan, Peter P. and Hall, Paul S. and Dong, Jihai}, + year={2023}, + month=jan, + pages={323–339} +} From d76926cbd6e654530ff7f2bcd004da926c3f4a3f Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 23 Jul 2024 12:09:07 -0800 Subject: [PATCH 03/16] Failing attempt to link to mle docs. --- docs/parameterizations_lateral.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 279a5749fe..2b8341b881 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -35,6 +35,8 @@ Mixed layer restratification from :cite:`fox-kemper2008` and :cite:`fox-kemper2008-2` is implemented in MOM_mixed_layer_restrat, which now also contains the mixed layer restratication comes from :cite: Bodner2023. + :ref:`namespacemom__mixed__layer__restrat_1section_mle` + Lateral diffusion ----------------- From 780d9318577c40c6ca27b7e2d44138fe1dac6ce1 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 23 Jul 2024 13:18:59 -0800 Subject: [PATCH 04/16] Two more references --- docs/zotero.bib | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index 997b1d24cb..ff09ada402 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2772,5 +2772,32 @@ @article{Bodner2023 author={Bodner, Abigail S. and Fox-Kemper, Baylor and Johnson, Leah and Van Roekel, Luke P. and McWilliams, James C. and Sullivan, Peter P. and Hall, Paul S. and Dong, Jihai}, year={2023}, month=jan, - pages={323–339} + pages={323-–339} } + +@incollection{Niiler1977, + author={P. P. Niiler and E. B. Kraus}, + title={One-dimensional models of the upper ocean}, + booktitle={Modeling and Prediction of the Upper Layers of the Ocean}, + editor={E. B. Kraus}, + publisher={Pergamon}, + address={New York}, + pages={143-–172}, + year={1977} +} + +@article{Oberhuber1993, + title={Simulation of the Atlantic Circulation with a Coupled Sea Ice-Mixed Layer-Isopycnal General Circulation Model. Part I: Model Description}, + volume={23}, + ISSN={1520-0485}, + url={http://dx.doi.org/10.1175/1520-0485(1993)023<0808:SOTACW>2.0.CO;2}, + DOI={10.1175/1520-0485(1993)023<0808:sotacw>2.0.co;2}, + number={5}, + journal={Journal of Physical Oceanography}, + publisher={American Meteorological Society}, + author={Oberhuber, Josef M.}, + year={1993}, + month=may, + pages={808–829} +} + From ec40a7ca0b74d0e0f7c46302ce23adb7b0de8147 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 24 Jul 2024 10:25:28 -0800 Subject: [PATCH 05/16] Working on documentation links --- docs/parameterizations_lateral.rst | 29 +++- docs/zotero.bib | 159 +++++++++++++++++- .../lateral/MOM_hor_visc.F90 | 14 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 2 +- .../lateral/MOM_load_love_numbers.F90 | 10 +- .../lateral/MOM_mixed_layer_restrat.F90 | 22 +-- .../lateral/MOM_self_attr_load.F90 | 8 +- .../lateral/MOM_spherical_harmonics.F90 | 6 +- .../lateral/MOM_thickness_diffuse.F90 | 4 +- .../lateral/MOM_tidal_forcing.F90 | 2 +- 10 files changed, 222 insertions(+), 34 deletions(-) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 2b8341b881..3a444098b2 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -8,6 +8,8 @@ Lateral viscosity Laplacian and bi-harmonic viscosities with linear and Smagorinsky options are implemented in MOM_hor_visc. + :ref:`namespacemom__hor__visc_1section_horizontal_viscosity` + Gent-McWilliams/TEM/isopycnal height diffusion ---------------------------------------------- @@ -20,7 +22,7 @@ scaling. A model of sub-grid scale Mesoscale Eddy Kinetic Energy (MEKE) is implement in MOM_MEKE and the associated diffusivity added in MOM_thickness_diffuse. See :cite:`jansen2015` and :cite:`marshall2010`. - :ref:`namespacemom__meke_1section_MEKE` + :ref:`namespacemom__meke_1section_MEKE` Backscatter ----------- @@ -37,15 +39,38 @@ which now also contains the mixed layer restratication comes from :cite: Bodner2 :ref:`namespacemom__mixed__layer__restrat_1section_mle` +Interface filtering +------------------- + +For layer mode, one can filter the interface thicknesses: + + :ref:`namespacemom__interface__filter_1section_interface_filter` + Lateral diffusion ----------------- See :ref:`Horizontal_Diffusion`. +See also :ref:`namespacemom__lateral__mixing__coeffs_1section_Resolution_Function` + Tidal forcing ------------- Astronomical tidal forcings and self-attraction and loading are implement in MOM_tidal_forcing. Tides can also be added via an open boundary tidal specification, -see [OBC wiki page](https://github.com/NOAA-GFDL/MOM6-examples/wiki/Open-Boundary-Conditions). +see `OBC wiki page `_. + +The Love numbers are stored internally in MOM_load_love_numbers: + + :ref:`namespacemom__load__love__numbers_1section_Love_numbers` + +While the self attraction and loading is computed in MOM_self_attr_load: + + :ref:`namespaceself__attr__load_1section_SAL` + +Spherical Harmonics +------------------- + +The model needs spherical, computed in MOM_spherical_harmonics: + :ref:`namespacemom__spherical__harmonics_1section_spherical_harmonics` diff --git a/docs/zotero.bib b/docs/zotero.bib index ff09ada402..ededf3ac55 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2524,7 +2524,7 @@ @article{visbeck1997 } @article{visbeck1996, - author = {Viscbeck, M. and J.C. Marshall and H. Jones}, + author = {Visbeck, M. and J.C. Marshall and H. Jones}, year = {1996}, title = {Dynamics of isolated convective regions in the ocean}, journal = {J. Phys. Oceanogr.}, @@ -2801,3 +2801,160 @@ @article{Oberhuber1993 pages={808–829} } +@article{Smith2003, + title={Anisotropic horizontal viscosity for ocean models}, + volume={5}, + ISSN={1463-5003}, + url={http://dx.doi.org/10.1016/s1463-5003(02)00016-1}, + DOI={10.1016/s1463-5003(02)00016-1}, + number={2}, + journal={Ocean Modelling}, + publisher={Elsevier BV}, + author={Smith, Richard D. and McWilliams, James C.}, + year={2003}, + month=jan, + pages={129–156} +} + +@article{Large2001, + title={Equatorial Circulation of a Global Ocean Climate Model with Anisotropic Horizontal Viscosity}, + volume={31}, + ISSN={1520-0485}, + url={http://dx.doi.org/10.1175/1520-0485(2001)031<0518:ECOAGO>2.0.CO;2}, + DOI={10.1175/1520-0485(2001)031<0518:ecoago>2.0.co;2}, + number={2}, + journal={Journal of Physical Oceanography}, + publisher={American Meteorological Society}, + author={Large, William G. and Danabasoglu, Gokhan and McWilliams, James C. and Gent, Peter R. and Bryan, Frank O.}, + year={2001}, + month=feb, + pages={518–536} +} + +@inproceedings{Smagorinsky1993, + author={Joseph Smagorinsky}, + year={1993}, + title={Some historical remarks on the use of non-linear viscosities}, + booktitle={Large Eddy Simulation of Complex Engineering and Geophysical Flows}, + note={Proceedings of an International Workshop in Large Eddy Simulation}, + address={Cambridge, UK}, + publisher={Cambridge University Press}, + pages={1--34} +} + +@article{Barton2022, + title={Global Barotropic Tide Modeling Using Inline Self‐Attraction and Loading in MPAS‐Ocean}, + volume={14}, + ISSN={1942-2466}, + url={http://dx.doi.org/10.1029/2022MS003207}, + DOI={10.1029/2022ms003207}, + number={11}, + journal={Journal of Advances in Modeling Earth Systems}, + publisher={American Geophysical Union (AGU)}, + author={Barton, Kristin N. and Pal, Nairita and Brus, Steven R. and Petersen, Mark R. and Arbic, Brian K. and Engwirda, Darren and Roberts, Andrew F. and Westerink, Joannes J. and Wirasaet, Damrongsak and Schindelegger, Michael}, + year={2022}, + month=nov +} + +@article{Brus2023, + title={Scalable self attraction and loading calculations for unstructured ocean tide models}, + volume={182}, + ISSN={1463-5003}, + url={http://dx.doi.org/10.1016/j.ocemod.2023.102160}, + DOI={10.1016/j.ocemod.2023.102160}, + journal={Ocean Modelling}, + publisher={Elsevier BV}, + author={Brus, Steven R. and Barton, Kristin N. and Pal, Nairita and Roberts, Andrew F. and Engwirda, Darren and Petersen, Mark R. and Arbic, Brian K. and Wirasaet, Damrongsak and Westerink, Joannes J. and Schindelegger, Michael}, + year={2023}, + month=apr, + pages={102160} +} + +@article{Blewitt2003, + title={Self‐consistency in reference frames, geocenter definition, and surface loading of the solid Earth}, + volume={108}, + ISSN={0148-0227}, + url={http://dx.doi.org/10.1029/2002JB002082}, + DOI={10.1029/2002jb002082}, + number={B2}, + journal={Journal of Geophysical Research: Solid Earth}, + publisher={American Geophysical Union (AGU)}, + author={Blewitt, Geoffrey}, + year={2003}, + month=feb +} + +@article{Wang2012, + author={Wang, H., Xiang, L., Jia, L., Jiang, L., Wang, Z., Hu, B. and Gao, P.}, + year={2012}, + title={Load Love numbers and Green's functions +for elastic Earth models PREM, iasp91, ak135, and modified models with refined crustal structure from Crust 2.0}, + journal={Computers & Geosciences}, + volume={49}, + pages={190--199} +} + +@article{Jansen2015, + title={Parameterization of eddy fluxes based on a mesoscale energy budget}, + volume={92}, + ISSN={1463-5003}, + url={http://dx.doi.org/10.1016/j.ocemod.2015.05.007}, + DOI={10.1016/j.ocemod.2015.05.007}, + journal={Ocean Modelling}, + publisher={Elsevier BV}, + author={Jansen, Malte F. and Adcroft, Alistair J. and Hallberg, Robert and Held, Isaac M.}, + year={2015}, + month=aug, + pages={28–-41} +} + +@inproceedings{Hallberg2003, + title={The ability of large-scale ocean models to accept parameterizations of boundary mixing, and a description of a refined bulk mixed-layer model}, + author={Robert Hallberg}, + year={2003}, + booktitle={Internal Gravity Waves and Small-Scale Turbulence: Proc.‘Aha Huliko ‘a Hawaiian Winter Workshop}, + pages={187--203} +} + +@article{1978, + volume={290}, + ISSN={2054-0272}, + url={http://dx.doi.org/10.1098/rsta.1978.0083}, + DOI={10.1098/rsta.1978.0083}, + number={1368}, + journal={Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences}, + publisher={The Royal Society}, + year={1978}, + month=nov, + pages={235-–266} +} + +@article{Arbic2004, + title={The accuracy of surface elevations in forward global barotropic and baroclinic tide models}, + volume={51}, + ISSN={0967-0645}, + url={http://dx.doi.org/10.1016/j.dsr2.2004.09.014}, + DOI={10.1016/j.dsr2.2004.09.014}, + number={25–26}, + journal={Deep Sea Research Part II: Topical Studies in Oceanography}, + publisher={Elsevier BV}, + author={Arbic, Brian K. and Garner, Stephen T. and Hallberg, Robert W. and Simmons, Harper L.}, + year={2004}, + month=dec, + pages={3069-–3101} +} + +@article{Schaeffer2013, + title={Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations}, + volume={14}, + ISSN={1525-2027}, + url={http://dx.doi.org/10.1002/ggge.20071}, + DOI={10.1002/ggge.20071}, + number={3}, + journal={Geochemistry, Geophysics, Geosystems}, + publisher={American Geophysical Union (AGU)}, + author={Schaeffer, Nathanaël}, + year={2013}, + month=mar, + pages={751-–758} +} diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 2a4181804b..ea06795654 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -3193,25 +3193,25 @@ subroutine hor_visc_end(CS) end subroutine hor_visc_end !> \namespace mom_hor_visc !! +!! \section section_horizontal_viscosity Horizontal viscosity in MOM +!! !! This module contains the subroutine horizontal_viscosity() that calculates the !! effects of horizontal viscosity, including parameterizations of the value of !! the viscosity itself. horizontal_viscosity() calculates the acceleration due to !! some combination of a biharmonic viscosity and a Laplacian viscosity. Either or !! both may use a coefficient that depends on the shear and strain of the flow. !! All metric terms are retained. The Laplacian is calculated as the divergence of -!! a stress tensor, using the form suggested by Smagorinsky (1993). The biharmonic +!! a stress tensor, using the form suggested by \cite Smagorinsky1993. The biharmonic !! is calculated by twice applying the divergence of the stress tensor that is !! used to calculate the Laplacian, but without the dependence on thickness in the !! first pass. This form permits a variable viscosity, and indicates no !! acceleration for either resting fluid or solid body rotation. !! -!! The form of the viscous accelerations is discussed extensively in Griffies and -!! Hallberg (2000), and the implementation here follows that discussion closely. -!! We use the notation of Smith and McWilliams (2003) with the exception that the +!! The form of the viscous accelerations is discussed extensively in \cite griffies2000, +!! and the implementation here follows that discussion closely. +!! We use the notation of \cite Smith2003 with the exception that the !! isotropic viscosity is \f$\kappa_h\f$. !! -!! \section section_horizontal_viscosity Horizontal viscosity in MOM -!! !! In general, the horizontal stress tensor can be written as !! \f[ !! {\bf \sigma} = @@ -3357,7 +3357,7 @@ end subroutine hor_visc_end !! !! \subsection section_anisotropic_viscosity Anisotropic viscosity !! -!! Large et al., 2001, proposed enhancing viscosity in a particular direction and the +!! \cite Large2001, proposed enhancing viscosity in a particular direction and the !! approach was generalized in Smith and McWilliams, 2003. We use the second form of their !! two coefficient anisotropic viscosity (section 4.3). We also replace their !! \f$A^\prime\f$ and $D$ such that \f$2A^\prime = 2 \kappa_h + D\f$ and diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 7315c82601..1a9afad897 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1728,7 +1728,7 @@ end subroutine VarMix_end !! \f] !! !! \todo Check this reference to Bob on/off paper. -!! The resolution function used in scaling diffusivities (Hallberg, 2010) is +!! The resolution function used in scaling diffusivities (\cite hallberg2013) is !! !! \f[ !! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p} diff --git a/src/parameterizations/lateral/MOM_load_love_numbers.F90 b/src/parameterizations/lateral/MOM_load_love_numbers.F90 index 3d573d894d..94a107c4c7 100644 --- a/src/parameterizations/lateral/MOM_load_love_numbers.F90 +++ b/src/parameterizations/lateral/MOM_load_love_numbers.F90 @@ -1452,15 +1452,17 @@ module MOM_load_love_numbers /), (/4, lmax+1/)) !< Load Love numbers !> \namespace mom_load_love_numbers +!! \section section_Love_numbers The Love numbers +!! !! This module serves the sole purpose of storing load Love number. The Love numbers are used for the spherical harmonic !! self-attraction and loading (SAL) calculation in MOM_self_attr_load module. This separate module ensures readability !! of the SAL module. !! !! Variable Love_Data stores the Love numbers up to degree 1440. From left to right: degree, h, l, and k. Data in this !! module is imported from SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean developed by Los Alamos -!! National Laboratory and University of Michigan [Barton et al. (2022) and Brus et al. (2022)]. The load Love numbers -!! are from Wang et al. (2012), which are in the center of mass of total Earth system reference frame (CM). When used, -!! Love numbers with degree<2 should be converted to center of mass solid Earth reference frame (CF) [Blewitt (2003)], +!! National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2022]. The load Love numbers +!! are from \cite Wang2012, which are in the center of mass of total Earth system reference frame (CM). When used, +!! Love numbers with degree<2 should be converted to center of mass solid Earth reference frame (CF) [\cite Blewitt2003], !! as in subroutine calc_love_scaling in MOM_tidal_forcing module. !! !! References: @@ -1483,4 +1485,4 @@ module MOM_load_love_numbers !! for elastic Earth models PREM, iasp91, ak135, and modified models with refined crustal structure from Crust 2.0. !! Computers & Geosciences, 49, pp.190-199. !! https://doi.org/10.1016/j.cageo.2012.06.022 -end module MOM_load_love_numbers \ No newline at end of file +end module MOM_load_love_numbers diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index d1e0d15293..18c219bc7f 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -1704,7 +1704,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, "mesoscale eddy kinetic energy to the large-scale "//& "geostrophic kinetic energy or 1 plus the square of the "//& "grid spacing over the deformation radius, as detailed "//& - "by Fox-Kemper et al. (2010)", units="nondim", default=0.0) + "by Fox-Kemper et al. (2011)", units="nondim", default=0.0) ! These parameters are only used in the OM4-era version of Fox-Kemper call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_Stanley_ML, & "If true, turn on Stanley SGS T variance parameterization "// & @@ -1750,7 +1750,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, call get_param(param_file, mdl, "MLE_DENSITY_DIFF", CS%MLE_density_diff, & "Density difference used to detect the mixed-layer "//& "depth used for the mixed-layer eddy parameterization "//& - "by Fox-Kemper et al. (2010)", units="kg/m3", default=0.03, scale=US%kg_m3_to_R) + "by Fox-Kemper et al. (2011)", units="kg/m3", default=0.03, scale=US%kg_m3_to_R) endif call get_param(param_file, mdl, "MLE_TAIL_DH", CS%MLE_tail_dh, & "Fraction by which to extend the mixed-layer restratification "//& @@ -1965,14 +1965,14 @@ end function test_answer !! \section section_mle Mixed-layer eddy parameterization module !! !! The subroutines in this module implement a parameterization of unresolved viscous -!! mixed layer restratification of the mixed layer as described in Fox-Kemper et -!! al., 2008, and whose impacts are described in Fox-Kemper et al., 2011. +!! mixed layer restratification of the mixed layer as described in \cite fox-kemper2008, +!! and whose impacts are described in \cite fox-kemper2011. !! This is derived in part from the older parameterization that is described in -!! Hallberg (Aha Hulikoa, 2003), which this new parameterization surpasses, which -!! in turn is based on the sub-inertial mixed layer theory of Young (JPO, 1994). +!! \cite Hallberg2003, which this new parameterization surpasses, which +!! in turn is based on the sub-inertial mixed layer theory of \cite Young1994. !! There is no net horizontal volume transport due to this parameterization, and !! no direct effect below the mixed layer. A revised of the parameterization by -!! Bodner et al., 2023, is also available as an option. +!! \cite Bodner2023, is also available as an option. !! !! This parameterization sets the restratification timescale to agree with !! high-resolution studies of mixed layer restratification. @@ -1986,8 +1986,8 @@ end function test_answer !! !! The parameterization is colloquially referred to as "sub-meso". !! -!! The original Fox-Kemper et al., (2008b) paper proposed a quasi-Stokes -!! advection described by the stream function (eq. 5 of Fox-Kemper et al., 2011): +!! The original \cite fox-kemper2008-2 paper proposed a quasi-Stokes +!! advection described by the stream function (eq. 5 of \cite fox-kemper2011): !! \f[ !! {\bf \Psi}_o = C_e \frac{ H^2 \nabla \bar{b} \times \hat{\bf z} }{ |f| } \mu(z) !! \f] @@ -2001,7 +2001,7 @@ end function test_answer !! \f$ \nabla \bar{b} \f$ is a depth mean buoyancy gradient averaged over the mixed layer. !! !! For use in coarse-resolution models, an upscaling of the buoyancy gradients and adaption for the equator -!! leads to the following parameterization (eq. 6 of Fox-Kemper et al., 2011): +!! leads to the following parameterization (eq. 6 of \cite fox-kemper2011): !! \f[ !! {\bf \Psi} = C_e \Gamma_\Delta \frac{\Delta s}{l_f} \frac{ H^2 \nabla \bar{b} \times \hat{\bf z} } !! { \sqrt{ f^2 + \tau^{-2}} } \mu(z) @@ -2057,7 +2057,7 @@ end function test_answer !! available parameters. !! MLE_USE_PBL_MLD must be True to use the B23 modification. !! -!! Bodner et al., 2023, (B23) use an expression for the frontal width which changes the scaling from \f$ H^2 \f$ +!! \cite Bodner2023, (B23) use an expression for the frontal width which changes the scaling from \f$ H^2 \f$ !! to \f$ h H^2 \f$: !! \f[ !! {\bf \Psi} = C_r \frac{\Delta s |f| \bar{h} \bar{H}^2 \nabla \bar{b} \times \hat{\bf z} } diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index e7f8a73ab2..df22c9494c 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -245,19 +245,21 @@ end subroutine SAL_end !> \namespace self_attr_load !! +!! \section section_SAL Self attraction and loading +!! !! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height (SSH) !! (rather, it should be bottom pressure anomaly). SAL is primarily used for fast evolving processes like tides or !! storm surges, but the effect applies to all motions. !! -!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied (Accad and Pekeris 1978) and the SAL is simply +!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied (\cite Accad1978) and the SAL is simply !! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH . For tides, the scalar !! approximation can also be used to iterate the SAL to convergence [see USE_PREVIOUS_TIDES in MOM_tidal_forcing, -!! Arbic et al. (2004)]. +!! \cite Arbic2004]. !! !! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate SAL. !! Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is set by !! SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean -!! developed by Los Alamos National Laboratory and University of Michigan [Barton et al. (2022) and Brus et al. (2023)]. +!! developed by Los Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. !! !! References: !! diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 index 26258e6b8e..1782b5ae67 100644 --- a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 +++ b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 @@ -334,6 +334,8 @@ end function order2index !> \namespace mom_spherical_harmonics !! +!! \section section_spherical_harmonics Spherical harmonics +!! !! This module contains the subroutines to calculate spherical harmonic transforms (SHT), namely, forward transform !! of a two-dimensional field into a given number of spherical harmonic modes and its inverse transform. This module !! is primarily used to but not limited to calculate self-attraction and loading (SAL) term, which is mostly relevant to @@ -341,8 +343,8 @@ end function order2index !! Currently, the transforms are for t-cell fields only. !! !! This module is stemmed from SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean developed by Los -!! Alamos National Laboratory and University of Michigan [Barton et al. (2022) and Brus et al. (2023)]. The algorithm -!! for forward and inverse transforms loosely follows Schaeffer (2013). +!! Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. The algorithm +!! for forward and inverse transforms loosely follows \cite Schaeffer2013. !! !! In forward transform, a two-dimensional physical field can be projected into a series of spherical harmonics. The !! spherical harmonic coefficient of degree n and order m for a field \f$f(\theta, \phi)\f$ is calculated as follows: diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index d735248417..d884a61aee 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -2460,7 +2460,7 @@ end subroutine thickness_diffuse_end !! since the quantity \f$\frac{M^2}{\sqrt{N^4 + M^4}}\f$ is bounded between $-1$ and $1$ and does not change sign !! if \f$N^2<0\f$. !! -!! Optionally, the method of Ferrari et al, 2010, can be used to obtain the streamfunction which solves the +!! Optionally, the method of \cite ferrari2010, can be used to obtain the streamfunction which solves the !! vertically elliptic equation: !! \f[ !! \gamma_F \partial_z c^2 \partial_z \psi - N_*^2 \psi = -( 1 + \gamma_F ) \kappa_h N_*^2 \frac{M^2}{\sqrt{N^4+M^4}} @@ -2478,7 +2478,7 @@ end subroutine thickness_diffuse_end !! \f$ r(\Delta x,L_d) \f$ is a function of the local resolution (ratio of grid-spacing, \f$\Delta x\f$, !! to deformation radius, \f$L_d\f$). The length \f$L_s\f$ is provided by the mom_lateral_mixing_coeffs module !! (enabled with USE_VARIABLE_MIXING=True and the term \f$\f$ is the vertical average slope -!! times the buoyancy frequency prescribed by Visbeck et al., 1996. +!! times the buoyancy frequency prescribed by \cite visbeck1996. !! !! The result of the above expression is subsequently bounded by minimum and maximum values, including an upper !! diffusivity consistent with numerical stability (\f$ \kappa_{cfl} \f$ is calculated internally). diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index d53c3321f9..740977c7b9 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -765,7 +765,7 @@ end subroutine tidal_forcing_end !! are provided. TIDAL_SAL_FROM_FILE can be set to read the phase and !! amplitude of the tidal SAL. USE_PREVIOUS_TIDES may be useful in !! combination with the scalar approximation to iterate the SAL to -!! convergence (for details, see Arbic et al., 2004, DSR II). With +!! convergence (for details, see \cite Arbic2004). With !! TIDAL_SAL_FROM_FILE or USE_PREVIOUS_TIDES, a list of input files !! must be provided to describe each constituent's properties from !! a previous solution. The online SAL calculations that are functions From 0e17b7389fb3073487bd0317c7fa1b3b28fb0447 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 24 Jul 2024 11:03:53 -0800 Subject: [PATCH 06/16] Added Young reference --- docs/zotero.bib | 9 +++++++++ src/parameterizations/lateral/MOM_load_love_numbers.F90 | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index ededf3ac55..8ea5c05c38 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2958,3 +2958,12 @@ @article{Schaeffer2013 month=mar, pages={751-–758} } + +@article{Young1994, + author={Young, W.}, + title={The subinertial mixed layer approximation}, + journal={J. Phys. Oceanogr.}, + volume={24}, + pages={1812--1826}, + year={1994} +} diff --git a/src/parameterizations/lateral/MOM_load_love_numbers.F90 b/src/parameterizations/lateral/MOM_load_love_numbers.F90 index 94a107c4c7..b0251f5fb0 100644 --- a/src/parameterizations/lateral/MOM_load_love_numbers.F90 +++ b/src/parameterizations/lateral/MOM_load_love_numbers.F90 @@ -1460,7 +1460,7 @@ module MOM_load_love_numbers !! !! Variable Love_Data stores the Love numbers up to degree 1440. From left to right: degree, h, l, and k. Data in this !! module is imported from SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean developed by Los Alamos -!! National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2022]. The load Love numbers +!! National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. The load Love numbers !! are from \cite Wang2012, which are in the center of mass of total Earth system reference frame (CM). When used, !! Love numbers with degree<2 should be converted to center of mass solid Earth reference frame (CF) [\cite Blewitt2003], !! as in subroutine calc_love_scaling in MOM_tidal_forcing module. From 369d716ba8da35ef66282cdf4133624271c82355 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 24 Jul 2024 11:20:11 -0800 Subject: [PATCH 07/16] Fix line length in comment --- src/parameterizations/lateral/MOM_load_love_numbers.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_load_love_numbers.F90 b/src/parameterizations/lateral/MOM_load_love_numbers.F90 index b0251f5fb0..c8bd8bbc3c 100644 --- a/src/parameterizations/lateral/MOM_load_love_numbers.F90 +++ b/src/parameterizations/lateral/MOM_load_love_numbers.F90 @@ -1462,8 +1462,8 @@ module MOM_load_love_numbers !! module is imported from SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean developed by Los Alamos !! National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. The load Love numbers !! are from \cite Wang2012, which are in the center of mass of total Earth system reference frame (CM). When used, -!! Love numbers with degree<2 should be converted to center of mass solid Earth reference frame (CF) [\cite Blewitt2003], -!! as in subroutine calc_love_scaling in MOM_tidal_forcing module. +!! Love numbers with degree<2 should be converted to center of mass solid Earth reference frame (CF) +!! [\cite Blewitt2003], as in subroutine calc_love_scaling in MOM_tidal_forcing module. !! !! References: !! From 86b4c85e5ecd02cc2ec5361f84a65bbed6bda3cc Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 24 Jul 2024 11:40:05 -0800 Subject: [PATCH 08/16] Delete duplicate entry --- docs/forcing.rst | 12 ++++++------ docs/zotero.bib | 11 ----------- .../vertical/MOM_bulk_mixed_layer.F90 | 2 +- 3 files changed, 7 insertions(+), 18 deletions(-) diff --git a/docs/forcing.rst b/docs/forcing.rst index 8496317608..21539d46d0 100644 --- a/docs/forcing.rst +++ b/docs/forcing.rst @@ -1,23 +1,23 @@ Forcing ======= Data Override -------- +------------- When running MOM6 with the Flexible Modelling System (FMS) coupler, forcing can be specified by a `data_table` file. This is particularly useful when running MOM6 with a data atmosphere, as paths to the relevent atmospheric forcing products (eg. JRA55-do or ERA5) can be provided here. Each item in the data table must be separated by a new line, and contains the following information: | ``gridname``: The component of the model this data applies to. eg. `atm` `ocn` `lnd` `ice`. | ``fieldname_code``: The field name according to the model component. eg. `salt` -| ``fieldname_file``: The name of the field within the source file. +| ``fieldname_file``: The name of the field within the source file. | ``file_name``: Path to the source file. | ``interpol_method``: Interpolation method eg. `bilinear` -| ``factor``: A scalar by which to multiply the field ahead of passing it onto the model. This is a quick way to do unit conversions for example. +| ``factor``: A scalar by which to multiply the field ahead of passing it onto the model. This is a quick way to do unit conversions for example. +| -| The data table is commonly formatted by specifying each of the fields in the order listed above, with a new line for each entry. Example Format: "ATM", "t_bot", "t2m", "./INPUT/2t_ERA5.nc", "bilinear", 1.0 -A `yaml` format is also possible if you prefer. This is outlined in the `FMS data override `_ github page, along with other details. +A `yaml` format is also possible if you prefer. This is outlined in the `FMS data override `_ github page, along with other details. Speficying a constant value: Rather than overriding with data from a file, one can also set a field to constant. To do this, pass empty strings to `fieldname_file` and `file_name`. The `factor` now corresponds to the override value. For example, the following sets the temperature at the bottom of the atmosphere to 290 Kelvin. @@ -26,7 +26,7 @@ Speficying a constant value: "ATM", "t_bot", "", "", "bilinear", 290.0 Which units do I need? - For configurations using SIS2 and MOM, a list of available surface flux variables along with the expected units can be found in the `flux_exchange `_ file. + For configurations using SIS2 and MOM, a list of available surface flux variables along with the expected units can be found in the `flux_exchange `_ file. .. toctree:: diff --git a/docs/zotero.bib b/docs/zotero.bib index 8ea5c05c38..c58089236b 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2775,17 +2775,6 @@ @article{Bodner2023 pages={323-–339} } -@incollection{Niiler1977, - author={P. P. Niiler and E. B. Kraus}, - title={One-dimensional models of the upper ocean}, - booktitle={Modeling and Prediction of the Upper Layers of the Ocean}, - editor={E. B. Kraus}, - publisher={Pergamon}, - address={New York}, - pages={143-–172}, - year={1977} -} - @article{Oberhuber1993, title={Simulation of the Atlantic Circulation with a Coupled Sea Ice-Mixed Layer-Isopycnal General Circulation Model. Part I: Model Description}, volume={23}, diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 561ace60a7..4c39d9c11e 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -4256,7 +4256,7 @@ end function EF4 !! !! This file contains the subroutine (bulkmixedlayer) that !! implements a Kraus-Turner-like bulk mixed layer, based on the work -!! of various people, as described in the review paper by \cite Niiler1977, +!! of various people, as described in the review paper by \cite niiler1977, !! with particular attention to the form proposed by \cite Oberhuber1993, !! with an extension to a refined bulk mixed layer as described in !! Hallberg (\cite muller2003). The physical processes portrayed in From fb40a4ace2b8ecf5c72f3b994e43290b1b6042d2 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 24 Jul 2024 12:01:46 -0800 Subject: [PATCH 09/16] Another duplicate bib entry --- docs/zotero.bib | 2 +- src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index c58089236b..fe5dd37909 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2775,7 +2775,7 @@ @article{Bodner2023 pages={323-–339} } -@article{Oberhuber1993, +@article{Oberhuber1993a, title={Simulation of the Atlantic Circulation with a Coupled Sea Ice-Mixed Layer-Isopycnal General Circulation Model. Part I: Model Description}, volume={23}, ISSN={1520-0485}, diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 4c39d9c11e..792dae1411 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -4257,7 +4257,7 @@ end function EF4 !! This file contains the subroutine (bulkmixedlayer) that !! implements a Kraus-Turner-like bulk mixed layer, based on the work !! of various people, as described in the review paper by \cite niiler1977, -!! with particular attention to the form proposed by \cite Oberhuber1993, +!! with particular attention to the form proposed by \cite Oberhuber1993a, !! with an extension to a refined bulk mixed layer as described in !! Hallberg (\cite muller2003). The physical processes portrayed in !! this subroutine include convective adjustment and mixed layer entrainment From 4fa86fdbcc042e2ca34076df276cd745c20a4773 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 24 Jul 2024 21:32:44 -0800 Subject: [PATCH 10/16] Another duplicate ref --- docs/zotero.bib | 19 +++---------------- src/parameterizations/lateral/MOM_MEKE.F90 | 6 +++--- .../lateral/MOM_hor_visc.F90 | 4 ++-- .../lateral/MOM_load_love_numbers.F90 | 2 +- 4 files changed, 9 insertions(+), 22 deletions(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index fe5dd37909..abaeea4c72 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2873,8 +2873,9 @@ @article{Blewitt2003 month=feb } -@article{Wang2012, - author={Wang, H., Xiang, L., Jia, L., Jiang, L., Wang, Z., Hu, B. and Gao, P.}, +@article{Wang2012-2, + author={Wang, H. and Xiang, L. and Jia, L. and Jiang, L. and Wang, Z. and Hu, B. + and Gao, P.}, year={2012}, title={Load Love numbers and Green's functions for elastic Earth models PREM, iasp91, ak135, and modified models with refined crustal structure from Crust 2.0}, @@ -2883,20 +2884,6 @@ @article{Wang2012 pages={190--199} } -@article{Jansen2015, - title={Parameterization of eddy fluxes based on a mesoscale energy budget}, - volume={92}, - ISSN={1463-5003}, - url={http://dx.doi.org/10.1016/j.ocemod.2015.05.007}, - DOI={10.1016/j.ocemod.2015.05.007}, - journal={Ocean Modelling}, - publisher={Elsevier BV}, - author={Jansen, Malte F. and Adcroft, Alistair J. and Hallberg, Robert and Held, Isaac M.}, - year={2015}, - month=aug, - pages={28–-41} -} - @inproceedings{Hallberg2003, title={The ability of large-scale ocean models to accept parameterizations of boundary mixing, and a description of a refined bulk mixed-layer model}, author={Robert Hallberg}, diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index f3bd9e8934..c3859c3bd1 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -1892,7 +1892,7 @@ end subroutine MEKE_end !! \f$ U_b \f$ is a constant background bottom velocity scale and is !! typically not used (i.e. set to zero). !! -!! Following Jansen et al., 2015, the projection of eddy energy on to the bottom +!! Following \cite jansen2015, the projection of eddy energy on to the bottom !! is given by the ratio of bottom energy to column mean energy: !! \f[ !! \gamma_b^2 = \frac{E_b}{E} = \gamma_{d0} @@ -1924,12 +1924,12 @@ end subroutine MEKE_end !! \f[ \kappa_M = \gamma_\kappa \sqrt{ \gamma_t^2 U_e^2 A_\Delta } \f] !! !! where \f$ A_\Delta \f$ is the area of the grid cell. -!! Following Jansen et al., 2015, we now use +!! Following \cite jansen2015, we now use !! !! \f[ \kappa_M = \gamma_\kappa l_M \sqrt{ \gamma_t^2 U_e^2 } \f] !! !! where \f$ \gamma_\kappa \in [0,1] \f$ is a non-dimensional factor and, -!! following Jansen et al., 2015, \f$\gamma_t^2\f$ is the ratio of barotropic +!! following \cite jansen2015, \f$\gamma_t^2\f$ is the ratio of barotropic !! eddy energy to column mean eddy energy given by !! \f[ !! \gamma_t^2 = \frac{E_t}{E} = \left( 1 + c_{t} \frac{L_d}{L_f} \right)^{-\frac{1}{4}} diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index ea06795654..4a901beab0 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -3195,9 +3195,9 @@ end subroutine hor_visc_end !! !! \section section_horizontal_viscosity Horizontal viscosity in MOM !! -!! This module contains the subroutine horizontal_viscosity() that calculates the +!! This module contains the subroutine horizontal_viscosity that calculates the !! effects of horizontal viscosity, including parameterizations of the value of -!! the viscosity itself. horizontal_viscosity() calculates the acceleration due to +!! the viscosity itself. Subroutine horizontal_viscosity calculates the acceleration due to !! some combination of a biharmonic viscosity and a Laplacian viscosity. Either or !! both may use a coefficient that depends on the shear and strain of the flow. !! All metric terms are retained. The Laplacian is calculated as the divergence of diff --git a/src/parameterizations/lateral/MOM_load_love_numbers.F90 b/src/parameterizations/lateral/MOM_load_love_numbers.F90 index c8bd8bbc3c..8faf3aafab 100644 --- a/src/parameterizations/lateral/MOM_load_love_numbers.F90 +++ b/src/parameterizations/lateral/MOM_load_love_numbers.F90 @@ -1461,7 +1461,7 @@ module MOM_load_love_numbers !! Variable Love_Data stores the Love numbers up to degree 1440. From left to right: degree, h, l, and k. Data in this !! module is imported from SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean developed by Los Alamos !! National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. The load Love numbers -!! are from \cite Wang2012, which are in the center of mass of total Earth system reference frame (CM). When used, +!! are from \cite Wang2012-2, which are in the center of mass of total Earth system reference frame (CM). When used, !! Love numbers with degree<2 should be converted to center of mass solid Earth reference frame (CF) !! [\cite Blewitt2003], as in subroutine calc_love_scaling in MOM_tidal_forcing module. !! From d93f8bf6b2b8bdb5d661553fbb82dca2eca21c33 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 25 Jul 2024 14:13:27 -0800 Subject: [PATCH 11/16] Clean out another duplicate ref. --- docs/ocean.bib | 12 ------------ docs/zotero.bib | 10 +++++----- src/parameterizations/lateral/MOM_hor_visc.F90 | 10 +++++----- .../lateral/MOM_interface_filter.F90 | 2 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 6 +++--- .../lateral/MOM_mixed_layer_restrat.F90 | 15 ++++++++------- .../lateral/MOM_self_attr_load.F90 | 4 ++-- 7 files changed, 24 insertions(+), 35 deletions(-) diff --git a/docs/ocean.bib b/docs/ocean.bib index ec35116efc..33107fbae1 100644 --- a/docs/ocean.bib +++ b/docs/ocean.bib @@ -10,18 +10,6 @@ @article{Adcroft2004 journal = {Ocean Modelling} } -@article{Adcroft2019, - doi = {10.1029/2019ms001726}, - year = 2019, - publisher = {American Geophysical Union ({AGU})}, - volume = {11}, - number = {10}, - pages = {3167--3211}, - author = {A. Adcroft and W. Anderson and V. Balaji and C. Blanton and M. Bushuk and C. O. Dufour and J. P. Dunne and S. M. Griffies and R. Hallberg and M. J. Harrison and I. M. Held and M. F. Jansen and J. G. John and J. P. Krasting and A. R. Langenhorst and S. Legg and Z. Liang and C. McHugh and A. Radhakrishnan and B. G. Reichl and T. Rosati and B. L. Samuels and A. Shao and R. Stouffer and M. Winton and A. T. Wittenberg and B. Xiang and N. Zadeh and R. Zhang}, - title = {The {GFDL} Global Ocean and Sea Ice Model {OM}4.0: Model Description and Simulation Features}, - journal = {J. Adv. Mod. Earth Sys.} -} - @article{Campin2004, doi = {10.1016/s1463-5003(03)00009-x}, year = 2004, diff --git a/docs/zotero.bib b/docs/zotero.bib index abaeea4c72..1120771b8a 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -967,7 +967,7 @@ @article{bitz1999 pages = {15669--15677} } -@inproceedings{briegleb2007, +@incollection{briegleb2007, series = {Technical {Note}}, title = {A {Delta}-{Eddington} {Mutiple} {Scattering} {Parameterization} for {Solar} {Radiation} in the {Sea} {Ice} {Component} of the {Community} {Climate} {System} {Model} {\textbar} {OpenSky} {Repository}}, url = {http://opensky.ucar.edu/islandora/object/technotes:484}, @@ -2447,7 +2447,7 @@ @article{russell1981 doi = {10.1175/1520-0450(1981)020<1483:ANFDSF>2.0.CO;2} } -@inproceedings{huynh1997, +@incollection{huynh1997, title = {Schemes and constraints for advection}, booktitle = {Fifteenth International Conference on Numerical Methods in Fluid Dynamics}, @@ -2564,7 +2564,7 @@ @article{marshall2010 doi = {10.1016/j.ocemod.2010.02.001} } -@inproceedings{millero1978, +@incollection{millero1978, author = {Millero, F.J.}, title = {Freezing point of seawater}, note = {Annex 6}, @@ -2820,7 +2820,7 @@ @article{Large2001 pages={518–536} } -@inproceedings{Smagorinsky1993, +@incollection{Smagorinsky1993, author={Joseph Smagorinsky}, year={1993}, title={Some historical remarks on the use of non-linear viscosities}, @@ -2884,7 +2884,7 @@ @article{Wang2012-2 pages={190--199} } -@inproceedings{Hallberg2003, +@incollection{Hallberg2003, title={The ability of large-scale ocean models to accept parameterizations of boundary mixing, and a description of a refined bulk mixed-layer model}, author={Robert Hallberg}, year={2003}, diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 4a901beab0..94afd0c858 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2044,9 +2044,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, end subroutine horizontal_viscosity -!> Allocates space for and calculates static variables used by horizontal_viscosity(). +!> Allocates space for and calculates static variables used by horizontal_viscosity. !! hor_visc_init calculates and stores the values of a number of metric functions that -!! are used in horizontal_viscosity(). +!! are used in horizontal_viscosity. subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) type(time_type), intent(in) :: Time !< Current model time. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. @@ -3259,7 +3259,7 @@ end subroutine hor_visc_end !! \f} !! !! The viscosity \f$\kappa_h\f$ may either be a constant or variable. For example, -!! \f$\kappa_h\f$ may vary with the shear, as proposed by Smagorinsky (1993). +!! \f$\kappa_h\f$ may vary with the shear, as proposed by \cite Smagorinsky1993. !! !! The accelerations resulting form the divergence of the stress tensor are !! \f{eqnarray*}{ @@ -3357,8 +3357,8 @@ end subroutine hor_visc_end !! !! \subsection section_anisotropic_viscosity Anisotropic viscosity !! -!! \cite Large2001, proposed enhancing viscosity in a particular direction and the -!! approach was generalized in Smith and McWilliams, 2003. We use the second form of their +!! \cite Large2001 proposed enhancing viscosity in a particular direction and the +!! approach was generalized in \cite Smith2003. We use the second form of their !! two coefficient anisotropic viscosity (section 4.3). We also replace their !! \f$A^\prime\f$ and $D$ such that \f$2A^\prime = 2 \kappa_h + D\f$ and !! \f$\kappa_a = D\f$ so that \f$\kappa_h\f$ can be considered the isotropic diff --git a/src/parameterizations/lateral/MOM_interface_filter.F90 b/src/parameterizations/lateral/MOM_interface_filter.F90 index ea86b80594..a63a2fc141 100644 --- a/src/parameterizations/lateral/MOM_interface_filter.F90 +++ b/src/parameterizations/lateral/MOM_interface_filter.F90 @@ -480,7 +480,7 @@ end subroutine interface_filter_end !! filter, depending on the order of the filter, or to the slope for a Laplacian !! filter !! \f[ -!! \vec{\psi} = - \kappa_h {\nabla \eta - \eta_smooth} +!! \vec{\psi} = - \kappa_h {\nabla \eta - \eta_{smooth}} !! \f] !! !! The result of the above expression is subsequently bounded by minimum and maximum values, including a diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 1a9afad897..eb1cf7b1cf 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1734,8 +1734,8 @@ end subroutine VarMix_end !! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p} !! \f] !! -!! The resolution function can be applied independently to thickness diffusion (module mom_thickness_diffuse), -!! tracer diffusion (mom_tracer_hordiff) lateral viscosity (mom_hor_visc). +!! The resolution function can be applied independently to thickness diffusion \(module mom_thickness_diffuse\), +!! tracer diffusion \(mom_tracer_hordiff\) lateral viscosity \(mom_hor_visc\). !! !! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. !! Ocean Modelling, 71, pp 92-103. http://dx.doi.org/10.1016/j.ocemod.2013.08.007 @@ -1757,7 +1757,7 @@ end subroutine VarMix_end !! \section section_Vicbeck Visbeck diffusivity !! !! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997, -!! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse() but calculated in this module. +!! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse but calculated in this module. !! !! \f[ !! \kappa_h = \alpha_s L_s^2 S N diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 18c219bc7f..f5deea1f66 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -1971,15 +1971,15 @@ end function test_answer !! \cite Hallberg2003, which this new parameterization surpasses, which !! in turn is based on the sub-inertial mixed layer theory of \cite Young1994. !! There is no net horizontal volume transport due to this parameterization, and -!! no direct effect below the mixed layer. A revised of the parameterization by -!! \cite Bodner2023, is also available as an option. +!! no direct effect below the mixed layer. A revised version of the parameterization by +!! \cite Bodner2023 is also available as an option. !! !! This parameterization sets the restratification timescale to agree with !! high-resolution studies of mixed layer restratification. !! -!! The run-time parameter FOX_KEMPER_ML_RESTRAT_COEF is a non-dimensional number of +!! The run-time parameter FOX_KEMPER_ML_RESTRAT_COEF is a non-dimensional number of !! order a few tens, proportional to the ratio of the deformation radius or the -!! grid scale (whichever is smaller to the dominant horizontal length-scale of the +!! grid scale (whichever is smaller) to the dominant horizontal length-scale of the !! sub-meso-scale mixed layer instabilities. !! !! \subsection section_mle_nutshell "Sub-meso" in a nutshell @@ -2011,14 +2011,15 @@ end function test_answer !! \f$ \tau \f$ is a time-scale for mixing momentum across the mixed layer. !! \f$ l_f \f$ is thought to be of order hundreds of meters. !! -!! The upscaling factor \f$ \frac{\Delta s}{l_f} \f$ can be a global constant, model parameter FOX_KEMPER_ML_RESTRAT, -!! so that in practice the parameterization is: +!! The upscaling factor \f$ \frac{\Delta s}{l_f} \f$ can be a global constant, model parameter +!! FOX_KEMPER_ML_RESTRAT, so that in practice the parameterization is: !! \f[ !! {\bf \Psi} = C_e \Gamma_\Delta \frac{ H^2 \nabla \bar{b} \times \hat{\bf z} }{ \sqrt{ f^2 + \tau^{-2}} } \mu(z) !! \f] !! with non-unity \f$ \Gamma_\Delta \f$. !! !! \f$ C_e \f$ is hard-coded as 0.0625. \f$ \tau \f$ is calculated from the surface friction velocity \f$ u^* \f$. +!! !! \todo Explain expression for momentum mixing time-scale. !! !! | Symbol | Module parameter | @@ -2102,7 +2103,7 @@ end function test_answer !! | \f$ \tau_{h+} \f$ | MLE\%BLD_GROWING_TFILTER | !! | \f$ \tau_{h-} \f$ | MLE\%BLD_DECAYING_TFILTER | !! | \f$ \tau_{H+} \f$ | MLE\%MLD_GROWING_TFILTER | -!! | \f$ \tau_{H-} \f$ | MLE\%BLD_DECAYING_TFILTER | +!! | \f$ \tau_{H-} \f$ | MLE\%MLD_DECAYING_TFILTER | !! !! \subsection section_mle_ref References !! diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index df22c9494c..fd6c0c1e5f 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -251,12 +251,12 @@ end subroutine SAL_end !! (rather, it should be bottom pressure anomaly). SAL is primarily used for fast evolving processes like tides or !! storm surges, but the effect applies to all motions. !! -!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied (\cite Accad1978) and the SAL is simply +!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied [\cite Accad1978] and the SAL is simply !! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH . For tides, the scalar !! approximation can also be used to iterate the SAL to convergence [see USE_PREVIOUS_TIDES in MOM_tidal_forcing, !! \cite Arbic2004]. !! -!! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate SAL. +!! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate SAL. !! Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is set by !! SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean !! developed by Los Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. From 078ec30d80254f86423d893b530b50cc71ab7cc0 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 25 Jul 2024 14:25:43 -0800 Subject: [PATCH 12/16] Fix the Accad ref. --- docs/zotero.bib | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/docs/zotero.bib b/docs/zotero.bib index 1120771b8a..bbd2e30478 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2892,7 +2892,7 @@ @incollection{Hallberg2003 pages={187--203} } -@article{1978, +@article{Accad1978, volume={290}, ISSN={2054-0272}, url={http://dx.doi.org/10.1098/rsta.1978.0083}, @@ -2902,7 +2902,10 @@ @article{1978 publisher={The Royal Society}, year={1978}, month=nov, - pages={235-–266} + pages={235-–266}, + author={Accad, Y. and Pekeris, C.L.}, + title={Solution of the tidal equations for the M2 and S2 tides in the world oceans from a + knowledge of the tidal potential alone} } @article{Arbic2004, From da3ec99d2d99158e1695525f4e13b24514fd5a88 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 25 Jul 2024 15:21:45 -0800 Subject: [PATCH 13/16] still cleaning the docs --- docs/parameterizations_lateral.rst | 5 +---- .../lateral/MOM_self_attr_load.F90 | 15 ++++++++------- .../lateral/MOM_spherical_harmonics.F90 | 8 ++++---- 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 3a444098b2..ce223e5155 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -68,9 +68,6 @@ While the self attraction and loading is computed in MOM_self_attr_load: :ref:`namespaceself__attr__load_1section_SAL` -Spherical Harmonics -------------------- - -The model needs spherical, computed in MOM_spherical_harmonics: +The self attraction and loading needs spherical harmonics, computed in MOM_spherical_harmonics: :ref:`namespacemom__spherical__harmonics_1section_spherical_harmonics` diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index fd6c0c1e5f..f72b6f513e 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -251,14 +251,15 @@ end subroutine SAL_end !! (rather, it should be bottom pressure anomaly). SAL is primarily used for fast evolving processes like tides or !! storm surges, but the effect applies to all motions. !! -!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied [\cite Accad1978] and the SAL is simply -!! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH . For tides, the scalar -!! approximation can also be used to iterate the SAL to convergence [see USE_PREVIOUS_TIDES in MOM_tidal_forcing, -!! \cite Arbic2004]. +!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied (\cite Accad1978) and the SAL is simply +!! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH. +!! For tides, the scalar approximation can also be used to iterate the SAL to convergence [see +!! USE_PREVIOUS_TIDES in MOM_tidal_forcing, \cite Arbic2004]. !! -!! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate SAL. -!! Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is set by -!! SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean +!! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate +!! SAL. Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is +!! set by SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across +!! Scales (MPAS)-Ocean !! developed by Los Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. !! !! References: diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 index 1782b5ae67..4f9cee03aa 100644 --- a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 +++ b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 @@ -361,7 +361,7 @@ end function order2index !! \f[ !! f^m_n = \sum^{Nj}_{0}\sum^{Ni}_{0}f(i,j)Y^m_n(i,j)A(i,j)/r_e^2 !! \f] -!! where $A$ is the area of the cell and $r_e$ is the radius of the Earth. +!! where \f$A\f$ is the area of the cell and \f$r_e\f$ is the radius of the Earth. !! !! In inverse transform, the first N degree spherical harmonic coefficients are used to reconstruct a two-dimensional !! physical field: @@ -374,10 +374,10 @@ end function order2index !! array vectorization. !! !! The maximum degree of the spherical harmonics is a runtime parameter and the maximum used by all SHT applications. -!! At the moment, it is only decided by SAL_HARMONICS_DEGREE. +!! At the moment, it is only decided by SAL_HARMONICS_DEGREE. !! -!! The forward transforms involve a global summation. Runtime flag SHT_REPRODUCING_SUM controls whether this is done -!! in a bit-wise reproducing way or not. +!! The forward transforms involve a global summation. Runtime flag SHT_REPRODUCING_SUM controls +!! whether this is done in a bit-wise reproducing way or not. !! !! References: !! From b7fee4d14a1e042885ad5c308aa509e65d5ede52 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 25 Jul 2024 16:28:27 -0800 Subject: [PATCH 14/16] Fix case --- docs/parameterizations_lateral.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index ce223e5155..9f161d5a04 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -64,7 +64,7 @@ The Love numbers are stored internally in MOM_load_love_numbers: :ref:`namespacemom__load__love__numbers_1section_Love_numbers` -While the self attraction and loading is computed in MOM_self_attr_load: +while the self attraction and loading is computed in MOM_self_attr_load: :ref:`namespaceself__attr__load_1section_SAL` From 4f1ecf4dd6b6f895ac81d8da55b2dde651161663 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 25 Jul 2024 22:02:09 -0800 Subject: [PATCH 15/16] Cleaning up tides section --- docs/parameterizations_lateral.rst | 9 ++++++--- src/parameterizations/lateral/MOM_tidal_forcing.F90 | 12 +++++++----- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/docs/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 9f161d5a04..d175c7e8bb 100644 --- a/docs/parameterizations_lateral.rst +++ b/docs/parameterizations_lateral.rst @@ -56,9 +56,9 @@ See also :ref:`namespacemom__lateral__mixing__coeffs_1section_Resolution_Functio Tidal forcing ------------- -Astronomical tidal forcings and self-attraction and loading are implement in MOM_tidal_forcing. -Tides can also be added via an open boundary tidal specification, -see `OBC wiki page `_. +Astronomical tidal forcings and self-attraction and loading are implement in + + :ref:`namespacetidal__forcing_1section_tides` The Love numbers are stored internally in MOM_load_love_numbers: @@ -71,3 +71,6 @@ while the self attraction and loading is computed in MOM_self_attr_load: The self attraction and loading needs spherical harmonics, computed in MOM_spherical_harmonics: :ref:`namespacemom__spherical__harmonics_1section_spherical_harmonics` + +Tides can also be added via an open boundary tidal specification, +see `OBC wiki page `_. diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 740977c7b9..177becf84f 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -744,6 +744,8 @@ end subroutine tidal_forcing_end !> \namespace tidal_forcing !! +!! \section section_tides Tidal forcing +!! !! Code by Robert Hallberg, August 2005, based on C-code by Harper !! Simmons, February, 2003, in turn based on code by Brian Arbic. !! @@ -762,14 +764,14 @@ end subroutine tidal_forcing_end !! !! In addition, approaches to calculate self-attraction and loading !! due to tides (harmonics of astronomical forcing frequencies) -!! are provided. TIDAL_SAL_FROM_FILE can be set to read the phase and -!! amplitude of the tidal SAL. USE_PREVIOUS_TIDES may be useful in +!! are provided. TIDAL_SAL_FROM_FILE can be set to read the phase and +!! amplitude of the tidal SAL. USE_PREVIOUS_TIDES may be useful in !! combination with the scalar approximation to iterate the SAL to !! convergence (for details, see \cite Arbic2004). With -!! TIDAL_SAL_FROM_FILE or USE_PREVIOUS_TIDES, a list of input files -!! must be provided to describe each constituent's properties from +!! TIDAL_SAL_FROM_FILE or USE_PREVIOUS_TIDES, a list of input +!! files must be provided to describe each constituent's properties from !! a previous solution. The online SAL calculations that are functions !! of SSH (rather should be bottom pressure anmoaly), either a scalar !! approximation or with spherical harmonic transforms, are located in -!! MOM_self_attr_load. +!! MOM_self_attr_load. end module MOM_tidal_forcing From 1b39d07be2787787907db6aa382dc61e3484709d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 23 Jul 2024 14:38:26 -0400 Subject: [PATCH 16/16] *+Fix non-Boussinesq MASS_WEIGHT_IN_PGF bug Optionally corrected a bug (essentially a sign error) in the selection of where MASS_WEIGHT_IN_PRESSURE_GRADIENT is applied in non-Boussinesq test cases. The (incorrect) non-Boussinesq version of the test for where interfaces are outside of the range of hydrostatic consistency due to interactions with bathymetry was converted from the (correct) Boussinesq version without properly taking into account that pressure increases downward while height increases upward. This bug is repeated 12 times in 6 different specific volume integral routines, including in the analytical specific volume integrals with 4 equations of state. To accommodate these corrections as well as a future expansion of the mass weighting to apply near the surface under ice-shelves or icebergs, the previous optional logical argument (useMassWghtInterp) was replaced with a new optional integer argument (MassWghtInterp) that can be used to encode information about where this weighting is applied as well as whether to fix this bug. This combination of settings (MASS_WEIGHT_IN_PRESSURE_GRADIENT = True and non-Boussinesq) does not seem to be widely used yet (it is not used in the GFDL regression suite), so rather than preserving the old (incorrect) solutions by default, this bug is corrected by default. However, the previous answers can be recovered by setting the new runtime parameter MASS_WEIGHT_IN_PGF_NONBOUS_BUG to true. This parameter that is only used (and logged) for non-Boussinesq cases with MASS_WEIGHT_IN_PRESSURE_GRADIENT set to true. It is intended that this new runtime bug-fix parameter will be obsoleted with an aggressive schedule if it is not needed to recreate any production runs. In addition, there are two new diagnostics, MassWt_u and MassWt_v, that show the fractional (0 to 1) effects of the mass weighting in the pressure gradient forces at the velocity points, along with the new diagnostic subroutines diagnose_mass_weight_Z and diagnose_mass_weight_p that calculate these diagnostics by layer in code that replicates the 13 copies for various equations of state and vertical structures within layers. By default, this commit does change answers in non-Boussinesq cases that use MASS_WEIGHT_IN_PRESSURE_GRADIENT = True, and there is a new runtime parameter in such cases. There are also two new diagnostics. Answers are bitwise identical in all Boussinesq cases. --- src/core/MOM_PressureForce_FV.F90 | 60 ++++- src/core/MOM_density_integrals.F90 | 244 ++++++++++++++---- src/equation_of_state/MOM_EOS.F90 | 36 +-- src/equation_of_state/MOM_EOS_Wright.F90 | 50 ++-- src/equation_of_state/MOM_EOS_Wright_full.F90 | 50 ++-- src/equation_of_state/MOM_EOS_Wright_red.F90 | 50 ++-- src/equation_of_state/MOM_EOS_linear.F90 | 50 ++-- 7 files changed, 378 insertions(+), 162 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 5d0f4ff167..9c4f355692 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -20,6 +20,7 @@ module MOM_PressureForce_FV use MOM_density_integrals, only : int_density_dz_generic_plm, int_density_dz_generic_ppm use MOM_density_integrals, only : int_spec_vol_dp_generic_plm use MOM_density_integrals, only : int_density_dz_generic_pcm, int_spec_vol_dp_generic_pcm +use MOM_density_integrals, only : diagnose_mass_weight_Z, diagnose_mass_weight_p use MOM_ALE, only : TS_PLM_edge_values, TS_PPM_edge_values, ALE_CS implicit none ; private @@ -46,7 +47,7 @@ module MOM_PressureForce_FV type(time_type), pointer :: Time !< A pointer to the ocean model's clock. type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. - logical :: useMassWghtInterp !< Use mass weighting in T/S interpolation + integer :: MassWghtInterp !< A flag indicating whether and how to use mass weighting in T/S interpolation logical :: use_inaccurate_pgf_rho_anom !< If true, uses the older and less accurate !! method to calculate density anomalies, as used prior to !! March 2018. @@ -71,6 +72,8 @@ module MOM_PressureForce_FV integer :: id_rho_pgf = -1 !< Diagnostic identifier integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier integer :: id_p_stanley = -1 !< Diagnostic identifier + integer :: id_MassWt_u = -1 !< Diagnostic identifier + integer :: id_MassWt_v = -1 !< Diagnostic identifier type(SAL_CS), pointer :: SAL_CSp => NULL() !< SAL control structure type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure end type PressureForce_FV_CS @@ -146,6 +149,10 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! interfaces, divided by the grid spacing [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & inty_dza ! The change in inty_za through a layer [L2 T-2 ~> m2 s-2]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & + MassWt_u ! The fractional mass weighting at a u-point [nondim]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & + MassWt_v ! The fractional mass weighting at a v-point [nondim]. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). @@ -194,6 +201,10 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ alpha_ref = 1.0 / CS%Rho0 I_gEarth = 1.0 / GV%g_Earth + if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then + MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0 + endif + if (use_p_atm) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -263,7 +274,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), & p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, & tv%eqn_of_state, US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), & - useMassWghtInterp=CS%useMassWghtInterp) + MassWghtInterp=CS%MassWghtInterp) elseif ( CS%Recon_Scheme == 2 ) then call MOM_error(FATAL, "PressureForce_FV_nonBouss: "//& "int_spec_vol_dp_generic_ppm does not exist yet.") @@ -277,8 +288,11 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, & US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), & inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, & - useMassWghtInterp=CS%useMassWghtInterp) + MassWghtInterp=CS%MassWghtInterp) endif + if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) & + call diagnose_mass_weight_p(p(:,:,K), p(:,:,K+1), dp_neglect, p(:,:,nz+1), G%HI, & + MassWt_u(:,:,k), MassWt_v(:,:,k)) else alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -468,6 +482,8 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag) if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag) + if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag) + if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag) end subroutine PressureForce_FV_nonBouss @@ -543,6 +559,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! of salinity within each layer [S ~> ppt]. T_t, T_b ! Top and bottom edge values for linear reconstructions ! of temperature within each layer [C ~> degC]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & + MassWt_u ! The fractional mass weighting at a u-point [nondim]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & + MassWt_v ! The fractional mass weighting at a v-point [nondim]. real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & rho_pgf, rho_stanley_pgf ! Density [R ~> kg m-3] from EOS with and without SGS T variance ! in Stanley parameterization. @@ -591,6 +611,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm G_Rho0 = GV%g_Earth / GV%Rho0 rho_ref = CS%Rho0 + if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then + MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0 + endif + if (CS%tides_answer_date>20230630) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = -G%bathyT(i,j) @@ -744,20 +768,20 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), & intx_dpa(:,:,k), inty_dpa(:,:,k), & - useMassWghtInterp=CS%useMassWghtInterp, & + MassWghtInterp=CS%MassWghtInterp, & use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), & intx_dpa(:,:,k), inty_dpa(:,:,k), & - useMassWghtInterp=CS%useMassWghtInterp, Z_0p=G%Z_ref) + MassWghtInterp=CS%MassWghtInterp, Z_0p=G%Z_ref) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa(:,:,k), & intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, dz_neglect, & - CS%useMassWghtInterp, Z_0p=G%Z_ref) + CS%MassWghtInterp, Z_0p=G%Z_ref) endif if (GV%Z_to_H /= 1.0) then !$OMP parallel do default(shared) @@ -765,6 +789,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm intz_dpa(i,j,k) = intz_dpa(i,j,k)*GV%Z_to_H enddo ; enddo endif + if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) & + call diagnose_mass_weight_Z(e(:,:,K), e(:,:,K+1), dz_neglect, G%bathyT, G%HI, & + MassWt_u(:,:,k), MassWt_v(:,:,k)) else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -955,6 +982,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag) if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag) + if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag) + if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag) if (CS%id_rho_pgf>0) call post_data(CS%id_rho_pgf, rho_pgf, CS%diag) if (CS%id_rho_stanley_pgf>0) call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag) @@ -978,6 +1007,8 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale ! temperature variance [nondim] integer :: default_answer_date ! Global answer date + logical :: useMassWghtInterp ! If true, use near-bottom mass weighting for T and S + logical :: MassWghtInterp_NonBous_bug ! If true, use a buggy mass weighting when non-Boussinesq ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl ! This module's name. @@ -1014,10 +1045,20 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, & "If True, use the ALE algorithm (regridding/remapping). "//& "If False, use the layered isopycnal algorithm.", default=.false. ) - call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", CS%useMassWghtInterp, & + call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", useMassWghtInterp, & "If true, use mass weighting when interpolating T/S for "//& "integrals near the bathymetry in FV pressure gradient "//& "calculations.", default=.false.) + call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_NONBOUS_BUG", MassWghtInterp_NonBous_bug, & + "If true, use a masking bug in non-Boussinesq calculations with mass weighting "//& + "when interpolating T/S for integrals near the bathymetry in FV pressure "//& + "gradient calculations.", & + default=.false., do_not_log=(GV%Boussinesq .or. (.not.useMassWghtInterp))) + CS%MassWghtInterp = 0 + if (useMassWghtInterp) & + CS%MassWghtInterp = ibset(CS%MassWghtInterp, 0) ! Same as CS%MassWghtInterp + 1 + if ((.not.GV%Boussinesq) .and. MassWghtInterp_NonBous_bug) & + CS%MassWghtInterp = ibset(CS%MassWghtInterp, 3) ! Same as CS%MassWghtInterp + 8 call get_param(param_file, mdl, "USE_INACCURATE_PGF_RHO_ANOM", CS%use_inaccurate_pgf_rho_anom, & "If true, use a form of the PGF that uses the reference density "//& "in an inaccurate way. This is not recommended.", default=.false.) @@ -1068,6 +1109,11 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, 'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m) endif + CS%id_MassWt_u = register_diag_field('ocean_model', 'MassWt_u', diag%axesCuL, Time, & + 'The fractional mass weighting at u-point PGF calculations', 'nondim') + CS%id_MassWt_v = register_diag_field('ocean_model', 'MassWt_v', diag%axesCvL, Time, & + 'The fractional mass weighting at v-point PGF calculations', 'nondim') + CS%GFS_scale = 1.0 if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 6d80d4dd55..d96116ba0c 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -31,6 +31,7 @@ module MOM_density_integrals public int_spec_vol_dp_generic_plm public avg_specific_vol public find_depth_of_pressure_in_cell +public diagnose_mass_weight_Z, diagnose_mass_weight_p contains @@ -39,7 +40,7 @@ module MOM_density_integrals !! required for calculating the finite-volume form pressure accelerations in a !! Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the arrays real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -77,16 +78,16 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] if (EOS_quadrature(EOS)) then call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p=Z_0p) else call analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif end subroutine int_density_dz @@ -96,7 +97,7 @@ end subroutine int_density_dz !! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, use_inaccurate_form, Z_0p) + dz_neglect, MassWghtInterp, use_inaccurate_form, Z_0p) type(hor_index_type), intent(in) :: HI !< Horizontal index type for input variables. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature of the layer [C ~> degC] @@ -134,8 +135,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of !! density anomalies, as was used prior to March 2018. real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] @@ -190,13 +191,13 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & endif do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (do_massWeight) then if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - "bathyT must be present if useMassWghtInterp is present and true.") + "bathyT must be present if MassWghtInterp is present and true.") if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - "dz_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif + "dz_neglect must be present if MassWghtInterp is present and true.") + endif ! Set the loop ranges for equation of state calculations at various points. EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2) @@ -368,7 +369,7 @@ end subroutine int_density_dz_generic_pcm !! T and S are linear profiles. subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, dpa, & - intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, & + intz_dpa, intx_dpa, inty_dpa, MassWghtInterp, & use_inaccurate_form, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays @@ -409,8 +410,8 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of !! density anomalies, as was used prior to March 2018. real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] @@ -481,13 +482,9 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & I_Rho = 1.0 / rho_0 z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. - if (present(useMassWghtInterp)) then - if (useMassWghtInterp) massWeightToggle = 1. - endif + if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif use_rho_ref = .true. - if (present(use_inaccurate_form)) then - if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form - endif + if (present(use_inaccurate_form)) use_rho_ref = .not. use_inaccurate_form use_varT = .false. !ensure initialized use_covarTS = .false. @@ -773,7 +770,7 @@ end subroutine int_density_dz_generic_plm !! are parabolic profiles subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, & - dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, Z_0p) + dpa, intz_dpa, intx_dpa, inty_dpa, MassWghtInterp, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -813,8 +810,8 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! This subroutine calculates (by numerical quadrature) integrals of @@ -884,9 +881,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & I_Rho = 1.0 / rho_0 z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. - if (present(useMassWghtInterp)) then - if (useMassWghtInterp) massWeightToggle = 1. - endif + if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif ! In event PPM calculation is bypassed with use_PPM=False s6 = 0. @@ -1179,7 +1174,7 @@ end subroutine int_density_dz_generic_ppm !! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) + bathyP, dP_tiny, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index structure real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1215,17 +1210,17 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_tiny !< A minuscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals if (EOS_quadrature(EOS)) then call int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) + bathyP, dP_tiny, MassWghtInterp) else call analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) + bathyP, dP_tiny, MassWghtInterp) endif end subroutine int_specific_vol_dp @@ -1237,7 +1232,7 @@ end subroutine int_specific_vol_dp !! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, dza, & intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, useMassWghtInterp) + bathyP, dP_neglect, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature of the layer [C ~> degC] @@ -1273,9 +1268,9 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A minuscule pressure change with - !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + !! the same units as p_t [R L2 T-2 ~> Pa] + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -1308,6 +1303,7 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state @@ -1320,14 +1316,15 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh); endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh); endif - do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. + do_massWeight = .false. ; massWeight_bug = .false. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + if (do_massWeight) then if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& - "bathyP must be present if useMassWghtInterp is present and true.") + "bathyP must be present if MassWghtInterp is present and true.") if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& - "dP_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif + "dP_neglect must be present if MassWghtInterp is present and true.") + endif ! Set the loop ranges for equation of state calculations at various points. EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(ieh-ish+1) @@ -1366,8 +1363,11 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -1421,8 +1421,11 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect @@ -1478,7 +1481,7 @@ end subroutine int_spec_vol_dp_generic_pcm !! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, & dP_neglect, bathyP, HI, EOS, US, dza, & - intp_dza, intx_dza, inty_dza, useMassWghtInterp) + intp_dza, intx_dza, inty_dza, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T_t !< Potential temperature at the top of the layer [C ~> degC] @@ -1518,8 +1521,8 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, optional, intent(inout) :: inty_dza !< The integral in y of the difference between !! the geopotential anomaly at the top and bottom of the layer divided !! by the y grid spacing [L2 T-2 ~> m2 s-2] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -1556,6 +1559,7 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! 5 sub-column locations [L2 T-2 ~> m2 s-2] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state @@ -1563,8 +1567,9 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB - do_massWeight = .false. - if (present(useMassWghtInterp)) do_massWeight = useMassWghtInterp + do_massWeight = .false. ; massWeight_bug = .false. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set do n = 1, 5 ! Note that these are reversed from int_density_dz. wt_t(n) = 0.25 * real(n-1) @@ -1607,8 +1612,11 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! weighting. Note: To work in terrain following coordinates we could ! offset this distance by the layer thickness to replicate other models. hWght = 0.0 - if (do_massWeight) & - hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + if (do_massWeight .and. massWeight_bug) then + hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -1668,8 +1676,11 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! hydrostatic consistency. For large hWght we bias the interpolation ! of T,S along the top and bottom integrals, like thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect @@ -1726,6 +1737,133 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, end subroutine int_spec_vol_dp_generic_plm +!> Diagnose the fractional mass weighting in a layer that might be used with a Boussinesq calculation. +subroutine diagnose_mass_weight_Z(z_t, z_b, dz_neglect, bathyT, HI, MassWt_u, MassWt_v) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m] + real, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] + + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] + real, dimension(SZIB_(HI),SZJ_(HI)), & + optional, intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim] + real, dimension(SZI_(HI),SZJB_(HI)), & + optional, intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim] + + ! Local variables + real :: hWght ! A pressure-thickness below topography [Z ~> m] + real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m] + real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2] + integer :: Isq, Ieq, Jsq, Jeq, i, j + + Isq = HI%IscB ; Ieq = HI%IecB + Jsq = HI%JscB ; Jeq = HI%JecB + + if (present(MassWt_u)) then + do j=HI%jsc,HI%jec ; do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (hWght > 0.) then + hL = (z_t(i,j) - z_b(i,j)) + dz_neglect + hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + MassWt_u(I,j) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_u(I,j) = 0.0 + endif + enddo ; enddo + endif + + if (present(MassWt_v)) then + do J=Jsq,Jeq ; do i=HI%isc,HI%iec + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (hWght > 0.) then + hL = (z_t(i,j) - z_b(i,j)) + dz_neglect + hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + MassWt_v(i,J) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_v(i,J) = 0.0 + endif + enddo ; enddo + endif + +end subroutine diagnose_mass_weight_Z + + +!> Diagnose the fractional mass weighting in a layer that might be used with a non-Boussinesq calculation. +subroutine diagnose_mass_weight_p(p_t, p_b, dP_neglect, bathyP, HI, MassWt_u, MassWt_v) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] + real, intent(in) :: dP_neglect ! Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(SZIB_(HI),SZJ_(HI)), & + optional, intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim] + real, dimension(SZI_(HI),SZJB_(HI)), & + optional, intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim] + + ! Local variables + real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa] + real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa] + real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2] + integer :: Isq, Ieq, Jsq, Jeq, i, j + + Isq = HI%IscB ; Ieq = HI%IecB + Jsq = HI%JscB ; Jeq = HI%JecB + + if (present(MassWt_u)) then + do j=HI%jsc,HI%jec ; do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + if (hWght > 0.) then + hL = (p_b(i,j) - p_t(i,j)) + dP_neglect + hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + MassWt_u(I,j) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_u(I,j) = 0.0 + endif + enddo ; enddo + endif + + if (present(MassWt_v)) then + do J=Jsq,Jeq ; do i=HI%isc,HI%iec + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + if (hWght > 0.) then + hL = (p_b(i,j) - p_t(i,j)) + dP_neglect + hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + MassWt_v(i,J) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_v(i,J) = 0.0 + endif + enddo ; enddo + endif + +end subroutine diagnose_mass_weight_p + !> Find the depth at which the reconstructed pressure matches P_tgt subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, & rho_ref, G_e, EOS, US, P_b, z_out, z_tol) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 7a9de49573..0180cdd2c0 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1226,7 +1226,7 @@ end function EOS_domain !! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) + bathyP, dP_tiny, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1261,8 +1261,8 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_tiny !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals ! Local variables real :: dRdT_scale ! A factor to convert drho_dT to the desired units [R degC m3 C-1 kg-1 ~> 1] @@ -1280,20 +1280,20 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%kg_m3_to_R*EOS%Rho_T0_S0, & dRdT_scale*EOS%dRho_dT, dRdS_scale*EOS%dRho_dS, dza, & intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) + bathyP, dP_tiny, MassWghtInterp) case (EOS_WRIGHT) call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & + inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, & SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt) case (EOS_WRIGHT_FULL) call int_spec_vol_dp_wright_full(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & + inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, & SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt) case (EOS_WRIGHT_REDUCED) call int_spec_vol_dp_wright_red(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & + inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, & SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt) case default @@ -1306,7 +1306,7 @@ end subroutine analytic_int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1343,8 +1343,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! Local variables @@ -1366,11 +1366,11 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (dRdT_scale /= 1.0) .or. (dRdS_scale /= 1.0)) then call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & rho_scale*EOS%Rho_T0_S0, dRdT_scale*EOS%dRho_dT, dRdS_scale*EOS%dRho_dS, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp) else call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp) endif case (EOS_WRIGHT) rho_scale = EOS%kg_m3_to_R @@ -1378,12 +1378,12 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, rho_scale, pres_scale, & + dz_neglect, MassWghtInterp, rho_scale, pres_scale, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p) else call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, Z_0p=Z_0p) + dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif case (EOS_WRIGHT_FULL) rho_scale = EOS%kg_m3_to_R @@ -1391,12 +1391,12 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then call int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, rho_scale, pres_scale, & + dz_neglect, MassWghtInterp, rho_scale, pres_scale, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p) else call int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, Z_0p=Z_0p) + dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif case (EOS_WRIGHT_REDUCED) rho_scale = EOS%kg_m3_to_R @@ -1404,12 +1404,12 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then call int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, rho_scale, pres_scale, & + dz_neglect, MassWghtInterp, rho_scale, pres_scale, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p) else call int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, Z_0p=Z_0p) + dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif case default call MOM_error(FATAL, "No analytic integration option is available with this EOS!") diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 8b6d6495d1..e71ae1791d 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -380,7 +380,7 @@ end subroutine EoS_fit_range_buggy_Wright !! that assumes that |eps| < 0.34. subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, & - useMassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) + MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -417,8 +417,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure @@ -517,13 +517,13 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & endif ; endif do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if useMassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + ! if (do_massWeight) then + ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "bathyT must be present if MassWghtInterp is present and true.") + ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "dz_neglect must be present if MassWghtInterp is present and true.") + ! endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 al0_2d(i,j) = (a0 + a1s*T(i,j)) + a2s*S(i,j) @@ -640,7 +640,7 @@ end subroutine int_density_dz_wright !! that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, & - useMassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) + MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -678,8 +678,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure @@ -726,6 +726,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo @@ -762,14 +763,15 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. + do_massWeight = .false. ; massWeight_bug = .false. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set +! if (do_massWeight) then ! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if useMassWghtInterp is present and true.") +! "bathyP must be present if MassWghtInterp is present and true.") ! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif +! "dP_neglect must be present if MassWghtInterp is present and true.") +! endif ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) do j=jsh,jeh ; do i=ish,ieh @@ -794,8 +796,11 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -835,8 +840,11 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect diff --git a/src/equation_of_state/MOM_EOS_Wright_full.F90 b/src/equation_of_state/MOM_EOS_Wright_full.F90 index 31b82e6190..73956d18fd 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -394,7 +394,7 @@ end subroutine EoS_fit_range_Wright_full !! that assumes that |eps| < 0.34. subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, & - useMassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) + MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -431,8 +431,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure @@ -531,13 +531,13 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & endif ; endif do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if useMassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + ! if (do_massWeight) then + ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "bathyT must be present if MassWghtInterp is present and true.") + ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "dz_neglect must be present if MassWghtInterp is present and true.") + ! endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 al0_2d(i,j) = a0 + (a1s*T(i,j) + a2s*S(i,j)) @@ -653,7 +653,7 @@ end subroutine int_density_dz_wright_full !! that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, & - useMassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) + MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -691,8 +691,8 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure @@ -740,6 +740,7 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo @@ -776,14 +777,15 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. + do_massWeight = .false. ; massWeight_bug = .false. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set +! if (do_massWeight) then ! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if useMassWghtInterp is present and true.") +! "bathyP must be present if MassWghtInterp is present and true.") ! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif +! "dP_neglect must be present if MassWghtInterp is present and true.") +! endif ! alpha = (lambda + al0*(pressure + p0)) / (pressure + p0) do j=jsh,jeh ; do i=ish,ieh @@ -809,8 +811,11 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -851,8 +856,11 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect diff --git a/src/equation_of_state/MOM_EOS_Wright_red.F90 b/src/equation_of_state/MOM_EOS_Wright_red.F90 index 65bdb9e521..835e3ecd26 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -396,7 +396,7 @@ end subroutine EoS_fit_range_Wright_red !! that assumes that |eps| < 0.34. subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, & - useMassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) + MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -433,8 +433,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure @@ -533,13 +533,13 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & endif ; endif do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if useMassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + ! if (do_massWeight) then + ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "bathyT must be present if MassWghtInterp is present and true.") + ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "dz_neglect must be present if MassWghtInterp is present and true.") + ! endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 al0_2d(i,j) = a0 + (a1s*T(i,j) + a2s*S(i,j)) @@ -655,7 +655,7 @@ end subroutine int_density_dz_wright_red !! that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, & - useMassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) + MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -693,8 +693,8 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure @@ -742,6 +742,7 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo @@ -778,14 +779,15 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. + do_massWeight = .false. ; massWeight_bug = .false. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set +! if (do_massWeight) then ! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if useMassWghtInterp is present and true.") +! "bathyP must be present if MassWghtInterp is present and true.") ! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif +! "dP_neglect must be present if MassWghtInterp is present and true.") +! endif ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) do j=jsh,jeh ; do i=ish,ieh @@ -811,8 +813,11 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -853,8 +858,11 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + endif if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index 8984fbca88..4ecf525abc 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -258,7 +258,7 @@ end subroutine set_params_linear !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp) + bathyT, dz_neglect, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -300,8 +300,8 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to - !! interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals ! Local variables real :: rho_anom ! The density anomaly from rho_ref [R ~> kg m-3]. @@ -328,13 +328,13 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & js = HI%jsc ; je = HI%jec do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if useMassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + ! if (do_massWeight) then + ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "bathyT must be present if MassWghtInterp is present and true.") + ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& + ! "dz_neglect must be present if MassWghtInterp is present and true.") + ! endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dz = z_t(i,j) - z_b(i,j) @@ -429,7 +429,7 @@ end subroutine int_density_dz_linear !! model. Specific volume is assumed to vary linearly between adjacent points. subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, useMassWghtInterp) + bathyP, dP_neglect, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -471,8 +471,8 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting - !! to interpolate T/S for top and bottom integrals. + integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals ! Local variables real :: dRho_TS ! The density anomaly due to T and S [R ~> kg m-3] real :: alpha_anom ! The specific volume anomaly from 1/rho_ref [R-1 ~> m3 kg-1] @@ -488,6 +488,7 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_6 = 1.0/6.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]. integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo @@ -497,14 +498,15 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh) ; endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh) ; endif - do_massWeight = .false. - if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then - do_massWeight = .true. + do_massWeight = .false. ; massWeight_bug = .false. + if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set +! if (do_massWeight) then ! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if useMassWghtInterp is present and true.") +! "bathyP must be present if MassWghtInterp is present and true.") ! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if useMassWghtInterp is present and true.") - endif ; endif +! "dP_neglect must be present if MassWghtInterp is present and true.") +! endif do j=jsh,jeh ; do i=ish,ieh dp = p_b(i,j) - p_t(i,j) @@ -520,8 +522,11 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) + endif if (hWght <= 0.0) then dpL = p_b(i,j) - p_t(i,j) ; dpR = p_b(i+1,j) - p_t(i+1,j) @@ -565,8 +570,11 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & ! hydrostatic consistency. For large hWght we bias the interpolation of ! T & S along the top and bottom integrals, akin to thickness weighting. hWght = 0.0 - if (do_massWeight) & + if (do_massWeight .and. massWeight_bug) then hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then + hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) + endif if (hWght <= 0.0) then dpL = p_b(i,j) - p_t(i,j) ; dpR = p_b(i,j+1) - p_t(i,j+1)