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/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/parameterizations_lateral.rst b/docs/parameterizations_lateral.rst index 3a3266a2bb..d175c7e8bb 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 ----------- @@ -32,17 +34,43 @@ 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. + + :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). +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: + + :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` + +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/docs/zotero.bib b/docs/zotero.bib index c0f1ddccbb..bbd2e30478 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}, @@ -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.}, @@ -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}, @@ -2760,3 +2760,189 @@ @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} +} + +@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}, + 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} +} + +@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} +} + +@incollection{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-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}, + journal={Computers & Geosciences}, + volume={49}, + pages={190--199} +} + +@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}, + booktitle={Internal Gravity Waves and Small-Scale Turbulence: Proc.‘Aha Huliko ‘a Hawaiian Winter Workshop}, + pages={187--203} +} + +@article{Accad1978, + 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}, + 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, + 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} +} + +@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/core/MOM.F90 b/src/core/MOM.F90 index 8b34b3e6a3..9e358ab17e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -95,6 +95,7 @@ module MOM use MOM_forcing_type, only : homogenize_forcing, homogenize_mech_forcing use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end use MOM_grid, only : set_first_direction +use MOM_harmonic_analysis, only : HA_accum_FtF, HA_accum_FtSSH, harmonic_analysis_CS use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_hor_index, only : rotate_hor_index use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz @@ -389,6 +390,8 @@ module MOM !< Pointer to the control structure used for the mode-split RK2 dynamics type(MOM_dyn_split_RK2b_CS), pointer :: dyn_split_RK2b_CSp => NULL() !< Pointer to the control structure used for an alternate version of the mode-split RK2 dynamics + type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() + !< Pointer to the control structure for harmonic analysis type(thickness_diffuse_CS) :: thickness_diffuse_CSp !< Pointer to the control structure used for the isopycnal height diffusive transport. !! This is also common referred to as Gent-McWilliams diffusion @@ -911,6 +914,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS enddo ; enddo endif + if (associated(CS%HA_CSp)) call HA_accum_FtF(Time_Local, CS%HA_CSp) call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, & dt_therm_here, bbl_time_int, CS, & @@ -1020,6 +1024,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ssh(i,j) = CS%ssh_rint(i,j)*I_wt_ssh CS%ave_ssh_ibc(i,j) = ssh(i,j) enddo ; enddo + if (associated(CS%HA_CSp)) call HA_accum_FtSSH('ssh', ssh, Time_local, G, CS%HA_CSp) if (do_dyn) then call adjust_ssh_for_p_atm(CS%tv, G, GV, US, CS%ave_ssh_ibc, forces%p_surf_SSH, & CS%calc_rho_for_sea_lev) @@ -3247,13 +3252,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & allocate(eta(SZI_(G),SZJ_(G)), source=0.0) if (CS%use_alt_split) then call initialize_dyn_split_RK2b(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, & - G, GV, US, param_file, diag, CS%dyn_split_RK2b_CSp, restart_CSp, & + G, GV, US, param_file, diag, CS%dyn_split_RK2b_CSp, CS%HA_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) else call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, & - G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, restart_CSp, & + G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, CS%HA_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) 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_barotropic.F90 b/src/core/MOM_barotropic.F90 index a132200759..bc8fddbdde 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -16,6 +16,7 @@ module MOM_barotropic use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : mech_forcing use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, only : HA_accum_FtSSH, harmonic_analysis_CS use MOM_hor_index, only : hor_index_type use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, open_boundary_query @@ -289,6 +290,7 @@ module MOM_barotropic type(MOM_domain_type), pointer :: BT_Domain => NULL() !< Barotropic MOM domain type(hor_index_type), pointer :: debug_BT_HI => NULL() !< debugging copy of horizontal index_type type(SAL_CS), pointer :: SAL_CSp => NULL() !< Control structure for SAL + type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() !< Control structure for harmonic analysis logical :: module_is_initialized = .false. !< If true, module has been initialized integer :: isdw !< The lower i-memory limit for the wide halo arrays. @@ -2551,6 +2553,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, eta_out(i,j) = eta_wtd(i,j) * I_sum_wt_eta enddo ; enddo + ! Accumulator is updated at the end of every baroclinic time step. + ! Harmonic analysis will not be performed of a field that is not registered. + if (associated(CS%HA_CSp) .and. find_etaav) then + call HA_accum_FtSSH('ubt', ubt, CS%Time, G, CS%HA_CSp) + call HA_accum_FtSSH('vbt', vbt, CS%Time, G, CS%HA_CSp) + endif + if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post) if (id_clock_pass_post > 0) call cpu_clock_begin(id_clock_pass_post) if (G%nonblocking_updates) then @@ -4402,7 +4411,7 @@ end subroutine bt_mass_source !! barotropic calculation and initializes any barotropic fields that have not !! already been initialized. subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, & - restart_CS, calc_dtbt, BT_cont, SAL_CSp) + restart_CS, calc_dtbt, BT_cont, SAL_CSp, HA_CSp) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -4428,6 +4437,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, !! barotropic flow. type(SAL_CS), target, optional :: SAL_CSp !< A pointer to the control structure of the !! SAL module. + type(harmonic_analysis_CS), target, optional :: HA_CSp !< A pointer to the control structure of the + !! harmonic analysis module ! This include declares and sets the variable "version". # include "version_variable.h" @@ -4489,6 +4500,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, if (present(SAL_CSp)) then CS%SAL_CSp => SAL_CSp endif + if (present(HA_CSp)) then + CS%HA_CSp => HA_CSp + endif ! Read all relevant parameters and write them to the model log. call get_param(param_file, mdl, "SPLIT", CS%split, default=.true., do_not_log=.true.) 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/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index b545d99daa..11c3ff1873 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -47,6 +47,7 @@ module MOM_dynamics_split_RK2 use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, only : harmonic_analysis_CS use MOM_hor_index, only : hor_index_type use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS use MOM_hor_visc, only : hor_visc_init, hor_visc_end @@ -242,6 +243,8 @@ module MOM_dynamics_split_RK2 type(SAL_CS) :: SAL_CSp !> A pointer to the tidal forcing control structure type(tidal_forcing_CS) :: tides_CSp + !> A pointer to the harmonic analysis control structure + type(harmonic_analysis_CS) :: HA_CSp !> A pointer to the ALE control structure. type(ALE_CS), pointer :: ALE_CSp => NULL() @@ -482,7 +485,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f call cpu_clock_end(id_clock_pass) !--- end set up for group halo pass - ! PFu = d/dx M(h,T,S) ! pbce = dM/deta if (CS%begw == 0.0) call enable_averages(dt, Time_local, CS%diag) @@ -1308,7 +1310,7 @@ end subroutine remap_dyn_split_RK2_aux_vars !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & - diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & + diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil) @@ -1331,6 +1333,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p type(param_file_type), intent(in) :: param_file !< parameter file for parsing type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure + type(harmonic_analysis_CS), pointer :: HA_CSp !< A pointer to the control structure of the + !! harmonic analysis module type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure real, intent(in) :: dt !< time step [T ~> s] type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for @@ -1504,7 +1508,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) - if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp) + if (CS%use_tides) then + call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) + HA_CSp => CS%HA_CSp + else + HA_CSp => NULL() + endif call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & CS%SAL_CSp, CS%tides_CSp) call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp) @@ -1538,7 +1547,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, & - CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp) + CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp, CS%HA_CSp) if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .not. query_initialized(CS%diffv, "diffv", restart_CS)) then diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index c6642506f9..0220db7993 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -49,6 +49,7 @@ module MOM_dynamics_split_RK2b use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, only : harmonic_analysis_CS use MOM_hor_index, only : hor_index_type use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS use MOM_hor_visc, only : hor_visc_init, hor_visc_end @@ -239,6 +240,8 @@ module MOM_dynamics_split_RK2b type(SAL_CS) :: SAL_CSp !> A pointer to the tidal forcing control structure type(tidal_forcing_CS) :: tides_CSp + !> A pointer to the harmonic analysis control structure + type(harmonic_analysis_CS) :: HA_CSp !> A pointer to the ALE control structure. type(ALE_CS), pointer :: ALE_CSp => NULL() @@ -1233,7 +1236,7 @@ end subroutine remap_dyn_split_RK2b_aux_vars !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & - diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & + diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil) @@ -1256,6 +1259,8 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, type(param_file_type), intent(in) :: param_file !< parameter file for parsing type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics type(MOM_dyn_split_RK2b_CS), pointer :: CS !< module control structure + type(harmonic_analysis_CS), pointer :: HA_CSp !< A pointer to the control structure of the + !! harmonic analysis module type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure real, intent(in) :: dt !< time step [T ~> s] type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for @@ -1415,7 +1420,12 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) - if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp) + if (CS%use_tides) then + call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) + HA_CSp => CS%HA_CSp + else + HA_CSp => NULL() + endif call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & CS%SAL_CSp, CS%tides_CSp) call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp) @@ -1453,7 +1463,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, & CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, & - CS%SAL_CSp) + CS%SAL_CSp, CS%HA_CSp) flux_units = get_flux_units(GV) thickness_units = get_thickness_units(GV) diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 new file mode 100644 index 0000000000..76adad5c8e --- /dev/null +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -0,0 +1,452 @@ +!> Inline harmonic analysis (conventional) +module MOM_harmonic_analysis + +use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, get_date, increment_date, & + operator(+), operator(-), operator(<), operator(>), operator(>=) +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_file_parser, only : param_file_type, get_param +use MOM_io, only : file_exists, open_ASCII_file, READONLY_FILE, close_file, & + MOM_infra_file, vardesc, MOM_field, & + var_desc, create_MOM_file, SINGLE_FILE, MOM_write_field +use MOM_error_handler, only : MOM_mesg, MOM_error, NOTE + +implicit none ; private + +public HA_init, HA_register, HA_accum_FtF, HA_accum_FtSSH + +#include + +integer, parameter :: MAX_CONSTITUENTS = 10 !< The maximum number of tidal constituents + +!> The private control structure for storing the HA info of a particular field +type, private :: HA_type + character(len=16) :: key = "none" !< Name of the field of which harmonic analysis is to be performed + character(len=1) :: grid !< The grid on which the field is defined ('h', 'q', 'u', or 'v') + real :: old_time = -1.0 !< The time of the previous accumulating step [T ~> s] + real, allocatable :: ref(:,:) !< The initial field in arbitrary units [A] + real, allocatable :: FtSSH(:,:,:) !< Accumulator of (F' * SSH_in) in arbitrary units [A] + !>@{ Lower and upper bounds of input data + integer :: is, ie, js, je + !>@} +end type HA_type + +!> A linked list of control structures that store the HA info of different fields +type, private :: HA_node + type(HA_type) :: this !< Control structure of the current field in the list + type(HA_node), pointer :: next !< The list of other fields +end type HA_node + +!> The public control structure of the MOM_harmonic_analysis module +type, public :: harmonic_analysis_CS ; private + logical :: HAready = .true. !< If true, perform harmonic analysis + type(time_type) :: & + time_start, & !< Start time of harmonic analysis + time_end, & !< End time of harmonic analysis + time_ref !< Reference time (t = 0) used to calculate tidal forcing + real, dimension(MAX_CONSTITUENTS) :: & + freq, & !< The frequency of a tidal constituent [T-1 ~> s-1] + phase0 !< The phase of a tidal constituent at time 0 [rad] + real, allocatable :: FtF(:,:) !< Accumulator of (F' * F) for all fields [nondim] + integer :: nc !< The number of tidal constituents in use + integer :: length !< Number of fields of which harmonic analysis is to be performed + character(len=16) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent + character(len=255) :: path !< Path to directory where output will be written + type(unit_scale_type) :: US !< A dimensional unit scaling type + type(HA_node), pointer :: list => NULL() !< A linked list for storing the HA info of different fields +end type harmonic_analysis_CS + +contains + +!> This subroutine sets static variables used by this module and initializes CS%list. +!! THIS MUST BE CALLED AT THE END OF tidal_forcing_init. +subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name, CS) + type(time_type), intent(in) :: Time !< The current model time + type(time_type), intent(in) :: time_ref !< Reference time (t = 0) used to calculate tidal forcing + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + real, dimension(MAX_CONSTITUENTS), intent(in) :: freq !< The frequency of a tidal constituent [T-1 ~> s-1] + real, dimension(MAX_CONSTITUENTS), intent(in) :: phase0 !< The phase of a tidal constituent at time 0 [rad] + integer, intent(in) :: nc !< The number of tidal constituents in use + character(len=16), intent(in) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent + type(harmonic_analysis_CS), intent(out) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + type(HA_type) :: ha1 !< A temporary, null field used for initializing CS%list + real :: HA_start_time !< Start time of harmonic analysis [T ~> s] + real :: HA_end_time !< End time of harmonic analysis [T ~> s] + character(len=40) :: mdl="MOM_harmonic_analysis" !< This module's name + character(len=255) :: mesg + integer :: year, month, day, hour, minute, second + + ! Determine CS%time_start and CS%time_end + call get_param(param_file, mdl, "HA_START_TIME", HA_start_time, & + "Start time of harmonic analysis, in units of days after "//& + "the start of the current run segment. Must be smaller than "//& + "HA_END_TIME, otherwise harmonic analysis will not be performed. "//& + "If negative, |HA_START_TIME| determines the length of harmonic analysis, "//& + "and harmonic analysis will start |HA_START_TIME| days before HA_END_TIME, "//& + "or at the beginning of the run segment, whichever occurs later.", & + units="days", default=0.0, scale=86400.0*US%s_to_T) + call get_param(param_file, mdl, "HA_END_TIME", HA_end_time, & + "End time of harmonic analysis, in units of days after "//& + "the start of the current run segment. Must be positive "//& + "and smaller than the length of the currnet run segment, "//& + "otherwise harmonic analysis will not be performed.", & + units="days", default=0.0, scale=86400.0*US%s_to_T) + + if (HA_end_time <= 0.0) then + call MOM_mesg('MOM_harmonic_analysis: HA_END_TIME is zero or negative. '//& + 'Harmonic analysis will not be performed.') + CS%HAready = .false. ; return + endif + + if (HA_end_time <= HA_start_time) then + call MOM_mesg('MOM_harmonic_analysis: HA_END_TIME is smaller than or equal to HA_START_TIME. '//& + 'Harmonic analysis will not be performed.') + CS%HAready = .false. ; return + endif + + if (HA_start_time < 0.0) then + HA_start_time = HA_end_time + HA_start_time + if (HA_start_time <= 0.0) HA_start_time = 0.0 + endif + + CS%time_start = Time + real_to_time(US%T_to_s * HA_start_time) + CS%time_end = Time + real_to_time(US%T_to_s * HA_end_time) + + call get_date(Time, year, month, day, hour, minute, second) + write(mesg,*) "MOM_harmonic_analysis: run segment starts on ", year, month, day, hour, minute, second + call MOM_error(NOTE, trim(mesg)) + call get_date(CS%time_start, year, month, day, hour, minute, second) + write(mesg,*) "MOM_harmonic_analysis: harmonic analysis starts on ", year, month, day, hour, minute, second + call MOM_error(NOTE, trim(mesg)) + call get_date(CS%time_end, year, month, day, hour, minute, second) + write(mesg,*) "MOM_harmonic_analysis: harmonic analysis ends on ", year, month, day, hour, minute, second + call MOM_error(NOTE, trim(mesg)) + + ! Set path to directory where output will be written + call get_param(param_file, mdl, "HA_PATH", CS%path, & + "Path to output files for runtime harmonic analysis.", default="./") + + ! Populate some parameters of the control structure + CS%time_ref = time_ref + CS%freq = freq + CS%phase0 = phase0 + CS%nc = nc + CS%const_name = const_name + CS%length = 0 + CS%US = US + + allocate(CS%FtF(2*nc+1,2*nc+1), source=0.0) + + ! Initialize CS%list + allocate(CS%list) + CS%list%this = ha1 + nullify(CS%list%next) + +end subroutine HA_init + +!> This subroutine registers each of the fields on which HA is to be performed. +subroutine HA_register(key, grid, CS) + character(len=*), intent(in) :: key !< Name of the current field + character(len=1), intent(in) :: grid !< The grid on which the key field is defined + type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + type(HA_type) :: ha1 !< Control structure for the current field + type(HA_node), pointer :: tmp !< A temporary list to hold the current field + + if (.not. CS%HAready) return + + allocate(tmp) + ha1%key = trim(key) + ha1%grid = trim(grid) + tmp%this = ha1 + tmp%next => CS%list + CS%list => tmp + CS%length = CS%length + 1 + +end subroutine HA_register + +!> This subroutine accumulates the temporal basis functions in FtF. +!! The tidal constituents are those used in MOM_tidal_forcing, plus the mean (of zero frequency). +subroutine HA_accum_FtF(Time, CS) + type(time_type), intent(in) :: Time !< The current model time + type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + real :: now !< The relative time compared with tidal reference [T ~> s] + real :: cosomegat, sinomegat, ccosomegat, ssinomegat !< The components of the phase [nondim] + integer :: nc, c, icos, isin, cc, iccos, issin + + ! Exit the accumulator in the following cases + if (.not. CS%HAready) return + if (CS%length == 0) return + if (Time < CS%time_start) return + if (Time > CS%time_end) return + + nc = CS%nc + now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref) + + ! Accumulate FtF + CS%FtF(1,1) = CS%FtF(1,1) + 1.0 !< For the zero frequency + do c=1,nc + icos = 2*c + isin = 2*c+1 + cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) + sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) + CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat + CS%FtF(isin,1) = CS%FtF(isin,1) + sinomegat + CS%FtF(1,icos) = CS%FtF(icos,1) + CS%FtF(1,isin) = CS%FtF(isin,1) + do cc=c,nc + iccos = 2*cc + issin = 2*cc+1 + ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc)) + ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc)) + CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat + CS%FtF(icos,issin) = CS%FtF(icos,issin) + cosomegat * ssinomegat + CS%FtF(isin,iccos) = CS%FtF(isin,iccos) + sinomegat * ccosomegat + CS%FtF(isin,issin) = CS%FtF(isin,issin) + sinomegat * ssinomegat + enddo ! cc=c,nc + enddo ! c=1,nc + +end subroutine HA_accum_FtF + +!> This subroutine accumulates the temporal basis functions in FtSSH and then calls HA_write to compute +!! harmonic constants and write results. The tidal constituents are those used in MOM_tidal_forcing, plus the +!! mean (of zero frequency). +subroutine HA_accum_FtSSH(key, data, Time, G, CS) + character(len=*), intent(in) :: key !< Name of the current field + real, dimension(:,:), intent(in) :: data !< Input data of which harmonic analysis is to be performed [A] + type(time_type), intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + type(HA_type), pointer :: ha1 + type(HA_node), pointer :: tmp + real :: now !< The relative time compared with the tidal reference [T ~> s] + real :: dt !< The current time step size of the accumulator [T ~> s] + real :: cosomegat, sinomegat !< The components of the phase [nondim] + integer :: nc, i, j, k, c, icos, isin, is, ie, js, je + character(len=128) :: mesg + + ! Exit the accumulator in the following cases + if (.not. CS%HAready) return + if (CS%length == 0) return + if (Time < CS%time_start) return + if (Time > CS%time_end) return + + ! Loop through the full list to find the current field + tmp => CS%list + do k=1,CS%length + ha1 => tmp%this + if (trim(key) == trim(ha1%key)) exit + tmp => tmp%next + if (k == CS%length) return !< Do not perform harmonic analysis of a field that is not registered + enddo + + nc = CS%nc + now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref) + + ! Additional processing at the initial accumulating step + if (ha1%old_time < 0.0) then + ha1%old_time = now + + write(mesg,*) "MOM_harmonic_analysis: initializing accumulator, key = ", trim(ha1%key) + call MOM_error(NOTE, trim(mesg)) + + ! Get the lower and upper bounds of input data + ha1%is = LBOUND(data,1) ; is = ha1%is + ha1%ie = UBOUND(data,1) ; ie = ha1%ie + ha1%js = LBOUND(data,2) ; js = ha1%js + ha1%je = UBOUND(data,2) ; je = ha1%je + + allocate(ha1%ref(is:ie,js:je), source=0.0) + allocate(ha1%FtSSH(is:ie,js:je,2*nc+1), source=0.0) + ha1%ref(:,:) = data(:,:) + endif + + dt = now - ha1%old_time + ha1%old_time = now !< Keep track of time so we know when Time approaches CS%time_end + + is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + ! Accumulate FtF and FtSSH + do c=1,nc + icos = 2*c + isin = 2*c+1 + cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) + sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) + do j=js,je ; do i=is,ie + ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j)) + ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat + ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat + enddo ; enddo + enddo ! c=1,nc + + ! Compute harmonic constants and write output as Time approaches CS%time_end + ! This guarantees that HA_write will be called before Time becomes larger than CS%time_end + if (time_type_to_real(CS%time_end - Time) <= dt) then + call HA_write(ha1, Time, G, CS) + + write(mesg,*) "MOM_harmonic_analysis: harmonic analysis done, key = ", trim(ha1%key) + call MOM_error(NOTE, trim(mesg)) + + ! De-register the current field and deallocate memory + ha1%key = 'none' + deallocate(ha1%ref) + deallocate(ha1%FtSSH) + endif + +end subroutine HA_accum_FtSSH + +!> This subroutine computes the harmonic constants and write output for the current field +subroutine HA_write(ha1, Time, G, CS) + type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field + type(time_type), intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(harmonic_analysis_CS), intent(in) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + real, dimension(:,:,:), allocatable :: FtSSHw !< An array containing the harmonic constants [A] + integer :: year, month, day, hour, minute, second + integer :: nc, k, is, ie, js, je + + character(len=255) :: filename !< Output file name + type(MOM_infra_file) :: cdf !< The file handle for output harmonic constants + type(vardesc), allocatable :: cdf_vars(:) !< Output variable names + type(MOM_field), allocatable :: cdf_fields(:) !< Field type variables for the output fields + + nc = CS%nc ; is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0) + + ! Compute the harmonic coefficients + call HA_solver(ha1, nc, CS%FtF, FtSSHw) + + ! Output file name + call get_date(Time, year, month, day, hour, minute, second) + write(filename, '(a,"HA_",a,i0.4,i0.2,i0.2,".nc")') & + trim(CS%path), trim(ha1%key), year, month, day + + allocate(cdf_vars(2*nc+1)) + allocate(cdf_fields(2*nc+1)) + + ! Variable names + cdf_vars(1) = var_desc("z0", "m" ,"mean value", ha1%grid, '1', '1') + do k=1,nc + cdf_vars(2*k ) = var_desc(trim(CS%const_name(k))//"cos", "m", "cosine coefficient", ha1%grid, '1', '1') + cdf_vars(2*k+1) = var_desc(trim(CS%const_name(k))//"sin", "m", "sine coefficient", ha1%grid, '1', '1') + enddo + + ! Create output file + call create_MOM_file(cdf, trim(filename), cdf_vars, & + 2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=G) + + ! Write data + call MOM_write_field(cdf, cdf_fields(1), G%domain, FtSSHw(:,:,1), 0.0) + do k=1,nc + call MOM_write_field(cdf, cdf_fields(2*k ), G%domain, FtSSHw(:,:,2*k ), 0.0) + call MOM_write_field(cdf, cdf_fields(2*k+1), G%domain, FtSSHw(:,:,2*k+1), 0.0) + enddo + + call cdf%flush() + deallocate(cdf_vars) + deallocate(cdf_fields) + deallocate(FtSSHw) + +end subroutine HA_write + +!> This subroutine computes the harmonic constants (stored in FtSSHw) using the dot products of the temporal +!! basis functions accumulated in FtF, and the dot products of the SSH (or other fields) with the temporal basis +!! functions accumulated in FtSSH. The system is solved by Cholesky decomposition, +!! +!! FtF * FtSSHw = FtSSH, => FtFw * (FtFw' * FtSSHw) = FtSSH, +!! +!! where FtFw is a lower triangular matrix, and the prime denotes matrix transpose. +!! +subroutine HA_solver(ha1, nc, FtF, FtSSHw) + type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field + integer, intent(in) :: nc !< Number of harmonic constituents + real, dimension(:,:), intent(in) :: FtF !< Accumulator of (F' * F) for all fields [nondim] + real, dimension(:,:,:), allocatable, intent(out) :: FtSSHw !< Work array for Cholesky decomposition [A] + + ! Local variables + real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim] + real, dimension(:), allocatable :: tmp1 !< Temporary variable for Cholesky decomposition [nondim] + real, dimension(:,:), allocatable :: tmp2 !< Temporary variable for Cholesky decomposition [A] + real, dimension(:,:), allocatable :: FtFw !< Lower triangular matrix for Cholesky decomposition [nondim] + integer :: k, m, n, is, ie, js, je + + is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + allocate(tmp1(1:2*nc+1), source=0.0) + allocate(tmp2(is:ie,js:je), source=0.0) + allocate(FtFw(1:2*nc+1,1:2*nc+1), source=0.0) + allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0) + + ! Construct FtFw + FtFw(:,:) = 0.0 + do m=1,2*nc+1 + tmp0 = 0.0 + do k=1,m-1 + tmp0 = tmp0 + FtFw(m,k) * FtFw(m,k) + enddo + FtFw(m,m) = sqrt(FtF(m,m) - tmp0) + tmp1(m) = 1 / FtFw(m,m) + do k=m+1,2*nc+1 + tmp0 = 0.0 + do n=1,m-1 + tmp0 = tmp0 + FtFw(k,n) * FtFw(m,n) + enddo + FtFw(k,m) = (FtF(k,m) - tmp0) * tmp1(m) + enddo + enddo + + ! Solve for (FtFw' * FtSSHw) + FtSSHw(:,:,:) = ha1%FtSSH(:,:,:) + do k=1,2*nc+1 + tmp2(:,:) = 0.0 + do m=1,k-1 + tmp2(:,:) = tmp2(:,:) + FtFw(k,m) * FtSSHw(:,:,m) + enddo + FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k) + enddo + + ! Solve for FtSSHw + do k=2*nc+1,1,-1 + tmp2(:,:) = 0.0 + do m=k+1,2*nc+1 + tmp2(:,:) = tmp2(:,:) + FtSSHw(:,:,m) * FtFw(m,k) + enddo + FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k) + enddo + + deallocate(tmp1) + deallocate(tmp2) + deallocate(FtFw) + +end subroutine HA_solver + +!> \namespace harmonic_analysis +!! +!! This module computes the harmonic constants which can be used to reconstruct the tidal elevation (or other +!! fields) through SSH = F * x, where F is an nt-by-2*nc matrix (nt is the number of time steps and nc is the +!! number of tidal constituents) containing the cosine/sine functions for each frequency evaluated at each time +!! step, and x is a 2*nc-by-1 vector containing the constant coefficients of the sine/cosine for each constituent +!! (i.e., the harmonic constants). At each grid point, the harmonic constants are computed using least squares, +!! +!! (F' * F) * x = F' * SSH_in, +!! +!! where the prime denotes matrix transpose, and SSH_in is the sea surface height (or other fields) determined by +!! the model. The dot products (F' * F) and (F' * SSH_in) are computed by accumulating the sums as the model is +!! running and stored in the arrays FtF and FtSSH, respectively. The FtF matrix is inverted as needed before +!! computing and writing out the harmonic constants. +!! +!! Ed Zaron and William Xu (chengzhu.xu@oregonstate.edu), April 2024. + +end module MOM_harmonic_analysis + 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) 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 2a4181804b..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. @@ -3193,25 +3193,25 @@ subroutine hor_visc_end(CS) end subroutine hor_visc_end !> \namespace mom_hor_visc !! -!! This module contains the subroutine horizontal_viscosity() that calculates the +!! \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 +!! 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 -!! 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} = @@ -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 !! -!! Large et al., 2001, 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 7315c82601..eb1cf7b1cf 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1728,14 +1728,14 @@ 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} !! \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_load_love_numbers.F90 b/src/parameterizations/lateral/MOM_load_love_numbers.F90 index 3d573d894d..8faf3aafab 100644 --- a/src/parameterizations/lateral/MOM_load_love_numbers.F90 +++ b/src/parameterizations/lateral/MOM_load_love_numbers.F90 @@ -1452,16 +1452,18 @@ 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)], -!! as in subroutine calc_love_scaling in MOM_tidal_forcing module. +!! National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. The load Love numbers +!! 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. !! !! 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..f5deea1f66 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,29 +1965,29 @@ 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. +!! 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 !! !! 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) @@ -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 | @@ -2057,7 +2058,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} } @@ -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 e7f8a73ab2..f72b6f513e 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -245,19 +245,22 @@ 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 -!! 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)]. +!! 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 -!! developed by Los Alamos National Laboratory and University of Michigan [Barton et al. (2022) and Brus et al. (2023)]. +!! 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 26258e6b8e..4f9cee03aa 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: @@ -359,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: @@ -372,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: !! 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 1cd8a45a78..177becf84f 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -9,6 +9,8 @@ module MOM_tidal_forcing use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, & + only : HA_init, HA_register, harmonic_analysis_CS use MOM_io, only : field_exists, file_exists, MOM_read_data use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_unit_scaling, only : unit_scale_type @@ -231,12 +233,13 @@ end subroutine nodal_fu !! while fields like the background viscosities are 2-D arrays. !! ALLOC is a macro defined in MOM_memory.h for allocate or nothing with !! static memory. -subroutine tidal_forcing_init(Time, G, US, param_file, CS) +subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) type(time_type), intent(in) :: Time !< The current model time. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(tidal_forcing_CS), intent(inout) :: CS !< Tidal forcing control structure + type(harmonic_analysis_CS), optional, intent(out) :: HA_CS !< Control structure for harmonic analysis ! Local variables real, dimension(SZI_(G), SZJ_(G)) :: & @@ -251,6 +254,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS) logical :: use_M2, use_S2, use_N2, use_K2, use_K1, use_O1, use_P1, use_Q1 logical :: use_MF, use_MM logical :: tides ! True if a tidal forcing is to be used. + logical :: HA_ssh, HA_ubt, HA_vbt ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_tidal_forcing" ! This module's name. @@ -523,6 +527,19 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS) enddo endif + if (present(HA_CS)) then + call HA_init(Time, US, param_file, CS%time_ref, CS%nc, CS%freq, CS%phase0, CS%const_name, HA_CS) + call get_param(param_file, mdl, "HA_SSH", HA_ssh, & + "If true, perform harmonic analysis of sea serface height.", default=.false.) + if (HA_ssh) call HA_register('ssh', 'h', HA_CS) + call get_param(param_file, mdl, "HA_UBT", HA_ubt, & + "If true, perform harmonic analysis of zonal barotropic velocity.", default=.false.) + if (HA_ubt) call HA_register('ubt', 'u', HA_CS) + call get_param(param_file, mdl, "HA_VBT", HA_vbt, & + "If true, perform harmonic analysis of meridional barotropic velocity.", default=.false.) + if (HA_vbt) call HA_register('vbt', 'v', HA_CS) + endif + id_clock_tides = cpu_clock_id('(Ocean tides)', grain=CLOCK_MODULE) end subroutine tidal_forcing_init @@ -727,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. !! @@ -745,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 Arbic et al., 2004, DSR II). With -!! TIDAL_SAL_FROM_FILE or USE_PREVIOUS_TIDES, a list of input files -!! must be provided to describe each constituent's properties from +!! 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 !! 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 diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 561ace60a7..792dae1411 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -4256,8 +4256,8 @@ 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, +!! of various people, as described in the review paper by \cite niiler1977, +!! 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 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