From 53e936153c717c63e56a47b92f75b980537080e7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 14 May 2023 07:16:04 -0400 Subject: [PATCH 01/15] +*Redefine GV%Angstrom_H in non-Boussinesq mode Redefined GV%Angstrom_H in non-Boussinesq mode so that it is equal to GV%H_to_Z*GV%Angstrom_Z, just as it is in Boussinesq mode. This will change answers (slightly) in all cases with BOUSSINESQ = False. In addition, this commit adds the elements semi_Boussinesq, dZ_subroundoff, m2_s_to_HZ_T, HZ_T_to_m2_s and HZ_T_to_MKS to the verticalGrid_type. The first 3 new elements are used in rescaling vertical viscosities and diffusivities. The last two elements are set using the new runtime parameters SEMI_BOUSSINESQ and RHO_KV_CONVERT, which are only used or logged when BOUSSINESQ = False. All answers and output are identical in Boussinesq cases, but answers change and there are new runtime parameters in non-Boussinesq cases. --- src/core/MOM_verticalGrid.F90 | 55 +++++++++++++++++++++++++++++++---- 1 file changed, 49 insertions(+), 6 deletions(-) diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index d6003ca626..5e9b5c476c 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -41,12 +41,18 @@ module MOM_verticalGrid ! The following variables give information about the vertical grid. logical :: Boussinesq !< If true, make the Boussinesq approximation. + logical :: semi_Boussinesq !< If true, do non-Boussinesq pressure force calculations and + !! use mass-based "thicknesses, but use Rho0 to convert layer thicknesses + !! into certain height changes. This only applies if BOUSSINESQ is false. real :: Angstrom_H !< A one-Angstrom thickness in the model thickness units [H ~> m or kg m-2]. real :: Angstrom_Z !< A one-Angstrom thickness in the model depth units [Z ~> m]. real :: Angstrom_m !< A one-Angstrom thickness [m]. real :: H_subroundoff !< A thickness that is so small that it can be added to a thickness of !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2]. !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m. + real :: dZ_subroundoff !< A thickness in height units that is so small that it can be added to a + !! vertical distance of Angstrom_Z or 1e-17 m without changing it at the bit + !! level [Z ~> m]. This is the height equivalent of H_subroundoff. real, allocatable, dimension(:) :: & g_prime, & !< The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2]. Rlay !< The target coordinate value (potential density) in each layer [R ~> kg m-3]. @@ -74,8 +80,17 @@ module MOM_verticalGrid !! thickness units [H R-1 Z-1 ~> m3 kg-2 or 1]. real :: H_to_MKS !< A constant that translates thickness units to its MKS unit !! (m or kg m-2) based on GV%Boussinesq [m H-1 ~> 1] or [kg m-2 H-1 ~> 1] + real :: m2_s_to_HZ_T !< The combination of conversion factors that converts kinematic viscosities + !! in m2 s-1 to the internal units of the kinematic (in Boussinesq mode) + !! or dynamic viscosity [H Z s T-1 m-2 ~> 1 or kg m-3] + real :: HZ_T_to_m2_s !< The combination of conversion factors that converts the viscosities from + !! their internal representation into a kinematic viscosity in m2 s-1 + !! [T m2 H-1 Z-1 s-1 ~> 1 or m3 kg-1] + real :: HZ_T_to_MKS !< The combination of conversion factors that converts the viscosities from + !! their internal representation into their unnscaled MKS units + !! (m2 s-1 or Pa s), depending on whether the model is Boussinesq + !! [T m2 H-1 Z-1 s-1 ~> 1] or [T Pa s H-1 Z-1 ~> 1] - real :: m_to_H_restart = 1.0 !< A copy of the m_to_H that is used in restart files. end type verticalGrid_type contains @@ -91,6 +106,8 @@ subroutine verticalGridInit( param_file, GV, US ) ! Local variables integer :: nk, H_power real :: H_rescale_factor ! The integer power of 2 by which thicknesses are rescaled [nondim] + real :: rho_Kv ! The density used convert input kinematic viscosities into dynamic viscosities + ! when in non-Boussinesq mode [R ~> kg m-3] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=16) :: mdl = 'MOM_verticalGrid' @@ -114,6 +131,17 @@ subroutine verticalGridInit( param_file, GV, US ) units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "BOUSSINESQ", GV%Boussinesq, & "If true, make the Boussinesq approximation.", default=.true.) + call get_param(param_file, mdl, "SEMI_BOUSSINESQ", GV%semi_Boussinesq, & + "If true, do non-Boussinesq pressure force calculations and use mass-based "//& + "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//& + "height changes. This only applies if BOUSSINESQ is false.", & + default=.true., do_not_log=GV%Boussinesq) + if (GV%Boussinesq) GV%semi_Boussinesq = .true. + call get_param(param_file, mdl, "RHO_KV_CONVERT", Rho_Kv, & + "The density used to convert input kinematic viscosities into dynamic "//& + "viscosities in non-BOUSSINESQ mode, and similarly for vertical diffusivities.", & + units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R, & + do_not_log=GV%Boussinesq) call get_param(param_file, mdl, "ANGSTROM", GV%Angstrom_Z, & "The minimum layer thickness, usually one-Angstrom.", & units="m", default=1.0e-10, scale=US%m_to_Z) @@ -156,26 +184,41 @@ subroutine verticalGridInit( param_file, GV, US ) GV%H_to_kg_m2 = US%R_to_kg_m3*GV%Rho0 * GV%H_to_m GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2 GV%m_to_H = 1.0 / GV%H_to_m - GV%Angstrom_H = GV%m_to_H * US%Z_to_m*GV%Angstrom_Z GV%H_to_MKS = GV%H_to_m + GV%m2_s_to_HZ_T = GV%m_to_H * US%m_to_Z * US%T_to_s else GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2 GV%m_to_H = US%R_to_kg_m3*GV%Rho0 * GV%kg_m2_to_H GV%H_to_m = GV%H_to_kg_m2 / (US%R_to_kg_m3*GV%Rho0) - GV%Angstrom_H = US%Z_to_m*GV%Angstrom_Z * 1000.0*GV%kg_m2_to_H GV%H_to_MKS = GV%H_to_kg_m2 + GV%m2_s_to_HZ_T = US%R_to_kg_m3*rho_Kv * GV%kg_m2_to_H * US%m_to_Z * US%T_to_s endif - GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H,GV%m_to_H*1e-17) - GV%H_to_Pa = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth * GV%H_to_kg_m2 GV%H_to_Z = GV%H_to_m * US%m_to_Z GV%Z_to_H = US%Z_to_m * GV%m_to_H + + GV%Angstrom_H = GV%Z_to_H * GV%Angstrom_Z GV%Angstrom_m = US%Z_to_m * GV%Angstrom_Z + GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H, GV%m_to_H*1e-17) + GV%dZ_subroundoff = 1e-20 * max(GV%Angstrom_Z, US%m_to_Z*1e-17) + + GV%H_to_Pa = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth * GV%H_to_kg_m2 + GV%H_to_RZ = GV%H_to_kg_m2 * US%kg_m3_to_R * US%m_to_Z GV%RZ_to_H = GV%kg_m2_to_H * US%R_to_kg_m3 * US%Z_to_m -! Log derivative values. + GV%HZ_T_to_m2_s = 1.0 / GV%m2_s_to_HZ_T + GV%HZ_T_to_MKS = GV%H_to_MKS * US%Z_to_m * US%s_to_T + + ! Note based on the above that for both Boussinsq and non-Boussinesq cases that: + ! GV%Rho0 = GV%Z_to_H * GV%H_to_RZ + ! 1.0/GV%Rho0 = GV%H_to_Z * GV%RZ_to_H + ! This is exact for power-of-2 scaling of the units, regardless of the value of Rho0, but + ! the first term on the right hand side is invertable in Boussinesq mode, but the second + ! is invertable when non-Boussinesq. + + ! Log derivative values. call log_param(param_file, mdl, "M to THICKNESS", GV%m_to_H*H_rescale_factor, units="H m-1") call log_param(param_file, mdl, "M to THICKNESS rescaled by 2^-n", GV%m_to_H, units="2^n H m-1") call log_param(param_file, mdl, "THICKNESS to M rescaled by 2^n", GV%H_to_m, units="2^-n m H-1") From 1faa9ab08ff5bb6594b0c719afc7e3b56eab7966 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 17 May 2023 14:48:51 -0400 Subject: [PATCH 02/15] +Set_interp_answer_date and REGRIDDING_ANSWER_DATE Add the ability to set the answer date for the regridding code, including the addition of the new subroutine set_interp_answer_date and the new runtime parameter REGRIDDING_ANSWER_DATE to specify the code vintage to use with state- dependent vertical coordinates. There is also new optional argument to set_regrid_params. By default, all answers are bitwise identical, but there are new or modified public interfaces and there is a new entry in some MOM_parameter_doc files. --- src/ALE/MOM_regridding.F90 | 14 ++++++++++++-- src/ALE/regrid_interp.F90 | 17 ++++++++++++----- 2 files changed, 24 insertions(+), 7 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 74b7bc784a..9da4e95b24 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -23,7 +23,7 @@ module MOM_regridding use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA use regrid_consts, only : REGRIDDING_ARBITRARY, REGRIDDING_SIGMA_SHELF_ZSTAR use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_ADAPTIVE -use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap +use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap, set_interp_answer_date use coord_zlike, only : init_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column, end_coord_zlike use coord_sigma, only : init_coord_sigma, sigma_CS, set_sigma_params, build_sigma_column, end_coord_sigma @@ -212,6 +212,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. logical :: remap_answers_2018 integer :: remap_answer_date ! The vintage of the remapping expressions to use. + integer :: regrid_answer_date ! The vintage of the regridding expressions to use. real :: tmpReal, P_Ref real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z). real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha @@ -291,6 +292,13 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m "If both REMAPPING_2018_ANSWERS and REMAPPING_ANSWER_DATE are specified, the "//& "latter takes precedence.", default=default_answer_date) call set_regrid_params(CS, remap_answer_date=remap_answer_date) + call get_param(param_file, mdl, "REGRIDDING_ANSWER_DATE", regrid_answer_date, & + "The vintage of the expressions and order of arithmetic to use for regridding. "//& + "Values below 20190101 result in the use of older, less accurate expressions "//& + "that were in use at the end of 2018. Higher values result in the use of more "//& + "robust and accurate forms of mathematically equivalent expressions.", & + default=20181231) ! ### change to default=default_answer_date) + call set_regrid_params(CS, regrid_answer_date=regrid_answer_date) endif if (main_parameters .and. coord_is_state_dependent) then @@ -2233,7 +2241,7 @@ end function getCoordinateShortName subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, & interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, & compress_fraction, ref_pressure, & - integrate_downward_for_e, remap_answers_2018, remap_answer_date, & + integrate_downward_for_e, remap_answers_2018, remap_answer_date, regrid_answer_date, & adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin, adaptDrho0) type(regridding_CS), intent(inout) :: CS !< Regridding control structure logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells @@ -2252,6 +2260,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri !! that recover the remapping answers from 2018. Otherwise !! use more robust but mathematically equivalent expressions. integer, optional, intent(in) :: remap_answer_date !< The vintage of the expressions to use for remapping + integer, optional, intent(in) :: regrid_answer_date !< The vintage of the expressions to use for regridding real, optional, intent(in) :: adaptTimeRatio !< Ratio of the ALE timestep to the grid timescale [nondim]. real, optional, intent(in) :: adaptZoom !< Depth of near-surface zooming region [H ~> m or kg m-2]. real, optional, intent(in) :: adaptZoomCoeff !< Coefficient of near-surface zooming diffusivity [nondim]. @@ -2265,6 +2274,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri if (present(interp_scheme)) call set_interp_scheme(CS%interp_CS, interp_scheme) if (present(boundary_extrapolation)) call set_interp_extrap(CS%interp_CS, boundary_extrapolation) + if (present(regrid_answer_date)) call set_interp_answer_date(CS%interp_CS, regrid_answer_date) if (present(old_grid_weight)) then if (old_grid_weight<0. .or. old_grid_weight>1.) & diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index e119ce9d53..641ae7e6c2 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -33,14 +33,12 @@ module regrid_interp !! boundary cells logical :: boundary_extrapolation - !> The vintage of the expressions to use for remapping - integer :: answer_date = 20181231 - !### Changing this to 99991231 changes answers in rho and Hycom1 configurations. - !### There is no point where the value of answer_date is reset. + !> The vintage of the expressions to use for regridding + integer :: answer_date = 99991231 end type interp_CS_type public regridding_set_ppolys, build_and_interpolate_grid -public set_interp_scheme, set_interp_extrap +public set_interp_scheme, set_interp_extrap, set_interp_answer_date ! List of interpolation schemes integer, parameter :: INTERPOLATION_P1M_H2 = 0 !< O(h^2) @@ -547,4 +545,13 @@ subroutine set_interp_extrap(CS, extrap) CS%boundary_extrapolation = extrap end subroutine set_interp_extrap +!> Store the value of the answer_date in the interp_CS +subroutine set_interp_answer_date(CS, answer_date) + type(interp_CS_type), intent(inout) :: CS !< A control structure for regrid_interp + integer, intent(in) :: answer_date !< An integer encoding the vintage of + !! the expressions to use for regridding + + CS%answer_date = answer_date +end subroutine set_interp_answer_date + end module regrid_interp From e672b981b08c2e7af8f165122ba40bc10f4b4d19 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 25 May 2023 20:10:01 -0400 Subject: [PATCH 03/15] *+Revise non-Boussinesq find_coupling_coef calcs Restructure one of the find_coupling_coef calculations to draw out the stress-magnitude terms, in preparation for future steps to reduce the dependency on the Boussinesq reference density. Using a value of VERT_FRICTION_ANSWER_DATE that is below 20230601 recovers the previous answers with non-Boussinesq test cases, but this is irrelevant for Boussinesq test cases. This updated code is mathematically equivalent to the previous expressions but it does change answers at roundoff in non-Boussinesq cases for recent answer dates. There are modifications to some comments in MOM_parameter_doc files. --- .../vertical/MOM_vert_friction.F90 | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index ea6c7f112b..80fff62f21 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -139,8 +139,11 @@ module MOM_vert_friction integer :: answer_date !< The vintage of the order of arithmetic and expressions in the viscous !! calculations. Values below 20190101 recover the answers from the end !! of 2018, while higher values use expressions that do not use an - !! arbitrary and hard-coded maximum viscous coupling coefficient - !! between layers. + !! arbitrary and hard-coded maximum viscous coupling coefficient between + !! layers. In non-Boussinesq cases, values below 20230601 recover a + !! form of the viscosity within the mixed layer that breaks up the + !! magnitude of the wind stress with BULKMIXEDLAYER, DYNAMIC_VISCOUS_ML + !! or FIXED_DEPTH_LOTW_ML, but not LOTW_VISCOUS_ML_FLOOR. logical :: debug !< If true, write verbose checksums for debugging purposes. integer :: nkml !< The number of layers in the mixed layer. integer, pointer :: ntrunc !< The number of times the velocity has been @@ -1516,6 +1519,8 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real, dimension(SZIB_(G)) :: & u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1]. + tau_mag, & ! The magnitude of the wind stress at a velocity point including gustiness, + ! divided by the Boussinesq refernce density [Z2 T-2 ~> m2 s-2] absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1]. ! h_ml, & ! The mixed layer depth [H ~> m or kg m-2]. z_t, & ! The distance from the top, sometimes normalized @@ -1888,7 +1893,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*GV%H_to_Z ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer) ! and be further limited by rotation to give the natural Ekman length. - visc_ml = u_star(i) * CS%vonKar * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) + if (GV%Boussinesq .or. (CS%answer_date < 20230601)) then + visc_ml = u_star(i) * CS%vonKar * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) + else + tau_mag(i) = u_star(i)**2 + visc_ml = CS%vonKar * (temp1*tau_mag(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) + endif a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * GV%H_to_Z + 0.5*I_amax*visc_ml) ! Choose the largest estimate of a_cpl, but these could be changed to be additive. @@ -2180,7 +2190,9 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "The vintage of the order of arithmetic and expressions in the viscous "//& "calculations. Values below 20190101 recover the answers from the end of 2018, "//& "while higher values use expressions that do not use an arbitrary hard-coded "//& - "maximum viscous coupling coefficient between layers. "//& + "maximum viscous coupling coefficient between layers. Values below 20230601 "//& + "recover a form of the viscosity within the mixed layer that breaks up the "//& + "magnitude of the wind stress in some non-Boussinesq cases. "//& "If both VERT_FRICTION_2018_ANSWERS and VERT_FRICTION_ANSWER_DATE are "//& "specified, the latter takes precedence.", default=default_answer_date) From edb22ec5825c315b5498d3b5064b1af3175dd714 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 29 Apr 2023 11:00:54 -0400 Subject: [PATCH 04/15] +Code to calculate layer averaged specific volumes Add routines to calculate and store the layer-averaged specific volume, along with code to do the unit testing of this new capability. The new public interfaces include avg_specific_vol, average_specific_vol, avg_spec_vol_Wright, avg_spec_vol_Wright_full, avg_spec_vol_Wright_red and avg_spec_vol_linear. There is also a new optional argument to test_EOS_consistency to control whether these new capabilties are tested for a particular equation of state. All answers are bitwise identical, and the new capabilities pass the unit testing for self consistency. --- src/core/MOM_density_integrals.F90 | 34 ++++- src/equation_of_state/MOM_EOS.F90 | 144 ++++++++++++++++-- src/equation_of_state/MOM_EOS_Wright.F90 | 39 ++++- src/equation_of_state/MOM_EOS_Wright_full.F90 | 38 +++++ src/equation_of_state/MOM_EOS_Wright_red.F90 | 38 +++++ src/equation_of_state/MOM_EOS_linear.F90 | 26 +++- 6 files changed, 307 insertions(+), 12 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index e1fb3d3278..9fed528e71 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -4,12 +4,13 @@ module MOM_density_integrals ! This file is part of MOM6. See LICENSE.md for the license. use MOM_EOS, only : EOS_type -use MOM_EOS, only : EOS_quadrature +use MOM_EOS, only : EOS_quadrature, EOS_domain use MOM_EOS, only : analytic_int_density_dz use MOM_EOS, only : analytic_int_specific_vol_dp use MOM_EOS, only : calculate_density use MOM_EOS, only : calculate_spec_vol use MOM_EOS, only : calculate_specific_vol_derivs +use MOM_EOS, only : average_specific_vol use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_hor_index, only : hor_index_type use MOM_string_functions, only : uppercase @@ -28,6 +29,7 @@ module MOM_density_integrals public int_specific_vol_dp public int_spec_vol_dp_generic_pcm public int_spec_vol_dp_generic_plm +public avg_specific_vol public find_depth_of_pressure_in_cell contains @@ -1613,6 +1615,36 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t end subroutine find_depth_of_pressure_in_cell +!> Calculate the average in situ specific volume across layers +subroutine avg_specific_vol(T, S, p_t, dp, HI, EOS, SpV_avg, halo_size) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: T !< Potential temperature of the layer [C ~> degC] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: S !< Salinity of the layer [S ~> ppt] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa] + type(EOS_type), intent(in) :: EOS !< Equation of state structure + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [R-1 ~> m3 kg-1] + integer, optional, intent(in) :: halo_size !< The number of halo points in which to work. + + ! Local variables + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + integer :: jsh, jeh, j, halo + + halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0) + jsh = HI%jsc-halo ; jeh = HI%jec+halo + + EOSdom(:) = EOS_domain(HI, halo_size) + do j=jsh,jeh + call average_specific_vol(T(:,j), S(:,j), p_t(:,j), dp(:,j), SpV_avg(:,j), EOS, EOSdom) + enddo + +end subroutine avg_specific_vol !> Returns change in anomalous pressure change from top to non-dimensional !! position pos between z_t and z_b [R L2 T-2 ~> Pa] diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 276c4c3019..c68dc7b661 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -8,24 +8,25 @@ module MOM_EOS use MOM_EOS_linear, only : calculate_specvol_derivs_linear, int_density_dz_linear use MOM_EOS_linear, only : calculate_density_second_derivs_linear, EoS_fit_range_linear use MOM_EOS_linear, only : calculate_compress_linear, int_spec_vol_dp_linear +use MOM_EOS_linear, only : avg_spec_vol_linear use MOM_EOS_Wright, only : calculate_density_wright, calculate_spec_vol_wright use MOM_EOS_Wright, only : calculate_density_derivs_wright use MOM_EOS_Wright, only : calculate_specvol_derivs_wright, int_density_dz_wright use MOM_EOS_Wright, only : calculate_compress_wright, int_spec_vol_dp_wright use MOM_EOS_Wright, only : calculate_density_second_derivs_wright, calc_density_second_derivs_wright_buggy -use MOM_EOS_Wright, only : EoS_fit_range_Wright +use MOM_EOS_Wright, only : EoS_fit_range_Wright, avg_spec_vol_Wright use MOM_EOS_Wright_full, only : calculate_density_wright_full, calculate_spec_vol_wright_full use MOM_EOS_Wright_full, only : calculate_density_derivs_wright_full use MOM_EOS_Wright_full, only : calculate_specvol_derivs_wright_full, int_density_dz_wright_full use MOM_EOS_Wright_full, only : calculate_compress_wright_full, int_spec_vol_dp_wright_full use MOM_EOS_Wright_full, only : calculate_density_second_derivs_wright_full -use MOM_EOS_Wright_full, only : EoS_fit_range_Wright_full +use MOM_EOS_Wright_full, only : EoS_fit_range_Wright_full, avg_spec_vol_Wright_full use MOM_EOS_Wright_red, only : calculate_density_wright_red, calculate_spec_vol_wright_red use MOM_EOS_Wright_red, only : calculate_density_derivs_wright_red use MOM_EOS_Wright_red, only : calculate_specvol_derivs_wright_red, int_density_dz_wright_red use MOM_EOS_Wright_red, only : calculate_compress_wright_red, int_spec_vol_dp_wright_red use MOM_EOS_Wright_red, only : calculate_density_second_derivs_wright_red -use MOM_EOS_Wright_red, only : EoS_fit_range_Wright_red +use MOM_EOS_Wright_red, only : EoS_fit_range_Wright_red, avg_spec_vol_Wright_red use MOM_EOS_Jackett06, only : calculate_density_Jackett06, calculate_spec_vol_Jackett06 use MOM_EOS_Jackett06, only : calculate_density_derivs_Jackett06, calculate_specvol_derivs_Jackett06 use MOM_EOS_Jackett06, only : calculate_compress_Jackett06, calculate_density_second_derivs_Jackett06 @@ -68,6 +69,7 @@ module MOM_EOS public EOS_unit_tests public analytic_int_density_dz public analytic_int_specific_vol_dp +public average_specific_vol public calculate_compress public calculate_density public calculate_density_derivs @@ -1324,6 +1326,97 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) end subroutine calculate_compress_scalar +!> Calls the appropriate subroutine to calculate the layer averaged specific volume either using +!! Boole's rule quadrature or analytical and nearly-analytical averages in pressure. +subroutine average_specific_vol(T, S, p_t, dp, SpV_avg, EOS, dom, scale) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] + real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt] + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [R-1 ~> m3 kg-1] + type(EOS_type), intent(in) :: EOS !< Equation of state structure + integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking + !! into account that arrays start at 1. + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale + !! output specific volume in combination with + !! scaling stored in EOS [various] + + ! Local variables + real, dimension(size(T)) :: pres ! Layer-top pressure converted to [Pa] + real, dimension(size(T)) :: dpres ! Pressure change converted to [Pa] + real, dimension(size(T)) :: Ta ! Temperature converted to [degC] + real, dimension(size(T)) :: Sa ! Salinity converted to [ppt] + real :: T5(5) ! Temperatures at five quadrature points [C ~> degC] + real :: S5(5) ! Salinities at five quadrature points [S ~> ppt] + real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa] + real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] + real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] + real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] + integer :: i, n, is, ie, npts + + if (present(dom)) then + is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is + else + is = 1 ; ie = size(T) ; npts = 1 + ie - is + endif + + if (EOS%EOS_quadrature) then + do i=is,ie + do n=1,5 + T5(n) = T(i) ; S5(n) = S(i) + p5(n) = p_t(i) + 0.25*real(5-n)*dp(i) + enddo + call calculate_spec_vol(T5, S5, p5, a5, EOS) + + ! Use Boole's rule to estimate the average specific volume. + SpV_avg(i) = C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3)) + enddo + elseif ((EOS%RL2_T2_to_Pa == 1.0) .and. (EOS%C_to_degC == 1.0) .and. (EOS%S_to_ppt == 1.0)) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call avg_spec_vol_linear(T, S, p_t, dp, SpV_avg, is, npts, EOS%Rho_T0_S0, & + EOS%dRho_dT, EOS%dRho_dS) + case (EOS_WRIGHT) + call avg_spec_vol_wright(T, S, p_t, dp, SpV_avg, is, npts) + case (EOS_WRIGHT_FULL) + call avg_spec_vol_wright_full(T, S, p_t, dp, SpV_avg, is, npts) + case (EOS_WRIGHT_REDUCED) + call avg_spec_vol_wright_red(T, S, p_t, dp, SpV_avg, is, npts) + case default + call MOM_error(FATAL, "No analytic average specific volume option is available with this EOS!") + end select + else + do i=is,ie + pres(i) = EOS%RL2_T2_to_Pa * p_t(i) + dpres(i) = EOS%RL2_T2_to_Pa * dp(i) + Ta(i) = EOS%C_to_degC * T(i) + Sa(i) = EOS%S_to_ppt * S(i) + enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call avg_spec_vol_linear(Ta, Sa, pres, dpres, SpV_avg, is, npts, EOS%Rho_T0_S0, & + EOS%dRho_dT, EOS%dRho_dS) + case (EOS_WRIGHT) + call avg_spec_vol_wright(Ta, Sa, pres, dpres, SpV_avg, is, npts) + case (EOS_WRIGHT_FULL) + call avg_spec_vol_wright_full(Ta, Sa, pres, dpres, SpV_avg, is, npts) + case (EOS_WRIGHT_REDUCED) + call avg_spec_vol_wright_red(Ta, Sa, pres, dpres, SpV_avg, is, npts) + case default + call MOM_error(FATAL, "No analytic average specific volume option is available with this EOS!") + end select + endif + + spv_scale = EOS%R_to_kg_m3 + if (EOS%EOS_quadrature) spv_scale = 1.0 + if (present(scale)) spv_scale = spv_scale * scale + if (spv_scale /= 1.0) then ; do i=is,ie + SpV_avg(i) = spv_scale * SpV_avg(i) + enddo ; endif + +end subroutine average_specific_vol + !> Return the range of temperatures, salinities and pressures for which the equation of state that !! is being used has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -2057,13 +2150,13 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_FULL) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_FULL", & - rho_check=1027.55177447616*EOS_tmp%kg_m3_to_R) + rho_check=1027.55177447616*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_FULL EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_REDUCED) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_REDUCED", & - rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R) + rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_REDUCED EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail @@ -2076,7 +2169,7 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT", & - rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R) + rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail @@ -2126,7 +2219,7 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_LINEAR, Rho_T0_S0=1000.0, drho_dT=-0.2, dRho_dS=0.8) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "LINEAR", & - rho_check=1023.0*EOS_tmp%kg_m3_to_R) + rho_check=1023.0*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "LINEAR EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail @@ -2293,7 +2386,7 @@ end subroutine write_check_msg !> Test an equation of state for self-consistency and consistency with check values, returning false !! if it is consistent by all tests, and true if it fails any test. logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & - EOS_name, rho_check, spv_check, skip_2nd) result(inconsistent) + EOS_name, rho_check, spv_check, skip_2nd, avg_Sv_check) result(inconsistent) real, intent(in) :: T_test !< Potential temperature or conservative temperature [C ~> degC] real, intent(in) :: S_test !< Salinity or absolute salinity [S ~> ppt] real, intent(in) :: p_test !< Pressure [R L2 T-2 ~> Pa] @@ -2302,7 +2395,9 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & character(len=*), intent(in) :: EOS_name !< A name used in error messages to describe the EoS real, optional, intent(in) :: rho_check !< A check value for the density [R ~> kg m-3] real, optional, intent(in) :: spv_check !< A check value for the specific volume [R-1 ~> m3 kg-1] - logical, optional, intent(in) :: skip_2nd !< If present and true, do not check the 2nd derivatives. + logical, optional, intent(in) :: skip_2nd !< If present and true, do not check the 2nd derivatives. + logical, optional, intent(in) :: avg_Sv_check !< If present and true, compare analytical and numerical + !! quadrature estimates of the layer-averaged specific volume. ! Local variables real, dimension(-3:3,-3:3,-3:3) :: T ! Temperatures at the test value and perturbed points [C ~> degC] @@ -2329,6 +2424,8 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & ! temperature [R-1 C-1 ~> m3 kg-1 degC-1] real :: dSV_dS(1) ! The partial derivative of specific volume with salinity ! [R-1 S-1 ~> m3 kg-1 ppt-1] + real :: SpV_avg_a(1) ! The pressure-averaged specific volume determined analytically [R-1 ~> m3 kg-1] + real :: SpV_avg_q(1) ! The pressure-averaged specific volume determined via quadrature [R-1 ~> m3 kg-1] real :: drho_dS_dS ! Second derivative of density with respect to S [R S-2 ~> kg m-3 ppt-2] real :: drho_dS_dT ! Second derivative of density with respect to T and S [R S-1 C-1 ~> kg m-3 ppt-1 degC-1] real :: drho_dT_dT ! Second derivative of density with respect to T [R C-2 ~> kg m-3 degC-2] @@ -2370,13 +2467,17 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & real :: count_fac2 ! A factor in the roundoff estimates based on the factors in the numerator and ! denominator in the finite difference second derivative expression [nondim] character(len=200) :: mesg + type(EOS_type) :: EOS_tmp logical :: test_OK ! True if a particular test is consistent. logical :: OK ! True if all checks so far are consistent. logical :: test_2nd ! If true, do tests on the 2nd derivative calculations + logical :: test_avg_Sv ! If true, compare numerical and analytical estimates of the vertically + ! averaged specific volume integer :: order ! The order of accuracy of the centered finite difference estimates (2, 4 or 6). integer :: i, j, k, n test_2nd = .true. ; if (present(skip_2nd)) test_2nd = .not.skip_2nd + test_avg_Sv = .false. ; if (present(avg_Sv_check)) test_avg_Sv = avg_Sv_check dT = 0.1*EOS%degC_to_C ! Temperature perturbations [C ~> degC] dS = 0.5*EOS%ppt_to_S ! Salinity perturbations [S ~> ppt] @@ -2442,6 +2543,14 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, EOS) call calculate_compress(T(0,0,0), S(0,0,0), p(0,0,0), rho_tmp, drho_dp, EOS) + if (test_avg_Sv) then + EOS_tmp = EOS + call EOS_manual_init(EOS_tmp, EOS_quadrature=.false.) + call average_specific_vol(T(0:0,0,0), S(0:0,0,0), p(0:0,0,0), p(0:0,0,0), SpV_avg_a, EOS_tmp) + call EOS_manual_init(EOS_tmp, EOS_quadrature=.true.) + call average_specific_vol(T(0:0,0,0), S(0:0,0,0), p(0:0,0,0), p(0:0,0,0), SpV_avg_q, EOS_tmp) + endif + OK = .true. tol = 1000.0*epsilon(tol) @@ -2532,6 +2641,23 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & OK = OK .and. check_FD(drho_dS_dP, drho_dS_dP_fd, tol_here, verbose, trim(EOS_name)//" drho_dS_dP", order) endif + if (test_avg_Sv) then + tol_here = 0.5*tol*(abs(SpV_avg_a(1)) + abs(SpV_avg_q(1))) + test_OK = (abs(SpV_avg_a(1) - SpV_avg_q(1)) < tol_here) + if (verbose) then + write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2"), tol=",ES16.8)') & + SpV_avg_a(1), SpV_avg_q(1), SpV_avg_a(1) - SpV_avg_q(1), & + 2.0*(SpV_avg_a(1) - SpV_avg_q(1)) / (abs(SpV_avg_a(1)) + abs(SpV_avg_q(1)) + tiny(SpV_avg_a(1))), & + tol_here + if (verbose .and. .not.test_OK) then + call MOM_error(WARNING, "The values of "//trim(EOS_name)//" SpV_avg disagree. "//trim(mesg)) + elseif (verbose) then + call MOM_mesg("The values of "//trim(EOS_name)//" SpV_avg agree: "//trim(mesg)) + endif + endif + OK = OK .and. test_OK + endif + inconsistent = .not.OK contains diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 25ae9219a8..d8dee28aa2 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -10,7 +10,7 @@ module MOM_EOS_Wright public calculate_compress_wright, calculate_density_wright, calculate_spec_vol_wright public calculate_density_derivs_wright, calculate_specvol_derivs_wright public calculate_density_second_derivs_wright, calc_density_second_derivs_wright_buggy -public EoS_fit_range_Wright +public EoS_fit_range_Wright, avg_spec_vol_Wright public int_density_dz_wright, int_spec_vol_dp_wright !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to @@ -547,6 +547,42 @@ subroutine calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) enddo end subroutine calculate_compress_wright +!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine +!! the layer-average specific volumes. There are essentially no free assumptions, apart from a +!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. +subroutine avg_spec_vol_Wright(T, S, p_t, dp, SpV_avg, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! [degC]. + real, dimension(:), intent(in) :: S !< Salinity [PSU]. + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + + ! Local variables + real :: al0 ! A term in the Wright EOS [m3 kg-1] + real :: p0 ! A term in the Wright EOS [Pa] + real :: lambda ! A term in the Wright EOS [m2 s-2] + real :: eps2 ! The square of a nondimensional ratio [nondim] + real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1]. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim] + integer :: j + + ! alpha(j) = al0 + lambda / (pressure(j) + p0) + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + + I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j))) + eps2 = (0.5 * dp(j) * I_pterm)**2 + SpV_avg(j) = al0 + (lambda * I_pterm) * & + (1.0 + eps2*(C1_3 + eps2*(0.2 + eps2*(C1_7 + eps2*C1_9)))) + enddo +end subroutine avg_spec_vol_Wright + !> Return the range of temperatures, salinities and pressures for which the reduced-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -1066,6 +1102,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & enddo ; enddo ; endif end subroutine int_spec_vol_dp_wright + !> \namespace mom_eos_wright !! !! \section section_EOS_Wright Wright equation of state diff --git a/src/equation_of_state/MOM_EOS_Wright_full.F90 b/src/equation_of_state/MOM_EOS_Wright_full.F90 index 3f00a92cef..107ced3f5b 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -11,6 +11,7 @@ module MOM_EOS_Wright_full public calculate_density_derivs_wright_full, calculate_specvol_derivs_wright_full public calculate_density_second_derivs_wright_full, EoS_fit_range_Wright_full public int_density_dz_wright_full, int_spec_vol_dp_wright_full +public avg_spec_vol_Wright_full !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to !! a reference density, from salinity in practical salinity units ([PSU]), potential @@ -450,6 +451,42 @@ subroutine calculate_compress_wright_full(T, S, pressure, rho, drho_dp, start, n enddo end subroutine calculate_compress_wright_full +!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine +!! the layer-average specific volumes. There are essentially no free assumptions, apart from a +!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. +subroutine avg_spec_vol_Wright_full(T, S, p_t, dp, SpV_avg, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! [degC]. + real, dimension(:), intent(in) :: S !< Salinity [PSU]. + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + + ! Local variables + real :: al0 ! A term in the Wright EOS [m3 kg-1] + real :: p0 ! A term in the Wright EOS [Pa] + real :: lambda ! A term in the Wright EOS [m2 s-2] + real :: eps2 ! The square of a nondimensional ratio [nondim] + real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1]. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim] + integer :: j + + ! alpha(j) = al0 + lambda / (pressure(j) + p0) + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + + I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j))) + eps2 = (0.5 * dp(j) * I_pterm)**2 + SpV_avg(j) = al0 + (lambda * I_pterm) * & + (1.0 + eps2*(C1_3 + eps2*(0.2 + eps2*(C1_7 + eps2*C1_9)))) + enddo +end subroutine avg_spec_vol_Wright_full + !> Return the range of temperatures, salinities and pressures for which full-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -972,6 +1009,7 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & enddo ; enddo ; endif end subroutine int_spec_vol_dp_wright_full + !> \namespace mom_eos_wright_full !! !! \section section_EOS_Wright_full Wright equation of state diff --git a/src/equation_of_state/MOM_EOS_Wright_red.F90 b/src/equation_of_state/MOM_EOS_Wright_red.F90 index cf78ce2211..5553112274 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -11,6 +11,7 @@ module MOM_EOS_Wright_red public calculate_density_derivs_wright_red, calculate_specvol_derivs_wright_red public calculate_density_second_derivs_wright_red, EoS_fit_range_Wright_red public int_density_dz_wright_red, int_spec_vol_dp_wright_red +public avg_spec_vol_Wright_red !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to !! a reference density, from salinity in practical salinity units ([PSU]), potential @@ -450,6 +451,42 @@ subroutine calculate_compress_wright_red(T, S, pressure, rho, drho_dp, start, np enddo end subroutine calculate_compress_wright_red +!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine +!! the layer-average specific volumes. There are essentially no free assumptions, apart from a +!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. +subroutine avg_spec_vol_Wright_red(T, S, p_t, dp, SpV_avg, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! [degC]. + real, dimension(:), intent(in) :: S !< Salinity [PSU]. + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + + ! Local variables + real :: al0 ! A term in the Wright EOS [m3 kg-1] + real :: p0 ! A term in the Wright EOS [Pa] + real :: lambda ! A term in the Wright EOS [m2 s-2] + real :: eps2 ! The square of a nondimensional ratio [nondim] + real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1]. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim] + integer :: j + + ! alpha(j) = al0 + lambda / (pressure(j) + p0) + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + + I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j))) + eps2 = (0.5 * dp(j) * I_pterm)**2 + SpV_avg(j) = al0 + (lambda * I_pterm) * & + (1.0 + eps2*(C1_3 + eps2*(0.2 + eps2*(C1_7 + eps2*C1_9)))) + enddo +end subroutine avg_spec_vol_Wright_red + !> Return the range of temperatures, salinities and pressures for which the reduced-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -972,6 +1009,7 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & enddo ; enddo ; endif end subroutine int_spec_vol_dp_wright_red + !> \namespace mom_eos_wright_red !! !! \section section_EOS_Wright_red Wright equation of state diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index 1899103f5d..b1dacf2780 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -13,6 +13,7 @@ module MOM_EOS_linear public calculate_density_scalar_linear, calculate_density_array_linear public calculate_density_second_derivs_linear, EoS_fit_range_linear public int_density_dz_linear, int_spec_vol_dp_linear +public avg_spec_vol_linear ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -292,7 +293,7 @@ end subroutine calculate_specvol_derivs_linear !> This subroutine computes the in situ density of sea water (rho) !! and the compressibility (drho/dp == C_sound^-2) at the given !! salinity, potential temperature, and pressure. -subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,& +subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts, & Rho_T0_S0, dRho_dT, dRho_dS) real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface !! [degC]. @@ -318,6 +319,29 @@ subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,& enddo end subroutine calculate_compress_linear +!> Calculates the layer average specific volumes. +subroutine avg_spec_vol_linear(T, S, p_t, dp, SpV_avg, start, npts, Rho_T0_S0, dRho_dT, dRho_dS) + real, dimension(:), intent(in) :: T !< Potential temperature [degC] + real, dimension(:), intent(in) :: S !< Salinity [PSU] + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3] + real, intent(in) :: dRho_dT !< The derivative of density with temperature + !! [kg m-3 degC-1] + real, intent(in) :: dRho_dS !< The derivative of density with salinity + !! [kg m-3 ppt-1] + ! Local variables + integer :: j + + do j=start,start+npts-1 + SpV_avg(j) = 1.0 / (Rho_T0_S0 + (dRho_dT*T(j) + dRho_dS*S(j))) + enddo +end subroutine avg_spec_vol_linear + !> Return the range of temperatures, salinities and pressures for which the reduced-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. From e86b35adfd9bc3828aa9ae66d69d860e04358da7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 29 Apr 2023 11:04:47 -0400 Subject: [PATCH 05/15] +Add thickness_to_dz and calc_derived_thermo Added the new overloaded interface thickness_to_dz to convert the layer thicknesses in thickness units [H ~> m or kg m-2] into vertical distances in [Z ~> m], with variants that set full 3-d arrays or an i-/k- slice. Also added a field (SpV_avg) for the layer-averaged specific volume to the thermo_vars_ptr type and the new subroutine calc_derived_thermo to set it. This new subroutine is being called after halo updates to the temperatures and salinities. The new runtime parameter SEMI_BOUSSINESQ was added to determine whether tv%SpV_avg is allocated and used; it is stored in GV%semi_Boussinesq. Also added the new element GV%dZ_subroundoff to the verticalGrid_type as a counterpart to GV%H_subroundoff but in height units. All answers are bitwise identical, but there is a new runtime parameter in some MOM_parameter_doc files, new elements in a transparent type and a new public interface. --- src/core/MOM.F90 | 42 +++++- src/core/MOM_interface_heights.F90 | 124 +++++++++++++++++- src/core/MOM_variables.F90 | 3 + .../vertical/MOM_diabatic_driver.F90 | 10 +- src/tracer/MOM_offline_main.F90 | 8 ++ 5 files changed, 179 insertions(+), 8 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c61ed72e0c..af8481fd1c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -91,7 +91,7 @@ module MOM use MOM_grid, only : set_first_direction, rescale_grid_bathymetry 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 +use MOM_interface_heights, only : find_eta, calc_derived_thermo use MOM_interface_filter, only : interface_filter, interface_filter_init, interface_filter_end use MOM_interface_filter, only : interface_filter_CS use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end @@ -1400,6 +1400,12 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) call create_group_pass(pass_T_S, CS%tv%T, G%Domain, To_All+Omit_Corners, halo=1) call create_group_pass(pass_T_S, CS%tv%S, G%Domain, To_All+Omit_Corners, halo=1) call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass) + halo_sz = 1 + endif + + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz) endif endif @@ -1581,6 +1587,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call create_group_pass(pass_uv_T_S_h, h, G%Domain, halo=dynamics_stencil) call do_group_pass(pass_uv_T_S_h, G%Domain, clock=id_clock_pass) + ! Update derived thermodynamic quantities. + if (allocated(tv%SpV_avg)) then + call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil) + endif + if (CS%debug .and. CS%use_ALE_algorithm) then call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US) call hchksum(tv%T, "Post-ALE T", G%HI, haloshift=1, scale=US%C_to_degC) @@ -1623,13 +1634,19 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call cpu_clock_end(id_clock_adiabatic) if (associated(tv%T)) then - call create_group_pass(pass_T_S, tv%T, G%Domain, To_All+Omit_Corners, halo=1) - call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=1) + dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call create_group_pass(pass_T_S, tv%T, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil) + call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil) call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass) if (CS%debug) then if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, scale=US%C_to_degC) if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, scale=US%S_to_ppt) endif + + ! Update derived thermodynamic quantities. + if (allocated(tv%SpV_avg)) then + call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil) + endif endif endif ! endif for the block "if (.not.CS%adiabatic)" @@ -1676,6 +1693,8 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS type(time_type), pointer :: accumulated_time => NULL() type(time_type), pointer :: vertical_time => NULL() + integer :: dynamics_stencil ! The computational stencil for the calculations + ! in the dynamic core. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz ! 3D pointers @@ -1848,6 +1867,12 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS fluxes%fluxes_used = .true. + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil) + endif + if (last_iter) then accumulated_time = real_to_time(0.0) endif @@ -2817,6 +2842,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif endif + ! Allocate any derived equation of state fields. + if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then + allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0) + endif + if (use_ice_shelf .and. CS%debug) then call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0) call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0,scale=US%RZ_to_kg_m2) @@ -3103,6 +3133,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call do_group_pass(pass_uv_T_S_h, G%Domain) + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil) + endif + if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) @@ -3931,6 +3966,7 @@ subroutine MOM_end(CS) if (associated(CS%Hml)) deallocate(CS%Hml) if (associated(CS%tv%salt_deficit)) deallocate(CS%tv%salt_deficit) if (associated(CS%tv%frazil)) deallocate(CS%tv%frazil) + if (allocated(CS%tv%SpV_avg)) deallocate(CS%tv%SpV_avg) if (associated(CS%tv%T)) then DEALLOC_(CS%T) ; CS%tv%T => NULL() ; DEALLOC_(CS%S) ; CS%tv%S => NULL() diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 4f41cb074b..befeb1c2ad 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -3,7 +3,7 @@ module MOM_interface_heights ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_density_integrals, only : int_specific_vol_dp +use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol use MOM_error_handler, only : MOM_error, FATAL use MOM_EOS, only : calculate_density, EOS_type, EOS_domain use MOM_file_parser, only : log_version @@ -16,18 +16,26 @@ module MOM_interface_heights #include -public find_eta, dz_to_thickness, dz_to_thickness_simple +public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple +public calc_derived_thermo !> Calculates the heights of the free surface or all interfaces from layer thicknesses. interface find_eta module procedure find_eta_2d, find_eta_3d end interface find_eta -!> Calculates layer thickness in thickness units from geometric thicknesses in height units. +!> Calculates layer thickness in thickness units from geometric distance between the +!! interfaces around that layer in height units. interface dz_to_thickness module procedure dz_to_thickness_tv, dz_to_thickness_EoS end interface dz_to_thickness +!> Converts layer thickness in thickness units into the vertical distance between the +!! interfaces around a layer in height units. +interface thickness_to_dz + module procedure thickness_to_dz_3d, thickness_to_dz_jslice +end interface thickness_to_dz + contains !> Calculates the heights of all interfaces between layers, using the appropriate @@ -253,6 +261,45 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref) end subroutine find_eta_2d +!> Calculate derived thermodynamic quantities for re-use later. +subroutine calc_derived_thermo(tv, h, G, GV, US, halo) + type(ocean_grid_type), intent(in) :: 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 + type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various + !! thermodynamic variables, some of + !! which will be set here. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + integer, optional, intent(in) :: halo !< Width of halo within which to + !! calculate thicknesses + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa] + real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa] + integer :: i, j, k, is, ie, js, je, halos, nz + + halos = 0 ; if (present(halo)) halos = max(0,halo) + is = G%isc-halos ; ie = G%iec+halos ; js = G%jsc-halos ; je = G%jec+halos ; nz = GV%ke + + if (allocated(tv%Spv_avg) .and. associated(tv%eqn_of_state)) then + if (associated(tv%p_surf)) then + do j=js,je ; do i=is,ie ; p_t(i,j) = tv%p_surf(i,j) ; enddo ; enddo + else + do j=js,je ; do i=is,ie ; p_t(i,j) = 0.0 ; enddo ; enddo + endif + do k=1,nz + do j=js,je ; do i=is,ie + dp(i,j) = GV%g_Earth*GV%H_to_RZ*h(i,j,k) + enddo ; enddo + call avg_specific_vol(tv%T(:,:,k), tv%S(:,:,k), p_t, dp, G%HI, tv%eqn_of_state, tv%SpV_avg(:,:,k), halo) + if (k Converts thickness from geometric height units to thickness units, perhaps via an !! inversion of the integral of the density in pressure using variables stored in !! the thermo_var_ptrs type when in non-Boussinesq mode. @@ -428,4 +475,75 @@ subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode) end subroutine dz_to_thickness_simple +!> Converts layer thicknesses in thickness units to the vertical distance between edges in height +!! units, perhaps by multiplication by the precomputed layer-mean specific volume stored in an +!! array in the thermo_var_ptrs type when in non-Boussinesq mode. +subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size) + type(ocean_grid_type), intent(in) :: 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 + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2]. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m] + !! This is essentially intent out, but declared as intent + !! inout to preserve any initialized values in halo points. + integer, optional, intent(in) :: halo_size !< Width of halo within which to + !! calculate thicknesses + ! Local variables + integer :: i, j, k, is, ie, js, je, halo, nz + + halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) + is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke + + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do k=1,nz ; do j=js,je ; do i=is,ie + dz(i,j,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) + enddo ; enddo ; enddo + else + do k=1,nz ; do j=js,je ; do i=is,ie + dz(i,j,k) = GV%H_to_Z * h(i,j,k) + enddo ; enddo ; enddo + endif + +end subroutine thickness_to_dz_3d + + +!> Converts a vertical i- / k- slice of layer thicknesses in thickness units to the vertical +!! distance between edges in height units, perhaps by multiplication by the precomputed layer-mean +!! specific volume stored in an array in the thermo_var_ptrs type when in non-Boussinesq mode. +subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2]. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + real, dimension(SZI_(G),SZK_(GV)), & + intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m] + !! This is essentially intent out, but declared as intent + !! inout to preserve any initialized values in halo points. + integer, intent(in) :: j !< The second (j-) index of the input thicknesses to work with + integer, optional, intent(in) :: halo_size !< Width of halo within which to + !! calculate thicknesses + ! Local variables + integer :: i, k, is, ie, halo, nz + + halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) + is = G%isc-halo ; ie = G%iec+halo ; nz = GV%ke + + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do k=1,nz ; do i=is,ie + dz(i,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) + enddo ; enddo + else + do k=1,nz ; do i=is,ie + dz(i,k) = GV%H_to_Z * h(i,j,k) + enddo ; enddo + endif + +end subroutine thickness_to_dz_jslice + end module MOM_interface_heights diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 5efb02fe44..bec93376af 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -93,6 +93,9 @@ module MOM_variables logical :: S_is_absS = .false. !< If true, the salinity variable tv%S is !! actually the absolute salinity in units of [gSalt kg-1]. real :: min_salinity !< The minimum value of salinity when BOUND_SALINITY=True [S ~> ppt]. + real, allocatable, dimension(:,:,:) :: SpV_avg + !< The layer averaged in situ specific volume [R-1 ~> m3 kg-1]. + ! These arrays are accumulated fluxes for communication with other components. real, dimension(:,:), pointer :: frazil => NULL() !< The energy needed to heat the ocean column to the diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 59d5aaf60a..0fe08a06b2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -43,7 +43,7 @@ module MOM_diabatic_driver use MOM_grid, only : ocean_grid_type use MOM_int_tide_input, only : set_int_tide_input, int_tide_input_init use MOM_int_tide_input, only : int_tide_input_end, int_tide_input_CS, int_tide_input_type -use MOM_interface_heights, only : find_eta +use MOM_interface_heights, only : find_eta, calc_derived_thermo use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used @@ -1844,9 +1844,15 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Also changes: visc%Kd_shear and visc%Kv_shear if ((CS%halo_TS_diff > 0) .and. (CS%ML_mix_first > 0.0)) then if (associated(tv%T)) call pass_var(tv%T, G%Domain, halo=CS%halo_TS_diff, complete=.false.) - if (associated(tv%T)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.) + if (associated(tv%S)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.) call pass_var(h, G%domain, halo=CS%halo_TS_diff, complete=.true.) endif + + ! Update derived thermodynamic quantities. + if ((CS%ML_mix_first > 0.0) .and. allocated(tv%SpV_avg)) then + call calc_derived_thermo(tv, h, G, GV, US, halo=CS%halo_TS_diff) + endif + if (CS%debug) & call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff) if (CS%double_diffuse) then diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index ea6167a6b8..40dced9b20 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -22,6 +22,7 @@ module MOM_offline_main use MOM_file_parser, only : read_param, get_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : calc_derived_thermo use MOM_io, only : MOM_read_data, MOM_read_vector, CENTER use MOM_offline_aux, only : update_offline_from_arrays, update_offline_from_files use MOM_offline_aux, only : next_modulo_time, offline_add_diurnal_sw @@ -1025,6 +1026,7 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale) type(forcing), intent(inout) :: fluxes !< Pointers to forcing fields logical, intent(in ) :: do_ale !< True if using ALE ! Local variables + integer :: stencil integer :: i, j, k, is, ie, js, je, nz real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_start ! Initial thicknesses [H ~> m or kg m-2] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1086,6 +1088,12 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale) call pass_var(CS%tv%T, G%Domain) call pass_var(CS%tv%S, G%Domain) + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call calc_derived_thermo(CS%tv, CS%h_end, G, GV, US, halo=stencil) + endif + ! Update the read indices CS%ridx_snap = next_modulo_time(CS%ridx_snap,CS%numtime) CS%ridx_sum = next_modulo_time(CS%ridx_sum,CS%numtime) From e0021fc0679e2896d9466055ea87707fcb76d166 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Fri, 21 Apr 2023 09:35:06 -0400 Subject: [PATCH 06/15] wave structure computation into wave_speeds wave_speeds now computes the wave structures (eigenvectors) for each mode speed (eigenvalue) similarly to the wave_speed (singular) function. This is a replacement for the MOM_wave_structure function, which could be removed in a subsequent PR. Additional arrays for mode strucures and integral quantities are passed as output hence this is a breaking change for the call to wave_speeds. However it is only called once in diabatic_driver and is used exclusively for internal tides ray tracing. The dimensional solutions for the wave structures are now computed inside MOM_internal_tides, and new diagnostics are added. An out-of-bounds bug is also corrected for the computation of an averaged coriolis parameter. --- src/diagnostics/MOM_wave_speed.F90 | 257 ++++++++++++++++-- .../lateral/MOM_internal_tides.F90 | 153 ++++++++++- .../vertical/MOM_diabatic_driver.F90 | 5 +- 3 files changed, 384 insertions(+), 31 deletions(-) diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 9c8cd099f3..bb1b381c15 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -7,7 +7,7 @@ module MOM_wave_speed use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : log_version use MOM_grid, only : ocean_grid_type -use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h +use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h, interpolate_column use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -651,17 +651,33 @@ subroutine tdma6(n, a, c, lam, y) end subroutine tdma6 !> Calculates the wave speeds for the first few barolinic modes. -subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables - integer, intent(in) :: nmodes !< Number of modes - real, dimension(G%isd:G%ied,G%jsd:G%jed,nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1] - type(wave_speed_CS), intent(in) :: CS !< Wave speed control struct - logical, optional, intent(in) :: full_halos !< If true, do the calculation - !! over the entire data domain. +subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_struct_max, u_struct_bot, Nb, int_w2, & + int_U2, int_N2w2, full_halos) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + integer, intent(in) :: nmodes !< Number of modes + type(wave_speed_CS), intent(in) :: CS !< Wave speed control struct + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1,nmodes),intent(out) :: w_struct !< Wave Vertical profile [nondim] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV),nmodes),intent(out) :: u_struct !< Wave Horizontal profile [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_max !< Maximum of wave horizontal profile + !! [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_bot !< Bottom value of wave horizontal + !! profile [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Nb !< Bottom value of Brunt Vaissalla freqency + !! [T-1 ~> s-1] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_w2 !< depth-integrated + !! vertical profile squared [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_U2 !< depth-integrated + !! horizontal profile squared [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_N2w2 !< depth-integrated Brunt Vaissalla + !! frequency times vertical + !! profile squared [Z T-2 ~> m s-2] + logical, optional, intent(in) :: full_halos !< If true, do the calculation + !! over the entire data domain. ! Local variables real, dimension(SZK_(GV)+1) :: & @@ -672,7 +688,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) S_int, & ! Salinity interpolated to interfaces [S ~> ppt] H_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m] H_bot, & ! The distance of each filtered interface from the bottom [Z ~> m] - gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. + gprime, & ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. + N2 ! The Brunt Vaissalla freqency squared [T-2 ~> s-2] real, dimension(SZK_(GV),SZI_(G)) :: & Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m] Tf, & ! Layer temperatures after very thin layers are combined [C ~> degC] @@ -684,7 +701,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) Hc, & ! A column of layer thicknesses after convective instabilities are removed [Z ~> m] Tc, & ! A column of layer temperatures after convective instabilities are removed [C ~> degC] Sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt] - Rc ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] + Rc, & ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] + Hc_H ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2] real :: I_Htot ! The inverse of the total filtered thicknesses [Z ~> m] real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and its ! derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim]. @@ -737,7 +755,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) real :: tol_merge ! The fractional change in estimated wave speed that is allowed ! when deciding to merge layers in the calculation [nondim] integer :: kf(SZI_(G)) ! The number of active layers after filtering. - integer, parameter :: max_itt = 10 + integer, parameter :: max_itt = 30 logical :: use_EOS ! If true, density is calculated from T & S using the equation of state. logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed. logical :: merge ! If true, merge the current layer with the one above. @@ -749,6 +767,21 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) integer :: kc ! The number of layers in the column after merging integer :: sub, sub_it integer :: i, j, k, k2, itt, is, ie, js, je, nz, iint, m + real, dimension(SZK_(GV)+1) :: modal_structure !< Normalized model structure [nondim] + real, dimension(SZK_(GV)) :: modal_structure_fder !< Normalized model structure [Z-1 ~> m-1] + real :: mode_struct(SZK_(GV)+1) ! The mode structure [nondim], but it is also temporarily + ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6. + real :: mode_struct_fder(SZK_(GV)) ! The mode structure 1st derivative [nondim], but it is also temporarily + ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6. + real :: mode_struct_sq(SZK_(GV)+1) ! The square of mode structure [nondim] + real :: mode_struct_fder_sq(SZK_(GV)) ! The square of mode structure 1st derivative [Z-2 ~> m-2] + + + real :: ms_min, ms_max ! The minimum and maximum mode structure values returned from tdma6 [L2 T-2 ~> m2 s-2] + real :: ms_sq ! The sum of the square of the values returned from tdma6 [L4 T-4 ~> m4 s-4] + real :: w2avg ! A total for renormalization + real, parameter :: a_int = 0.5 ! Integral total for normalization + real :: renorm ! Normalization factor is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -777,9 +810,17 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) endif cg1_min2 = CS%min_speed2 - ! Zero out all wave speeds. Values over land or for columns that are too weakly stratified + ! Zero out all local values. Values over land or for columns that are too weakly stratified ! are not changed from this zero value. cn(:,:,:) = 0.0 + u_struct_max(:,:,:) = 0.0 + u_struct_bot(:,:,:) = 0.0 + Nb(:,:) = 0.0 + int_w2(:,:,:) = 0.0 + int_N2w2(:,:,:) = 0.0 + int_U2(:,:,:) = 0.0 + u_struct(:,:,:,:) = 0.0 + w_struct(:,:,:,:) = 0.0 min_h_frac = tol_Hfrac / real(nz) !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,CS,min_h_frac,use_EOS, & @@ -1010,8 +1051,13 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) ! Calculate Igu, Igl, depth, and N2 at each interior interface ! [excludes surface (K=1) and bottom (K=kc+1)] + Igl(:) = 0. + Igu(:) = 0. + N2(:) = 0. + do K=2,kc Igl(K) = 1.0/(gprime(K)*Hc(k)) ; Igu(K) = 1.0/(gprime(K)*Hc(k-1)) + N2(K) = US%L_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1))) if (better_est) then speed2_tot = speed2_tot + gprime(K)*((H_top(K) * H_bot(K)) * I_Htot) else @@ -1019,9 +1065,21 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) endif enddo + ! Set stratification for surface and bottom (setting equal to nearest interface for now) + N2(1) = N2(2) ; N2(kc+1) = N2(kc) + ! set bottom stratification + Nb(i,j) = sqrt(N2(kc+1)) + ! Under estimate the first eigenvalue (overestimate the speed) to start with. lam_1 = 1.0 / speed2_tot + ! init and first guess for mode structure + mode_struct(:) = 0. + mode_struct_fder(:) = 0. + mode_struct(2:kc) = 1. ! Uniform flow, first guess + modal_structure(:) = 0. + modal_structure_fder(:) = 0. + ! Find the first eigen value do itt=1,max_itt ! calculate the determinant of (A-lam_1*I) @@ -1039,11 +1097,89 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) lam_1 = lam_1 + dlam endif + call tdma6(kc-1, Igu(2:kc), Igl(2:kc), lam_1, mode_struct(2:kc)) + ! Note that tdma6 changes the units of mode_struct to [L2 T-2 ~> m2 s-2] + ! apply BC + mode_struct(1) = 0. + mode_struct(kc+1) = 0. + + ! renormalization of the integral of the profile + w2avg = 0.0 + do k=1,kc + w2avg = w2avg + 0.5*(mode_struct(K)**2+mode_struct(K+1)**2)*Hc(k) ![Z L4 T-4] + enddo + renorm = sqrt(htot(i)*a_int/w2avg) ![L-2 T-2] + do K=1,kc+1 ; mode_struct(K) = renorm * mode_struct(K) ; enddo + ! after renorm, mode_struct is again [nondim] + if (abs(dlam) < tol_solve*lam_1) exit enddo if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1) + ! sign of wave structure is irrelevant, flip to positive if needed + if (mode_struct(2)<0.) then + mode_struct(2:kc) = -1. * mode_struct(2:kc) + endif + + ! vertical derivative of w at interfaces lives on the layer points + do k=1,kc + mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / Hc(k) + enddo + + ! boundary condition for derivative is no-gradient + do k=kc+1,nz + mode_struct_fder(k) = mode_struct_fder(kc) + enddo + + ! now save maximum value and bottom value + u_struct_bot(i,j,1) = mode_struct_fder(kc) + u_struct_max(i,j,1) = maxval(abs(mode_struct_fder(1:kc))) + + ! Calculate terms for vertically integrated energy equation + do k=1,kc + mode_struct_fder_sq(k) = mode_struct_fder(k)**2 + enddo + do K=1,kc+1 + mode_struct_sq(K) = mode_struct(K)**2 + enddo + + ! sum over layers for quantities defined on layer + do k=1,kc + int_U2(i,j,1) = int_U2(i,j,1) + mode_struct_fder_sq(k) * Hc(k) + enddo + + ! vertical integration with Trapezoidal rule for values at interfaces + do K=1,kc + int_w2(i,j,1) = int_w2(i,j,1) + 0.5*(mode_struct_sq(K)+mode_struct_sq(K+1)) * Hc(k) + int_N2w2(i,j,1) = int_N2w2(i,j,1) + 0.5*(mode_struct_sq(K)*N2(K) + & + mode_struct_sq(K+1)*N2(K+1)) * Hc(k) + enddo + + ! Note that remapping_core_h requires that the same units be used + ! for both the source and target grid thicknesses, here [H ~> m or kg m-2]. + do k = 1,kc + Hc_H(k) = GV%Z_to_H * Hc(k) + enddo + + ! for w (diag) interpolate onto all interfaces + call interpolate_column(kc, Hc_H(1:kc), mode_struct(1:kc+1), & + nz, h(i,j,:), modal_structure(:), .false.) + + ! for u (remap) onto all layers + call remapping_core_h(CS%remapping_CS, kc, Hc_H(1:kc), mode_struct_fder(1:kc), & + nz, h(i,j,:), modal_structure_fder(:), & + GV%H_subroundoff, GV%H_subroundoff) + + ! write the wave structure + do k=1,nz+1 + w_struct(i,j,k,1) = modal_structure(k) + enddo + + do k=1,nz + u_struct(i,j,k,1) = modal_structure_fder(k) + enddo + ! Find other eigen values if c1 is of significant magnitude, > cn_thresh nrootsfound = 0 ! number of extra roots found (not including 1st root) if ((nmodes > 1) .and. (kc >= nmodes+1) .and. (cn(i,j,1) > CS%c1_thresh)) then @@ -1128,16 +1264,105 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) ! Use Newton's method to find the roots within the identified windows do m=1,nrootsfound ! loop over the root-containing widows (excluding 1st mode) lam_n = xbl(m) ! first guess is left edge of window + + ! init and first guess for mode structure + mode_struct(:) = 0. + mode_struct_fder(:) = 0. + mode_struct(2:kc) = 1. ! Uniform flow, first guess + modal_structure(:) = 0. + modal_structure_fder(:) = 0. + do itt=1,max_itt ! calculate the determinant of (A-lam_n*I) call tridiag_det(Igu, Igl, 2, kc, lam_n, det, ddet, row_scale=c2_scale) ! Use Newton's method to find a new estimate of lam_n dlam = - det / ddet lam_n = lam_n + dlam + + call tdma6(kc-1, Igu(2:kc), Igl(2:kc), lam_n, mode_struct(2:kc)) + ! Note that tdma6 changes the units of mode_struct to [L2 T-2 ~> m2 s-2] + ! apply BC + mode_struct(1) = 0. + mode_struct(kc+1) = 0. + + ! renormalization of the integral of the profile + w2avg = 0.0 + do k=1,kc + w2avg = w2avg + 0.5*(mode_struct(K)**2+mode_struct(K+1)**2)*Hc(k) + enddo + renorm = sqrt(htot(i)*a_int/w2avg) + do K=1,kc+1 ; mode_struct(K) = renorm * mode_struct(K) ; enddo + if (abs(dlam) < tol_solve*lam_1) exit enddo ! itt-loop + ! calculate nth mode speed if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n) + + ! sign is irrelevant, flip to positive if needed + if (mode_struct(2)<0.) then + mode_struct(2:kc) = -1. * mode_struct(2:kc) + endif + + ! derivative of vertical profile (i.e. dw/dz) is evaluated at the layer point + do k=1,kc + mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / Hc(k) + enddo + + ! boundary condition for 1st derivative is no-gradient + do k=kc+1,nz + mode_struct_fder(k) = mode_struct_fder(kc) + enddo + + ! now save maximum value and bottom value + u_struct_bot(i,j,m) = mode_struct_fder(kc) + u_struct_max(i,j,m) = maxval(abs(mode_struct_fder(1:kc))) + + ! Calculate terms for vertically integrated energy equation + do k=1,kc + mode_struct_fder_sq(k) = mode_struct_fder(k)**2 + enddo + do K=1,kc+1 + mode_struct_sq(K) = mode_struct(K)**2 + enddo + + ! sum over layers for integral of quantities defined at layer points + do k=1,kc + int_U2(i,j,m) = int_U2(i,j,m) + mode_struct_fder_sq(k) * Hc(k) + enddo + + ! vertical integration with Trapezoidal rule for quantities on interfaces + do K=1,kc + int_w2(i,j,m) = int_w2(i,j,m) + 0.5*(mode_struct_sq(K)+mode_struct_sq(K+1)) * Hc(k) + int_N2w2(i,j,m) = int_N2w2(i,j,m) + 0.5*(mode_struct_sq(K)*N2(K) + & + mode_struct_sq(K+1)*N2(K+1)) * Hc(k) + enddo + + ! Note that remapping_core_h requires that the same units be used + ! for both the source and target grid thicknesses, here [H ~> m or kg m-2]. + do k = 1,kc + Hc_H(k) = GV%Z_to_H * Hc(k) + enddo + + ! for w (diag) interpolate onto all interfaces + call interpolate_column(kc, Hc_H(1:kc), mode_struct(1:kc+1), & + nz, h(i,j,:), modal_structure(:), .false.) + + ! for u (remap) onto all layers + call remapping_core_h(CS%remapping_CS, kc, Hc_H(1:kc), mode_struct_fder(1:kc), & + nz, h(i,j,:), modal_structure_fder(:), & + GV%H_subroundoff, GV%H_subroundoff) + + ! write the wave structure + ! note that m=1 solves for 2nd mode,... + do k=1,nz+1 + w_struct(i,j,k,m+1) = modal_structure(k) + enddo + + do k=1,nz + u_struct(i,j,k,m+1) = modal_structure_fder(k) + enddo + enddo ! n-loop endif ! if nmodes>1 .and. kc>nmodes .and. c1>c1_thresh endif ! if more than 2 layers diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 6dda4c1b1c..d3f202339a 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -34,7 +34,7 @@ module MOM_internal_tides public get_lowmode_loss !> This control structure has parameters for the MOM_internal_tides module -type, public :: int_tide_CS ; private +type, public :: int_tide_CS logical :: do_int_tides !< If true, use the internal tide code. integer :: nFreq = 0 !< The number of internal tide frequency bands integer :: nMode = 1 !< The number of internal tide vertical modes @@ -95,6 +95,20 @@ module MOM_internal_tides !! summed over angle, frequency and mode [R Z3 T-3 ~> W m-2] real, allocatable, dimension(:,:) :: tot_allprocesses_loss !< Energy loss rates due to all processes, !! summed over angle, frequency and mode [R Z3 T-3 ~> W m-2] + real, allocatable, dimension(:,:,:,:) :: w_struct !< Vertical structure of vertical velocity (normalized) + !! for each frequency and each mode [nondim] + real, allocatable, dimension(:,:,:,:) :: u_struct !< Vertical structure of horizontal velocity (normalized and + !! divided by layer thicknesses) for each frequency and each mode [Z-1 ~> m-1] + real, allocatable, dimension(:,:,:) :: u_struct_max !< Maximum of u_struct, + !! for each mode [Z-1 ~> m-1] + real, allocatable, dimension(:,:,:) :: u_struct_bot !< Bottom value of u_struct, + !! for each mode [Z-1 ~> m-1] + real, allocatable, dimension(:,:,:) :: int_w2 !< Vertical integral of w_struct squared, + !! for each mode [Z ~> m] + real, allocatable, dimension(:,:,:) :: int_U2 !< Vertical integral of u_struct squared, + !! for each mode [Z-1 ~> m-1] + real, allocatable, dimension(:,:,:) :: int_N2w2 !< Depth-integrated Brunt Vaissalla freqency times + !! vertical profile squared, for each mode [Z T-2 ~> m s-2] real :: q_itides !< fraction of local dissipation [nondim] real :: En_sum !< global sum of energy for use in debugging, in MKS units [J] type(time_type), pointer :: Time => NULL() !< A pointer to the model's clock. @@ -126,7 +140,6 @@ module MOM_internal_tides type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate the !! timing of diagnostic output. - type(wave_structure_CS) :: wave_struct !< Wave structure control structure !>@{ Diag handles ! Diag handles relevant to all modes, frequencies, and angles @@ -148,6 +161,12 @@ module MOM_internal_tides integer, allocatable, dimension(:,:) :: & id_En_ang_mode, & id_itidal_loss_ang_mode + integer, allocatable, dimension(:) :: & + id_Ustruct_mode, & + id_Wstruct_mode, & + id_int_w2_mode, & + id_int_U2_mode, & + id_int_N2w2_mode !>@} end type int_tide_CS @@ -205,6 +224,10 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & real :: I_D_here ! The inverse of the local depth [Z-1 ~> m-1] real :: I_rho0 ! The inverse fo the Boussinesq density [R-1 ~> m3 kg-1] real :: freq2 ! The frequency squared [T-2 ~> s-2] + real :: PE_term ! total potential energy of profile [R Z ~> kg m-2] + real :: KE_term ! total kinetic energy of profile [R Z ~> kg m-2] + real :: U_mag ! rescaled magnitude of horizontal profile [L Z T-1 ~> m2 s-1] + real :: W0 ! rescaled magnitude of vertical profile [Z T-1 ~> m s-1] real :: c_phase ! The phase speed [L T-1 ~> m s-1] real :: loss_rate ! An energy loss rate [T-1 ~> s-1] real :: Fr2_max ! The column maximum internal wave Froude number squared [nondim] @@ -222,6 +245,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nAngle = CS%NAngle + nzm = GV%ke I_rho0 = 1.0 / GV%Rho0 cn_subRO = 1e-30*US%m_s_to_L_T en_subRO = 1e-30*US%W_m2_to_RZ3_T3*US%s_to_T @@ -229,6 +253,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! initialize local arrays drag_scale(:,:) = 0. Ub(:,:,:,:) = 0. + Umax(:,:,:,:) = 0. ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.********************** ! This is wrong, of course, but it works reasonably in some cases. @@ -417,15 +442,43 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! First, find velocity profiles if (CS%apply_wave_drag .or. CS%apply_Froude_drag) then do m=1,CS%NMode ; do fr=1,CS%Nfreq - ! Calculate modal structure for given mode and frequency - call wave_structure(h, tv, G, GV, US, cn(:,:,m), m, CS%frequency(fr), & - CS%wave_struct, tot_En_mode(:,:,fr,m), full_halos=.true.) - ! Pick out near-bottom and max horizontal baroclinic velocity values at each point + + ! compute near-bottom and max horizontal baroclinic velocity values at each point do j=jsd,jed ; do i=isd,ied id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging - nzm = CS%wave_struct%num_intfaces(i,j) - Ub(i,j,fr,m) = CS%wave_struct%Uavg_profile(i,j,nzm) - Umax(i,j,fr,m) = maxval(CS%wave_struct%Uavg_profile(i,j,1:nzm)) + + ! Calculate wavenumber magnitude + freq2 = CS%frequency(fr)**2 + + f2 = (0.25*(G%CoriolisBu(I,J) + G%CoriolisBu(max(I-1,1),max(J-1,1)) + & + G%CoriolisBu(I,max(J-1,1)) + G%CoriolisBu(max(I-1,1),J)))**2 + Kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subRO**2) + + + ! Back-calculate amplitude from energy equation + if ( (G%mask2dT(i,j) > 0.5) .and. (freq2*Kmag2 > 0.0)) then + ! Units here are [R Z ~> kg m-2] + KE_term = 0.25*GV%Rho0*( ((freq2 + f2) / (freq2*Kmag2))*US%L_to_Z**2*CS%int_U2(i,j,m) + & + CS%int_w2(i,j,m) ) + PE_term = 0.25*GV%Rho0*( CS%int_N2w2(i,j,m) / freq2 ) + + if (KE_term + PE_term > 0.0) then + W0 = sqrt( tot_En_mode(i,j,fr,m) / (KE_term + PE_term) ) + else + !call MOM_error(WARNING, "MOM internal tides: KE + PE <= 0.0; setting to W0 to 0.0") + W0 = 0.0 + endif + + U_mag = W0 * sqrt((freq2 + f2) / (2.0*freq2*Kmag2)) + ! scaled maximum tidal velocity + Umax(i,j,fr,m) = abs(U_mag * CS%u_struct_max(i,j,m)) + ! scaled bottom tidal velocity + Ub(i,j,fr,m) = abs(U_mag * CS%u_struct_bot(i,j,m)) + else + Umax(i,j,fr,m) = 0. + Ub(i,j,fr,m) = 0. + endif + enddo ; enddo ! i-loop, j-loop enddo ; enddo ! fr-loop, m-loop endif ! apply_wave or _Froude_drag (Ub or Umax needed) @@ -454,7 +507,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Pick out maximum baroclinic velocity values; calculate Fr=max(u)/cg do m=1,CS%NMode ; do fr=1,CS%Nfreq freq2 = CS%frequency(fr)**2 - do j=jsd,jed ; do i=isd,ied + do j=js,je ; do i=is,ie id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging ! Calculate horizontal phase velocity magnitudes f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & @@ -463,7 +516,6 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & c_phase = 0.0 if (Kmag2 > 0.0) then c_phase = sqrt(freq2/Kmag2) - nzm = CS%wave_struct%num_intfaces(i,j) Fr2_max = (Umax(i,j,fr,m) / c_phase)**2 ! Dissipate energy if Fr>1; done here with an arbitrary time scale if (Fr2_max > 1.0) then @@ -635,6 +687,26 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & call post_data(CS%id_Ub_mode(fr,m), Ub(:,:,fr,m), CS%diag) endif ; enddo ; enddo + do m=1,CS%NMode ; if (CS%id_Ustruct_mode(m) > 0) then + call post_data(CS%id_Ustruct_mode(m), CS%u_struct(:,:,:,m), CS%diag) + endif ; enddo + + do m=1,CS%NMode ; if (CS%id_Wstruct_mode(m) > 0) then + call post_data(CS%id_Wstruct_mode(m), CS%w_struct(:,:,:,m), CS%diag) + endif ; enddo + + do m=1,CS%NMode ; if (CS%id_int_w2_mode(m) > 0) then + call post_data(CS%id_int_w2_mode(m), CS%int_w2(:,:,m), CS%diag) + endif ; enddo + + do m=1,CS%NMode ; if (CS%id_int_U2_mode(m) > 0) then + call post_data(CS%id_int_U2_mode(m), CS%int_U2(:,:,m), CS%diag) + endif ; enddo + + do m=1,CS%NMode ; if (CS%id_int_N2w2_mode(m) > 0) then + call post_data(CS%id_int_N2w2_mode(m), CS%int_N2w2(:,:,m), CS%diag) + endif ; enddo + ! Output 2-D horizontal phase velocity for each frequency and mode do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_cp_mode(fr,m) > 0) then call post_data(CS%id_cp_mode(fr,m), CS%cp(:,:,fr,m), CS%diag) @@ -2226,7 +2298,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) ! nominal ocean depth, or a negative value for no limit [nondim] real :: period_1 ! The period of the gravest modeled mode [T ~> s] integer :: num_angle, num_freq, num_mode, m, fr - integer :: isd, ied, jsd, jed, a, id_ang, i, j + integer :: isd, ied, jsd, jed, a, id_ang, i, j, nz type(axes_grp) :: axes_ang ! This include declares and sets the variable "version". # include "version_variable.h" @@ -2241,6 +2313,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) character(len=80) :: rough_var ! Input file variable names isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + nz = GV%ke use_int_tides = .false. call read_param(param_file, "INTERNAL_TIDES", use_int_tides) @@ -2407,6 +2480,13 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%tot_itidal_loss(isd:ied,jsd:jed), source=0.0) allocate(CS%tot_Froude_loss(isd:ied,jsd:jed), source=0.0) allocate(CS%tot_residual_loss(isd:ied,jsd:jed), source=0.0) + allocate(CS%u_struct_bot(isd:ied,jsd:jed,num_mode), source=0.0) + allocate(CS%u_struct_max(isd:ied,jsd:jed,num_mode), source=0.0) + allocate(CS%int_w2(isd:ied,jsd:jed,num_mode), source=0.0) + allocate(CS%int_U2(isd:ied,jsd:jed,num_mode), source=0.0) + allocate(CS%int_N2w2(isd:ied,jsd:jed,num_mode), source=0.0) + allocate(CS%w_struct(isd:ied,jsd:jed,1:nz+1,num_mode), source=0.0) + allocate(CS%u_struct(isd:ied,jsd:jed,1:nz,num_mode), source=0.0) ! Compute the fixed part of the bottom drag loss from baroclinic modes call get_param(param_file, mdl, "H2_FILE", h2_file, & @@ -2593,6 +2673,11 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%id_allprocesses_loss_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_itidal_loss_ang_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_Ub_mode(CS%nFreq,CS%nMode), source=-1) + allocate(CS%id_Ustruct_mode(CS%nMode), source=-1) + allocate(CS%id_Wstruct_mode(CS%nMode), source=-1) + allocate(CS%id_int_w2_mode(CS%nMode), source=-1) + allocate(CS%id_int_U2_mode(CS%nMode), source=-1) + allocate(CS%id_int_N2w2_mode(CS%nMode), source=-1) allocate(CS%id_cp_mode(CS%nFreq,CS%nMode), source=-1) allocate(angles(CS%NAngle), source=0.0) @@ -2656,8 +2741,42 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) enddo ; enddo - ! Initialize wave_structure (not sure if this should be here - BDM) - call wave_structure_init(Time, G, GV, param_file, diag, CS%wave_struct) + + do m=1,CS%nMode + + ! Register 3-D internal tide horizonal velocity profile for each mode + write(var_name, '("Itide_Ustruct","_mode",i1)') m + write(var_descript, '("horizonal velocity profile for mode ",i1)') m + CS%id_Ustruct_mode(m) = register_diag_field('ocean_model', var_name, & + diag%axesTl, Time, var_descript, 'm-1', conversion=US%m_to_L) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + + ! Register 3-D internal tide vertical velocity profile for each mode + write(var_name, '("Itide_Wstruct","_mode",i1)') m + write(var_descript, '("vertical velocity profile for mode ",i1)') m + CS%id_Wstruct_mode(m) = register_diag_field('ocean_model', var_name, & + diag%axesTi, Time, var_descript, '[]') + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + + write(var_name, '("Itide_int_w2","_mode",i1)') m + write(var_descript, '("integral of w2 for mode ",i1)') m + CS%id_int_w2_mode(m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'm', conversion=US%Z_to_m) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + + write(var_name, '("Itide_int_U2","_mode",i1)') m + write(var_descript, '("integral of U2 for mode ",i1)') m + CS%id_int_U2_mode(m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'm-1', conversion=US%m_to_L) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + + write(var_name, '("Itide_int_N2w2","_mode",i1)') m + write(var_descript, '("integral of N2w2 for mode ",i1)') m + CS%id_int_N2w2_mode(m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'm s-2', conversion=US%Z_to_m*US%s_to_T**2) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + + enddo end subroutine internal_tides_init @@ -2670,6 +2789,12 @@ subroutine internal_tides_end(CS) if (allocated(CS%id_En_mode)) deallocate(CS%id_En_mode) if (allocated(CS%id_Ub_mode)) deallocate(CS%id_Ub_mode) if (allocated(CS%id_cp_mode)) deallocate(CS%id_cp_mode) + if (allocated(CS%id_Ustruct_mode)) deallocate(CS%id_Ustruct_mode) + if (allocated(CS%id_Wstruct_mode)) deallocate(CS%id_Wstruct_mode) + if (allocated(CS%id_int_w2_mode)) deallocate(CS%id_int_w2_mode) + if (allocated(CS%id_int_U2_mode)) deallocate(CS%id_int_U2_mode) + if (allocated(CS%id_int_N2w2_mode)) deallocate(CS%id_int_N2w2_mode) + end subroutine internal_tides_end end module MOM_internal_tides diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 0fe08a06b2..b0d04e434c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -396,7 +396,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%uniform_test_cg > 0.0) then do m=1,CS%nMode ; cn_IGW(:,:,m) = CS%uniform_test_cg ; enddo else - call wave_speeds(h, tv, G, GV, US, CS%nMode, cn_IGW, CS%wave_speed, full_halos=.true.) + call wave_speeds(h, tv, G, GV, US, CS%nMode, cn_IGW, CS%wave_speed, CS%int_tide%w_struct, & + CS%int_tide%u_struct, CS%int_tide%u_struct_max, CS%int_tide%u_struct_bot, & + CS%int_tide_input%Nb, CS%int_tide%int_w2, CS%int_tide%int_U2, CS%int_tide%int_N2w2, & + full_halos=.true.) endif call propagate_int_tide(h, tv, cn_IGW, CS%int_tide_input%TKE_itidal_input, CS%int_tide_input%tideamp, & From e45b983eee130982a448cac211a5d0141f352046 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 24 Apr 2023 22:01:09 -0400 Subject: [PATCH 07/15] remove wave_structure broken code --- src/diagnostics/MOM_wave_structure.F90 | 793 ------------------ .../lateral/MOM_internal_tides.F90 | 1 - 2 files changed, 794 deletions(-) delete mode 100644 src/diagnostics/MOM_wave_structure.F90 diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 deleted file mode 100644 index 80d23eeb75..0000000000 --- a/src/diagnostics/MOM_wave_structure.F90 +++ /dev/null @@ -1,793 +0,0 @@ -!> Vertical structure functions for first baroclinic mode wave speed -module MOM_wave_structure - -! This file is part of MOM6. See LICENSE.md for the license. - -! By Benjamin Mater & Robert Hallberg, 2015 - -! The subroutine in this module calculates the vertical structure -! functions of the first baroclinic mode internal wave speed. -! Calculation of interface values is the same as done in -! MOM_wave_speed by Hallberg, 2008. - -use MOM_debugging, only : isnan => is_NaN -use MOM_checksums, only : chksum0, hchksum -use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl -use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type -use MOM_EOS, only : calculate_density_derivs -use MOM_error_handler, only : MOM_error, FATAL, WARNING -use MOM_file_parser, only : log_version, param_file_type, get_param -use MOM_grid, only : ocean_grid_type -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type -use regrid_solvers, only : solve_diag_dominant_tridiag - -implicit none ; private - -#include - -public wave_structure, wave_structure_init - -! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional -! consistency testing. These are noted in comments with units like Z, H, L, and T, along with -! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units -! vary with the Boussinesq approximation, the Boussinesq variant is given first. - -!> The control structure for the MOM_wave_structure module -type, public :: wave_structure_CS ; !private - logical :: initialized = .false. !< True if this control structure has been initialized. - type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to - !! regulate the timing of diagnostic output. - real, allocatable, dimension(:,:,:) :: w_strct - !< Vertical structure of vertical velocity (normalized) [nondim]. - real, allocatable, dimension(:,:,:) :: u_strct - !< Vertical structure of horizontal velocity (normalized and - !! divided by layer thicknesses) [Z-1 ~> m-1]. - real, allocatable, dimension(:,:,:) :: W_profile - !< Vertical profile of w_hat(z), where - !! w(x,y,z,t) = w_hat(z)*exp(i(kx+ly-freq*t)) is the full time- - !! varying vertical velocity with w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1]. - real, allocatable, dimension(:,:,:) :: Uavg_profile - !< Vertical profile of the magnitude of horizontal velocity, - !! (u^2+v^2)^0.5, averaged over a period [L T-1 ~> m s-1]. - real, allocatable, dimension(:,:,:) :: z_depths - !< Depths of layer interfaces [Z ~> m]. - real, allocatable, dimension(:,:,:) :: N2 - !< Squared buoyancy frequency at each interface [T-2 ~> s-2]. - integer, allocatable, dimension(:,:):: num_intfaces - !< Number of layer interfaces (including surface and bottom) [nondim]. - ! logical :: int_tide_source_test !< If true, apply an arbitrary generation site for internal tide testing - ! integer :: int_tide_source_i !< I Location of generation site - ! integer :: int_tide_source_j !< J Location of generation site - logical :: debug !< debugging prints - -end type wave_structure_CS - -contains - -!> This subroutine determines the internal wave velocity structure for any mode. -!! -!! This subroutine solves for the eigen vector [vertical structure, e(k)] associated with -!! the first baroclinic mode speed [i.e., smallest eigen value (lam = 1/c^2)] of the -!! system d2e/dz2 = -(N2/cn2)e, or (A-lam*I)e = 0, where A = -(1/N2)(d2/dz2), lam = 1/c^2, -!! and I is the identity matrix. 2nd order discretization in the vertical lets this system -!! be represented as -!! -!! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0 -!! -!! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving -!! -!! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0 -!! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0 -!! -!! where, upon noting N2 = reduced gravity/layer thickness, we get -!! Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1)) -!! -!! The eigen value for this system is approximated using "wave_speed." This subroutine uses -!! these eigen values (mode speeds) to estimate the corresponding eigen vectors (velocity -!! structure) using the "inverse iteration with shift" method. The algorithm is -!! -!! Pick a starting vector reasonably close to mode structure and with unit magnitude, b_guess -!! For n=1,2,3,... -!! Solve (A-lam*I)e = e_guess for e -!! Set e_guess=e/|e| and repeat, with each iteration refining the estimate of e -subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halos) - type(ocean_grid_type), intent(in) :: 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 - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various - !! thermodynamic variables. - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: cn !< The (non-rotational) mode internal - !! gravity wave speed [L T-1 ~> m s-1]. - integer, intent(in) :: ModeNum !< Mode number - real, intent(in) :: freq !< Intrinsic wave frequency [T-1 ~> s-1]. - type(wave_structure_CS), intent(inout) :: CS !< Wave structure control struct - real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: En !< Internal wave energy density [R Z3 T-2 ~> J m-2] - logical, optional, intent(in) :: full_halos !< If true, do the calculation - !! over the entire computational domain. - ! Local variables - real, dimension(SZK_(GV)+1) :: & - dRho_dT, & !< Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] - dRho_dS, & !< Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1] - pres, & !< Interface pressure [R L2 T-2 ~> Pa] - T_int, & !< Temperature interpolated to interfaces [C ~> degC] - S_int, & !< Salinity interpolated to interfaces [S ~> ppt] - gprime !< The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. - real, dimension(SZK_(GV)) :: & - Igl, Igu !< The inverse of the reduced gravity across an interface times - !< the thickness of the layer below (Igl) or above (Igu) it [T2 L-2 ~> s2 m-2]. - real, dimension(SZK_(GV),SZI_(G)) :: & - Hf, & !< Layer thicknesses after very thin layers are combined [Z ~> m] - Tf, & !< Layer temperatures after very thin layers are combined [C ~> degC] - Sf, & !< Layer salinities after very thin layers are combined [S ~> ppt] - Rf !< Layer densities after very thin layers are combined [R ~> kg m-3] - real, dimension(SZK_(GV)) :: & - Hc, & !< A column of layer thicknesses after convective instabilities are removed [Z ~> m] - Tc, & !< A column of layer temperatures after convective instabilities are removed [C ~> degC] - Sc, & !< A column of layer salinities after convective instabilities are removed [S ~> ppt] - Rc !< A column of layer densities after convective instabilities are removed [R ~> kg m-3] - real, dimension(SZI_(G),SZJ_(G)) :: & - htot !< The vertical sum of the thicknesses [Z ~> m] - real :: lam !< inverse of wave speed squared [T2 L-2 ~> s2 m-2] - real :: min_h_frac !< fractional (per layer) minimum thickness [nondim] - real :: Z_to_pres !< A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] - real, dimension(SZI_(G)) :: & - hmin, & !< Thicknesses [Z ~> m] - H_here, & !< A thickness [Z ~> m] - HxT_here, & !< A layer integrated temperature [C Z ~> degC m] - HxS_here, & !< A layer integrated salinity [S Z ~> ppt m] - HxR_here !< A layer integrated density [R Z ~> kg m-2] - real :: I_Hnew !< The inverse of a new layer thickness [Z-1 ~> m-1] - real :: drxh_sum !< The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2] - real, parameter :: tol1 = 0.0001, tol2 = 0.001 ! Nondimensional tolerances [nondim] - real :: g_Rho0 !< G_Earth/Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. - ! real :: rescale, I_rescale - integer :: kf(SZI_(G)) - integer, parameter :: max_itt = 1 !< number of times to iterate in solving for eigenvector - real :: cg_subRO !< A tiny wave speed to prevent division by zero [L T-1 ~> m s-1] - real, parameter :: a_int = 0.5 !< value of normalized integral: \int(w_strct^2)dz = a_int [nondim] - real :: I_a_int !< inverse of a_int [nondim] - real :: f2 !< squared Coriolis frequency [T-2 ~> s-2] - real :: Kmag2 !< magnitude of horizontal wave number squared [L-2 ~> m-2] - real :: emag2 ! The sum of the squared magnitudes of the guesses [nondim] - real :: pi_htot ! The gravest vertical wavenumber in this column [Z-1 ~> m-1] - real :: renorm ! A renormalization factor [nondim] - logical :: use_EOS !< If true, density is calculated from T & S using an - !! equation of state. - - ! local representations of variables in CS; note, - ! not all rows will be filled if layers get merged! - real, dimension(SZK_(GV)+1) :: w_strct !< Vertical structure of vertical velocity (normalized) [nondim]. - real, dimension(SZK_(GV)+1) :: u_strct !< Vertical structure of horizontal velocity (normalized and - !! divided by layer thicknesses) [Z-1 ~> m-1]. - real, dimension(SZK_(GV)+1) :: W_profile !< Vertical profile of w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1]. - real, dimension(SZK_(GV)+1) :: Uavg_profile !< Vertical profile of the magnitude of - !! horizontal velocity [L T-1 ~> m s-1]. - real, dimension(SZK_(GV)+1) :: z_int !< Integrated depth [Z ~> m] - real, dimension(SZK_(GV)+1) :: N2 !< Squared buoyancy frequency at each interface [T-2 ~> s-2]. - real, dimension(SZK_(GV)+1) :: w_strct2 !< squared values [nondim] - real, dimension(SZK_(GV)+1) :: u_strct2 !< squared values [Z-2 ~> m-2] - real, dimension(SZK_(GV)) :: dz !< thicknesses of merged layers (same as Hc I hope) [Z ~> m] - ! real, dimension(SZK_(GV)+1) :: dWdz_profile !< profile of dW/dz times total depth [Z T-1 ~> m s-1] - real :: w2avg !< average of squared vertical velocity structure function [Z ~> m] - real :: int_dwdz2 !< Vertical integral of the square of u_strct [Z-1 ~> m-1] - real :: int_w2 !< Vertical integral of the square of w_strct [Z ~> m] - real :: int_N2w2 !< Vertical integral of N2 [Z T-2 ~> m s-2] - real :: KE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2] - real :: PE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2] - real :: W0 !< A vertical velocity magnitude [Z T-1 ~> m s-1] - real :: U_mag !< A horizontal velocity magnitude times the depth of the - !! ocean [Z L T-1 ~> m2 s-1] - real, dimension(SZK_(GV)-1) :: lam_z !< product of eigen value and gprime(k); one value for each - !< interface (excluding surface and bottom) [Z-1 ~> m-1] - real, dimension(SZK_(GV)-1) :: a_diag !< upper diagonal of tridiagonal matrix; one value for each - !< interface (excluding surface and bottom) [Z-1 ~> m-1] - real, dimension(SZK_(GV)-1) :: c_diag !< lower diagonal of tridiagonal matrix; one value for each - !< interface (excluding surface and bottom) [Z-1 ~> m-1] - real, dimension(SZK_(GV)-1) :: b_dom !< Matrix center diagonal offset from a_diag + c_diag; one value - !< for each interface (excluding surface and bottom) [Z-1 ~> m-1] - real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitude (for TDMA) [nondim] - real, dimension(SZK_(GV)-1) :: e_itt !< improved guess at eigen vector (from TDMA) [nondim] - real :: Pi ! 3.1415926535... [nondim] - integer :: i, j, k, k2, kc, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop - - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - I_a_int = 1/a_int - - if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_structure: "// & - "Module must be initialized before it is used.") - - if (present(full_halos)) then ; if (full_halos) then - is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed - endif ; endif - - Pi = (4.0*atan(1.0)) - - g_Rho0 = GV%g_Earth / GV%Rho0 - - !if (CS%debug) call chksum0(g_Rho0, "g/rho0 in wave struct", & - ! scale=(US%L_to_m**2)*US%m_to_Z*(US%s_to_T**2)*US%kg_m3_to_R) - - if (CS%debug) call chksum0(freq, "freq in wave_struct", scale=US%s_to_T) - - cg_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. - use_EOS = associated(tv%eqn_of_state) - - ! Simplifying the following could change answers at roundoff. - Z_to_pres = GV%Z_to_H * (GV%H_to_RZ * GV%g_Earth) - ! rescale = 1024.0**4 ; I_rescale = 1.0/rescale - - min_h_frac = tol1 / real(nz) - - do j=js,je - ! First merge very thin layers with the one above (or below if they are - ! at the top). This also transposes the row order so that columns can - ! be worked upon one at a time. - do i=is,ie ; htot(i,j) = 0.0 ; enddo - do k=1,nz ; do i=is,ie ; htot(i,j) = htot(i,j) + h(i,j,k)*GV%H_to_Z ; enddo ; enddo - - do i=is,ie - hmin(i) = htot(i,j)*min_h_frac ; kf(i) = 1 ; H_here(i) = 0.0 - HxT_here(i) = 0.0 ; HxS_here(i) = 0.0 ; HxR_here(i) = 0.0 - enddo - if (use_EOS) then - do k=1,nz ; do i=is,ie - if ((H_here(i) > hmin(i)) .and. (h(i,j,k)*GV%H_to_Z > hmin(i))) then - Hf(kf(i),i) = H_here(i) - Tf(kf(i),i) = HxT_here(i) / H_here(i) - Sf(kf(i),i) = HxS_here(i) / H_here(i) - kf(i) = kf(i) + 1 - - ! Start a new layer - H_here(i) = h(i,j,k)*GV%H_to_Z - HxT_here(i) = (h(i,j,k) * GV%H_to_Z) * tv%T(i,j,k) - HxS_here(i) = (h(i,j,k) * GV%H_to_Z) * tv%S(i,j,k) - else - H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z - HxT_here(i) = HxT_here(i) + (h(i,j,k) * GV%H_to_Z) * tv%T(i,j,k) - HxS_here(i) = HxS_here(i) + (h(i,j,k) * GV%H_to_Z) * tv%S(i,j,k) - endif - enddo ; enddo - do i=is,ie ; if (H_here(i) > 0.0) then - Hf(kf(i),i) = H_here(i) - Tf(kf(i),i) = HxT_here(i) / H_here(i) - Sf(kf(i),i) = HxS_here(i) / H_here(i) - endif ; enddo - else - do k=1,nz ; do i=is,ie - if ((H_here(i) > hmin(i)) .and. (h(i,j,k)*GV%H_to_Z > hmin(i))) then - Hf(kf(i),i) = H_here(i) ; Rf(kf(i),i) = HxR_here(i) / H_here(i) - kf(i) = kf(i) + 1 - - ! Start a new layer - H_here(i) = h(i,j,k)*GV%H_to_Z - HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k) - else - H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z - HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k) - endif - enddo ; enddo - do i=is,ie ; if (H_here(i) > 0.0) then - Hf(kf(i),i) = H_here(i) ; Rf(kf(i),i) = HxR_here(i) / H_here(i) - endif ; enddo - endif ! use_EOS? - - ! From this point, we can work on individual columns without causing memory - ! to have page faults. - do i=is,ie ; if (cn(i,j) > 0.0) then - !----for debugging, remove later---- - ig = i + G%idg_offset ; jg = j + G%jdg_offset - !if (ig == CS%int_tide_source_i .and. jg == CS%int_tide_source_j) then - !----------------------------------- - if (G%mask2dT(i,j) > 0.0) then - - gprime(:) = 0.0 ! init gprime - pres(:) = 0.0 ! init pres - lam = 1/(cn(i,j)**2) - - ! Calculate drxh_sum - if (use_EOS) then - pres(1) = 0.0 - do k=2,kf(i) - pres(k) = pres(k-1) + Z_to_pres*Hf(k-1,i) - T_int(k) = 0.5*(Tf(k,i)+Tf(k-1,i)) - S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) - enddo - call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & - tv%eqn_of_state, (/2,kf(i)/) ) - - ! Sum the reduced gravities to find out how small a density difference - ! is negligibly small. - drxh_sum = 0.0 - do k=2,kf(i) - drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & - max(0.0,dRho_dT(k)*(Tf(k,i)-Tf(k-1,i)) + & - dRho_dS(k)*(Sf(k,i)-Sf(k-1,i))) - enddo - else - drxh_sum = 0.0 - do k=2,kf(i) - drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & - max(0.0,Rf(k,i)-Rf(k-1,i)) - enddo - endif ! use_EOS? - - ! Find gprime across each internal interface, taking care of convective - ! instabilities by merging layers. - if (drxh_sum >= 0.0) then - ! Merge layers to eliminate convective instabilities or exceedingly - ! small reduced gravities. - if (use_EOS) then - kc = 1 - Hc(1) = Hf(1,i) ; Tc(1) = Tf(1,i) ; Sc(1) = Sf(1,i) - do k=2,kf(i) - if ((dRho_dT(k)*(Tf(k,i)-Tc(kc)) + dRho_dS(k)*(Sf(k,i)-Sc(kc))) * & - (Hc(kc) + Hf(k,i)) < 2.0 * tol2*drxh_sum) then - ! Merge this layer with the one above and backtrack. - I_Hnew = 1.0 / (Hc(kc) + Hf(k,i)) - Tc(kc) = (Hc(kc)*Tc(kc) + Hf(k,i)*Tf(k,i)) * I_Hnew - Sc(kc) = (Hc(kc)*Sc(kc) + Hf(k,i)*Sf(k,i)) * I_Hnew - Hc(kc) = (Hc(kc) + Hf(k,i)) - ! Backtrack to remove any convective instabilities above... Note - ! that the tolerance is a factor of two larger, to avoid limit how - ! far back we go. - do k2=kc,2,-1 - if ((dRho_dT(k2)*(Tc(k2)-Tc(k2-1)) + dRho_dS(k2)*(Sc(k2)-Sc(k2-1))) * & - (Hc(k2) + Hc(k2-1)) < tol2*drxh_sum) then - ! Merge the two bottommost layers. At this point kc = k2. - I_Hnew = 1.0 / (Hc(kc) + Hc(kc-1)) - Tc(kc-1) = (Hc(kc)*Tc(kc) + Hc(kc-1)*Tc(kc-1)) * I_Hnew - Sc(kc-1) = (Hc(kc)*Sc(kc) + Hc(kc-1)*Sc(kc-1)) * I_Hnew - Hc(kc-1) = (Hc(kc) + Hc(kc-1)) - kc = kc - 1 - else ; exit ; endif - enddo - else - ! Add a new layer to the column. - kc = kc + 1 - drho_dS(kc) = drho_dS(k) ; drho_dT(kc) = drho_dT(k) - Tc(kc) = Tf(k,i) ; Sc(kc) = Sf(k,i) ; Hc(kc) = Hf(k,i) - endif - enddo - ! At this point there are kc layers and the gprimes should be positive. - do k=2,kc ! Revisit this if non-Boussinesq. - gprime(k) = g_Rho0 * (dRho_dT(k)*(Tc(k)-Tc(k-1)) + & - dRho_dS(k)*(Sc(k)-Sc(k-1))) - enddo - else ! .not.use_EOS - ! Do the same with density directly... - kc = 1 - Hc(1) = Hf(1,i) ; Rc(1) = Rf(1,i) - do k=2,kf(i) - if ((Rf(k,i) - Rc(kc)) * (Hc(kc) + Hf(k,i)) < 2.0*tol2*drxh_sum) then - ! Merge this layer with the one above and backtrack. - Rc(kc) = (Hc(kc)*Rc(kc) + Hf(k,i)*Rf(k,i)) / (Hc(kc) + Hf(k,i)) - Hc(kc) = (Hc(kc) + Hf(k,i)) - ! Backtrack to remove any convective instabilities above... Note - ! that the tolerance is a factor of two larger, to avoid limit how - ! far back we go. - do k2=kc,2,-1 - if ((Rc(k2)-Rc(k2-1)) * (Hc(k2)+Hc(k2-1)) < tol2*drxh_sum) then - ! Merge the two bottommost layers. At this point kc = k2. - Rc(kc-1) = (Hc(kc)*Rc(kc) + Hc(kc-1)*Rc(kc-1)) / (Hc(kc) + Hc(kc-1)) - Hc(kc-1) = (Hc(kc) + Hc(kc-1)) - kc = kc - 1 - else ; exit ; endif - enddo - else - ! Add a new layer to the column. - kc = kc + 1 - Rc(kc) = Rf(k,i) ; Hc(kc) = Hf(k,i) - endif - enddo - ! At this point there are kc layers and the gprimes should be positive. - do k=2,kc ! Revisit this if non-Boussinesq. - gprime(k) = g_Rho0 * (Rc(k)-Rc(k-1)) - enddo - endif ! use_EOS? - - !-----------------NOW FIND WAVE STRUCTURE------------------------------------- - ! Construct and solve tridiagonal system for the interior interfaces - ! Note that kc = number of layers, - ! kc+1 = nzm = number of interfaces, - ! kc-1 = number of interior interfaces (excluding surface and bottom) - ! Also, note that "K" refers to an interface, while "k" refers to the layer below. - ! Need at least 3 layers (2 internal interfaces) to generate a matrix, also - ! need number of layers to be greater than the mode number - if (kc >= max(3, ModeNum + 1)) then - ! Set depth at surface - z_int(1) = 0.0 - ! Calculate Igu, Igl, depth, and N2 at each interior interface - ! [excludes surface (K=1) and bottom (K=kc+1)] - do K=2,kc - Igl(K) = 1.0/(gprime(K)*Hc(k)) ; Igu(K) = 1.0/(gprime(K)*Hc(k-1)) - z_int(K) = z_int(K-1) + Hc(k-1) - N2(K) = US%L_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1))) - enddo - ! Set stratification for surface and bottom (setting equal to nearest interface for now) - N2(1) = N2(2) ; N2(kc+1) = N2(kc) - ! Calcualte depth at bottom - z_int(kc+1) = z_int(kc)+Hc(kc) - ! check that thicknesses sum to total depth - if (abs(z_int(kc+1)-htot(i,j)) > 1.e-14*htot(i,j)) then - call MOM_error(FATAL, "wave_structure: mismatch in total depths") - endif - - ! Populate interior rows of tridiagonal matrix; must multiply through by - ! gprime to get tridiagonal matrix to the symmetrical form: - ! [-1/H(k-1)]e(k-1) + [1/H(k-1)+1/H(k)-lam_z]e(k) + [-1/H(k)]e(k+1) = 0, - ! where lam_z = lam*gprime is now a function of depth. - ! First, populate interior rows - - ! init the values in matrix: since number of layers is variable, values need to be reset - lam_z(:) = 0.0 - a_diag(:) = 0.0 - b_dom(:) = 0.0 - c_diag(:) = 0.0 - e_guess(:) = 0.0 - e_itt(:) = 0.0 - w_strct(:) = 0.0 - do K=3,kc-1 - row = K-1 ! indexing for TD matrix rows - lam_z(row) = lam*gprime(K) - a_diag(row) = gprime(K)*(-Igu(K)) - b_dom(row) = 2.0*gprime(K)*(Igu(K)+Igl(K)) - lam_z(row) - c_diag(row) = gprime(K)*(-Igl(K)) - enddo - if (CS%debug) then ; do row=2,kc-2 - if (isnan(lam_z(row)))then ; print *, "Wave_structure: lam_z(row) is NAN" ; endif - if (isnan(a_diag(row)))then ; print *, "Wave_structure: a(k) is NAN" ; endif - if (isnan(c_diag(row)))then ; print *, "Wave_structure: c(k) is NAN" ; endif - enddo ; endif - ! Populate top row of tridiagonal matrix - K=2 ; row = K-1 ; - lam_z(row) = lam*gprime(K) - a_diag(row) = 0.0 - b_dom(row) = gprime(K)*(Igu(K)+2.0*Igl(K)) - lam_z(row) - c_diag(row) = gprime(K)*(-Igl(K)) - ! Populate bottom row of tridiagonal matrix - K=kc ; row = K-1 - lam_z(row) = lam*gprime(K) - a_diag(row) = gprime(K)*(-Igu(K)) - b_dom(row) = gprime(K)*(2.0*Igu(K) + Igl(K)) - lam_z(row) - c_diag(row) = 0.0 - - ! Guess a normalized vector shape to start with (excludes surface and bottom) - emag2 = 0.0 - pi_htot = Pi / htot(i,j) - do K=2,kc - e_guess(K-1) = sin(pi_htot * z_int(K)) - emag2 = emag2 + e_guess(K-1)**2 - enddo - renorm = 1.0 / sqrt(emag2) - do K=2,kc ; e_guess(K-1) = renorm*e_guess(K-1) ; enddo - - ! Perform inverse iteration with tri-diag solver - do itt=1,max_itt - ! this solver becomes unstable very quickly - ! b_diag(1:kc-1) = b_dom(1:kc-1) - (a_diag(1:kc-1) + c_diag(1:kc-1)) - !call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), & - ! -lam_z(1:kc-1),e_guess(1:kc-1),"TDMA_T",e_itt) - - call solve_diag_dominant_tridiag( c_diag, b_dom, a_diag, e_guess, e_itt, kc-1 ) - ! Renormalize the guesses of the structure.- - emag2 = 0.0 - do K=2,kc ; emag2 = emag2 + e_itt(K-1)**2 ; enddo - renorm = 1.0 / sqrt(emag2) - do K=2,kc ; e_guess(K-1) = renorm*e_itt(K-1) ; enddo - - ! A test should be added here to evaluate convergence. - enddo ! itt-loop - do K=2,kc ; w_strct(K) = e_guess(K-1) ; enddo - w_strct(1) = 0.0 ! rigid lid at surface - w_strct(kc+1) = 0.0 ! zero-flux at bottom - - ! Check to see if solver worked - if (CS%debug) then - ig_stop = 0 ; jg_stop = 0 - if (isnan(sum(w_strct(1:kc+1)))) then - print *, "Wave_structure: w_strct has a NAN at ig=", ig, ", jg=", jg - if (iG%iec .or. jG%jec)then - print *, "This is occuring at a halo point." - endif - ig_stop = ig ; jg_stop = jg - endif - endif - - ! Normalize vertical structure function of w such that - ! \int(w_strct)^2dz = a_int (a_int could be any value, e.g., 0.5) - nzm = kc+1 ! number of layer interfaces after merging - !(including surface and bottom) - w2avg = 0.0 - do k=1,nzm-1 - dz(k) = Hc(k) - w2avg = w2avg + 0.5*(w_strct(K)**2+w_strct(K+1)**2)*dz(k) - enddo - ! correct renormalization: - renorm = sqrt(htot(i,j)*a_int/w2avg) - do K=1,kc+1 ; w_strct(K) = renorm * w_strct(K) ; enddo - - ! Calculate vertical structure function of u (i.e. dw/dz) - do K=2,nzm-1 - u_strct(K) = 0.5*((w_strct(K-1) - w_strct(K) )/dz(k-1) + & - (w_strct(K) - w_strct(K+1))/dz(k)) - enddo - u_strct(1) = (w_strct(1) - w_strct(2) )/dz(1) - u_strct(nzm) = (w_strct(nzm-1)- w_strct(nzm))/dz(nzm-1) - - ! Calculate wavenumber magnitude - f2 = (0.25*(G%CoriolisBu(I,J) + G%CoriolisBu(max(I-1,1),max(J-1,1)) + & - G%CoriolisBu(I,max(J-1,1)) + G%CoriolisBu(max(I-1,1),J)))**2 - Kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cg_subRO**2) - - ! Calculate terms in vertically integrated energy equation - int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_N2w2 = 0.0 - do K=1,nzm - u_strct2(K) = u_strct(K)**2 - w_strct2(K) = w_strct(K)**2 - enddo - ! vertical integration with Trapezoidal rule - do k=1,nzm-1 - int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(K)+u_strct2(K+1)) * dz(k) - int_w2 = int_w2 + 0.5*(w_strct2(K)+w_strct2(K+1)) * dz(k) - int_N2w2 = int_N2w2 + 0.5*(w_strct2(K)*N2(K)+w_strct2(K+1)*N2(K+1)) * dz(k) - enddo - - ! Back-calculate amplitude from energy equation - if (present(En) .and. (freq**2*Kmag2 > 0.0)) then - ! Units here are [R Z ~> kg m-2] - KE_term = 0.25*GV%Rho0*( ((freq**2 + f2) / (freq**2*Kmag2))*US%L_to_Z**2*int_dwdz2 + int_w2 ) - PE_term = 0.25*GV%Rho0*( int_N2w2 / freq**2 ) - if (En(i,j) >= 0.0) then - W0 = sqrt( En(i,j) / (KE_term + PE_term) ) - else - call MOM_error(WARNING, "wave_structure: En < 0.0; setting to W0 to 0.0") - print *, "En(i,j)=", En(i,j), " at ig=", ig, ", jg=", jg - W0 = 0.0 - endif - ! Calculate actual vertical velocity profile and derivative - U_mag = W0 * sqrt((freq**2 + f2) / (2.0*freq**2*Kmag2)) - do K=1,nzm - W_profile(K) = W0*w_strct(K) - ! dWdz_profile(K) = W0*u_strct(K) - ! Calculate average magnitude of actual horizontal velocity over a period - Uavg_profile(K) = abs(U_mag * u_strct(K)) - enddo - else - do K=1,nzm - W_profile(K) = 0.0 - ! dWdz_profile(K) = 0.0 - Uavg_profile(K) = 0.0 - enddo - endif - - ! Store values in control structure - do K=1,nzm - CS%w_strct(i,j,K) = w_strct(K) - CS%u_strct(i,j,K) = u_strct(K) - CS%W_profile(i,j,K) = W_profile(K) - CS%Uavg_profile(i,j,K) = Uavg_profile(K) - CS%z_depths(i,j,K) = z_int(K) - CS%N2(i,j,K) = N2(K) - enddo - CS%num_intfaces(i,j) = nzm - else - ! If not enough layers, default to zero - nzm = kc+1 - do K=1,nzm - CS%w_strct(i,j,K) = 0.0 - CS%u_strct(i,j,K) = 0.0 - CS%W_profile(i,j,K) = 0.0 - CS%Uavg_profile(i,j,K) = 0.0 - CS%z_depths(i,j,K) = 0.0 ! could use actual values - CS%N2(i,j,K) = 0.0 ! could use with actual values - enddo - CS%num_intfaces(i,j) = nzm - endif ! kc >= 3 and kc > ModeNum + 1? - endif ! drxh_sum >= 0? - !else ! if at test point - delete later - ! return ! if at test point - delete later - !endif ! if at test point - delete later - endif ! mask2dT > 0.0? - else - ! if cn=0.0, default to zero - nzm = nz+1 ! could use actual values - do K=1,nzm - CS%w_strct(i,j,K) = 0.0 - CS%u_strct(i,j,K) = 0.0 - CS%W_profile(i,j,K) = 0.0 - CS%Uavg_profile(i,j,K) = 0.0 - CS%z_depths(i,j,K) = 0.0 ! could use actual values - CS%N2(i,j,K) = 0.0 ! could use with actual values - enddo - CS%num_intfaces(i,j) = nzm - endif ; enddo ! if cn>0.0? ; i-loop - enddo ! j-loop - - if (CS%debug) call hchksum(CS%N2, 'N2 in wave_struct', G%HI, scale=US%s_to_T**2) - if (CS%debug) call hchksum(cn, 'cn in wave_struct', G%HI, scale=US%L_T_to_m_s) - if (CS%debug) call hchksum(CS%W_profile, 'Wprofile in wave_struct', G%HI, scale=US%Z_to_L*US%L_T_to_m_s) - if (CS%debug) call hchksum(CS%Uavg_profile, 'Uavg_profile in wave_struct', G%HI, scale=US%L_T_to_m_s) - -end subroutine wave_structure - -! The subroutine tridiag_solver is never used and could perhaps be deleted. - -!> Solves a tri-diagonal system Ax=y using either the standard -!! Thomas algorithm (TDMA_T) or its more stable variant that invokes the -!! "Hallberg substitution" (TDMA_H). -subroutine tridiag_solver(a, b, c, h, y, method, x) - real, dimension(:), intent(in) :: a !< lower diagonal with first entry equal to zero. - real, dimension(:), intent(in) :: b !< middle diagonal. - real, dimension(:), intent(in) :: c !< upper diagonal with last entry equal to zero. - real, dimension(:), intent(in) :: h !< vector of values that have already been added to b; used - !! for systems of the form (e.g. average layer thickness in vertical diffusion case): - !! [ -alpha(k-1/2) ] * e(k-1) + - !! [ alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k) + - !! [ -alpha(k+1/2) ] * e(k+1) = y(k) - !! where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)], - !! and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method. - real, dimension(:), intent(in) :: y !< vector of known values on right hand side. - character(len=*), intent(in) :: method !< A string describing the algorithm to use - real, dimension(:), intent(out) :: x !< vector of unknown values to solve for. - ! Local variables - integer :: nrow ! number of rows in A matrix -! real, allocatable, dimension(:,:) :: A_check ! for solution checking -! real, allocatable, dimension(:) :: y_check ! for solution checking - real, allocatable, dimension(:) :: c_prime, y_prime, q, alpha - ! intermediate values for solvers - real :: Q_prime, beta ! intermediate values for solver - integer :: k ! row (e.g. interface) index - - nrow = size(y) - allocate(c_prime(nrow)) - allocate(y_prime(nrow)) - allocate(q(nrow)) - allocate(alpha(nrow)) -! allocate(A_check(nrow,nrow)) -! allocate(y_check(nrow)) - - if (method == 'TDMA_T') then - ! Standard Thomas algoritim (4th variant). - ! Note: Requires A to be non-singular for accuracy/stability - c_prime(:) = 0.0 ; y_prime(:) = 0.0 - c_prime(1) = c(1)/b(1) ; y_prime(1) = y(1)/b(1) - - ! Forward sweep - do k=2,nrow-1 - c_prime(k) = c(k)/(b(k)-a(k)*c_prime(k-1)) - enddo - !print *, 'c_prime=', c_prime(1:nrow) - do k=2,nrow - y_prime(k) = (y(k)-a(k)*y_prime(k-1))/(b(k)-a(k)*c_prime(k-1)) - enddo - !print *, 'y_prime=', y_prime(1:nrow) - x(nrow) = y_prime(nrow) - - ! Backward sweep - do k=nrow-1,1,-1 - x(k) = y_prime(k)-c_prime(k)*x(k+1) - enddo - !print *, 'x=',x(1:nrow) - - ! Check results - delete later - !do j=1,nrow ; do i=1,nrow - ! if (i==j)then ; A_check(i,j) = b(i) - ! elseif (i==j+1)then ; A_check(i,j) = a(i) - ! elseif (i==j-1)then ; A_check(i,j) = c(i) - ! endif - !enddo ; enddo - !print *, 'A(2,1),A(2,2),A(1,2)=', A_check(2,1), A_check(2,2), A_check(1,2) - !y_check = matmul(A_check,x) - !if (all(y_check /= y))then - ! print *, "tridiag_solver: Uh oh, something's not right!" - ! print *, "y=", y - ! print *, "y_check=", y_check - !endif - - elseif (method == 'TDMA_H') then - ! Thomas algoritim (4th variant) w/ Hallberg substitution. - ! For a layered system where k is at interfaces, alpha{k+1/2} refers to - ! some property (e.g. inverse thickness for mode-structure problem) of the - ! layer below and alpha{k-1/2} refers to the layer above. - ! Here, alpha(k)=alpha{k+1/2} and alpha(k-1)=alpha{k-1/2}. - ! Strictly speaking, this formulation requires A to be a non-singular, - ! symmetric, diagonally dominant matrix, with h>0. - ! Need to add a check for these conditions. - do k=1,nrow-1 - if (abs(a(k+1)-c(k)) > 1.e-10*(abs(a(k+1))+abs(c(k)))) then - call MOM_error(FATAL, "tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H") - endif - enddo - alpha = -c - ! Alpha of the bottom-most layer is not necessarily zero. Therefore, - ! back out the value from the provided b(nrow and h(nrow) values - alpha(nrow) = b(nrow)-h(nrow)-alpha(nrow-1) - ! Prime other variables - beta = 1/b(1) - y_prime(:) = 0.0 ; q(:) = 0.0 - y_prime(1) = beta*y(1) ; q(1) = beta*alpha(1) - Q_prime = 1-q(1) - - ! Forward sweep - do k=2,nrow-1 - beta = 1/(h(k)+alpha(k-1)*Q_prime+alpha(k)) - if (isnan(beta))then ; print *, "Tridiag_solver: beta is NAN" ; endif - q(k) = beta*alpha(k) - y_prime(k) = beta*(y(k)+alpha(k-1)*y_prime(k-1)) - Q_prime = beta*(h(k)+alpha(k-1)*Q_prime) - enddo - if ((h(nrow)+alpha(nrow-1)*Q_prime+alpha(nrow)) == 0.0)then - call MOM_error(FATAL, "Tridiag_solver: this system is not stable.") ! ; overriding beta(nrow) - ! This has hard-coded dimensions: beta = 1/(1e-15) ! place holder for unstable systems - delete later - else - beta = 1/(h(nrow)+alpha(nrow-1)*Q_prime+alpha(nrow)) - endif - y_prime(nrow) = beta*(y(nrow)+alpha(nrow-1)*y_prime(nrow-1)) - x(nrow) = y_prime(nrow) - ! Backward sweep - do k=nrow-1,1,-1 - x(k) = y_prime(k)+q(k)*x(k+1) - enddo - !print *, 'yprime=',y_prime(1:nrow) - !print *, 'x=',x(1:nrow) - endif - - deallocate(c_prime,y_prime,q,alpha) -! deallocate(A_check,y_check) - -end subroutine tridiag_solver - -!> Allocate memory associated with the wave structure module and read parameters. -subroutine wave_structure_init(Time, G, GV, param_file, diag, CS) - type(time_type), intent(in) :: Time !< The current model time. - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time - !! parameters. - type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate - !! diagnostic output. - type(wave_structure_CS), intent(inout) :: CS !< Wave structure control struct - - ! This include declares and sets the variable "version". -# include "version_variable.h" - character(len=40) :: mdl = "MOM_wave_structure" ! This module's name. - integer :: isd, ied, jsd, jed, nz - - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke - - CS%initialized = .true. - - ! call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_TEST", CS%int_tide_source_test, & - ! "If true, apply an arbitrary generation site for internal tide testing", & - ! default=.false.) - ! if (CS%int_tide_source_test) then - ! call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_I", CS%int_tide_source_i, & - ! "I Location of generation site for internal tide", default=0) - ! call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_J", CS%int_tide_source_j, & - ! "J Location of generation site for internal tide", default=0) - ! endif - call get_param(param_file, mdl, "DEBUG", CS%debug, & - "debugging prints", default=.false.) - - CS%diag => diag - - ! Allocate memory for variable in control structure; note, - ! not all rows will be filled if layers get merged! - allocate(CS%w_strct(isd:ied,jsd:jed,nz+1)) - allocate(CS%u_strct(isd:ied,jsd:jed,nz+1)) - allocate(CS%W_profile(isd:ied,jsd:jed,nz+1)) - allocate(CS%Uavg_profile(isd:ied,jsd:jed,nz+1)) - allocate(CS%z_depths(isd:ied,jsd:jed,nz+1)) - allocate(CS%N2(isd:ied,jsd:jed,nz+1)) - allocate(CS%num_intfaces(isd:ied,jsd:jed)) - - ! Write all relevant parameters to the model log. - call log_version(param_file, mdl, version, "") - -end subroutine wave_structure_init - -end module MOM_wave_structure diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index d3f202339a..ec07939ee4 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -23,7 +23,6 @@ module MOM_internal_tides use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_wave_structure, only: wave_structure_init, wave_structure, wave_structure_CS implicit none ; private From d0f7b297be08eb07fd2781fb03c858a39e648f9c Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 22 May 2023 09:59:36 -0400 Subject: [PATCH 08/15] Autoconf: Better Unicode Python support in makedep MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The `open()` commands in `makedep` for reading Fortran source now includes an `errors=` argument for catching bytes outside of the file character set. Unknown characters are replaced with the "unknown" character (usually �) rather than raising an error. This avoids problems with Unicode characters and older Pythons which do not support them, as well as characters from legacy encodings which can cause errors in Unicode. Substitution does not break any behavior, since Unicode is only permitted inside of comment blocks and strings. This fixes several errors which were silent in `.testing` but were observed by some users which using autoconf to build their own executables. --- ac/makedep | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ac/makedep b/ac/makedep index 439679f17d..225a241b93 100755 --- a/ac/makedep +++ b/ac/makedep @@ -4,9 +4,10 @@ from __future__ import print_function import argparse import glob +import io import os import re -import sys # used only to get path to current script +import sys # Pre-compile re searches @@ -255,7 +256,7 @@ def scan_fortran_file(src_file): """Scan the Fortran file "src_file" and return lists of module defined, module used, and files included.""" module_decl, used_modules, cpp_includes, f90_includes, programs = [], [], [], [], [] - with open(src_file, 'r') as file: + with io.open(src_file, 'r', errors='replace') as file: lines = file.readlines() for line in lines: match = re_module.match(line.lower()) From b075794d660a71693a67a0ee1de85fc1c24059ed Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 22 May 2023 14:54:43 -0400 Subject: [PATCH 09/15] Autoconf: Fix Python test and allow configuration The AC_PATH_PROGS macros used in Python testing were incorrectly using AC_MSG_ERROR in places where a missing value for PYTHON should be if the executable was not found. It also did not permit for a configurable PYTHON variable, since the autodetect was always run, even if PYTHON were set. This has been updated so that Python autodetection only runs if PYTHON is unset. It also correctly reports a failed configuration if PYTHON is not found. (It does not, however, test of PYTHON is actually a Python interpreter, but we can deal with that at a later date.) --- ac/configure.ac | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/ac/configure.ac b/ac/configure.ac index 1c10c14495..7ea1870816 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -236,15 +236,21 @@ AC_COMPILE_IFELSE( ] ) +# Python interpreter test -# Verify that Python is available -AC_PATH_PROGS([PYTHON], [python python3 python2], [ - AC_MSG_ERROR([Could not find python.]) -]) AC_ARG_VAR([PYTHON], [Python interpreter command]) +AS_VAR_SET_IF([PYTHON], [ + AC_PATH_PROGS([PYTHON], ["$PYTHON"], [none]) +], [ + AC_PATH_PROGS([PYTHON], [python python3 python2], [none]) +]) +AS_VAR_IF([PYTHON], [none], [ + AC_MSG_ERROR([Python interpreter not found.]) +]) + -# Verify that makedep is available +# Makedep test AC_PATH_PROG([MAKEDEP], [makedep], [${srcdir}/ac/makedep]) AC_SUBST([MAKEDEP]) From cb4574b3f540bd9e6bdfced6e04a56c8cbd6c3bb Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 29 Mar 2023 10:09:37 -0400 Subject: [PATCH 10/15] Fix PGI runtime issue with class(*) - Some tests such as global_ALE_z crash under PGI (ncrc4.pgi20 or ncrc5.pgi227) with FATAL from PE 27: unsupported attribute type: get_variable_attribute_0d: file:INPUT/tideamp.nc- variable:GRID_X_T attribute: axis - PGI in general has issues with class(*) construct and in this case cannot recognize the axis argument to be a string. - This mod helps PGI recognize that the argument is a string. --- config_src/infra/FMS2/MOM_io_infra.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 8802761774..2c3a5b8ad3 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -1524,9 +1524,9 @@ subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t if (variable_exists(fileobj, trim(dim_names(i)))) then cartesian = "" if (variable_att_exists(fileobj, trim(dim_names(i)), "cartesian_axis")) then - call get_variable_attribute(fileobj, trim(dim_names(i)), "cartesian_axis", cartesian) + call get_variable_attribute(fileobj, trim(dim_names(i)), "cartesian_axis", cartesian(1:1)) elseif (variable_att_exists(fileobj, trim(dim_names(i)), "axis")) then - call get_variable_attribute(fileobj, trim(dim_names(i)), "axis", cartesian) + call get_variable_attribute(fileobj, trim(dim_names(i)), "axis", cartesian(1:1)) endif cartesian = adjustl(cartesian) if ((index(cartesian, "X") == 1) .or. (index(cartesian, "x") == 1)) is_x(i) = .true. From 273da2fb272e37860c9617be3857aa8236a34f66 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 9 Jun 2023 18:17:59 -0400 Subject: [PATCH 11/15] Use fileset rather than threading for decompositon MOM IO was using the `threading` flag rather than `fileset` to determine whether a file should be forced as single file rather than domain-decomposed. This patch applies the correct flag. --- src/framework/MOM_io.F90 | 8 +++++--- src/framework/MOM_io_file.F90 | 4 ++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 6bde678eb4..bebce6f502 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -338,12 +338,14 @@ subroutine create_MOM_file(IO_handle, filename, vars, novars, fields, & if (one_file) then if (domain_set) then call IO_handle%open(filename, action=OVERWRITE_FILE, & - MOM_domain=domain, threading=thread) + MOM_domain=domain, threading=thread, fileset=SINGLE_FILE) else - call IO_handle%open(filename, action=OVERWRITE_FILE, threading=thread) + call IO_handle%open(filename, action=OVERWRITE_FILE, threading=thread, & + fileset=SINGLE_FILE) endif else - call IO_handle%open(filename, action=OVERWRITE_FILE, MOM_domain=Domain) + call IO_handle%open(filename, action=OVERWRITE_FILE, MOM_domain=Domain, & + threading=thread, fileset=thread) endif ! Define the coordinates. diff --git a/src/framework/MOM_io_file.F90 b/src/framework/MOM_io_file.F90 index 6909e597ba..6eaa10f622 100644 --- a/src/framework/MOM_io_file.F90 +++ b/src/framework/MOM_io_file.F90 @@ -929,8 +929,8 @@ subroutine open_file_infra(handle, filename, action, MOM_domain, threading, file ! True if the domain is replaced with a single-file IO layout. use_single_file_domain = .false. - if (present(MOM_domain) .and. present(threading)) then - if (threading == SINGLE_FILE) & + if (present(MOM_domain) .and. present(fileset)) then + if (fileset == SINGLE_FILE) & use_single_file_domain = .true. endif From 6038735c43a08a18f3000ddf61f27821251f07b1 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 31 May 2023 14:05:39 -0400 Subject: [PATCH 12/15] FMS2 interpolation ID replaced with derived type All instances of an FMS ID to the internal interpolation content is replaced with a derived type containing additional metadata recording the field's origin filename and fieldname. This additional information is required in order to replicate the axis data from the field, which is no longer provided by FMS2. The abstraction of this type also allows us to either extend it or redefine it in other frameworks as needed in the future. This primarily affects the usage of the following functions: - init_external_field - time_interp_external - horiz_interp_and_extrap_tracer The following solvers are updated: - MOM_open_boundary - MOM_ice_shelf - MOM_oda_driver - MOM_MEKE - MOM_ALE_sponge - MOM_diabatic_aux Of these, OBC was the most significant. The integer handle (fid) was previously used to determine if each segment field was constant or (if negative) read from a file. After being replaced by the derived type, a new flag was added to make this determination. All of the coupled drivers have been modified, since they support time interpolation of T and S fields. - FMS - MCT - NUOPC The NUOPC driver also includes modifications to its CFC11 and CFC12 fields. Changes to the MOM CFC modules replaces an `id == -1`-like test, which is not used by the derived type. This check has been removed, and we now solely rely on the `present(cfc_handle)` test. While this could change behavior, there does not seem to be any scenario where init_external_field would return -1 but would be passed to the function. (But I may eat these words.) --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 15 +++-- .../mct_cap/mom_surface_forcing_mct.F90 | 15 +++-- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 28 +++++---- config_src/infra/FMS1/MOM_interp_infra.F90 | 54 +++++++++++------- config_src/infra/FMS2/MOM_interp_infra.F90 | 57 +++++++++++-------- src/core/MOM_open_boundary.F90 | 31 +++++----- src/framework/MOM_horizontal_regridding.F90 | 11 ++-- src/framework/MOM_interpolate.F90 | 27 +++++---- src/ice_shelf/MOM_ice_shelf.F90 | 17 +++--- src/ocean_data_assim/MOM_oda_driver.F90 | 15 ++--- src/parameterizations/lateral/MOM_MEKE.F90 | 7 ++- .../vertical/MOM_ALE_sponge.F90 | 32 +++++------ .../vertical/MOM_diabatic_aux.F90 | 3 +- src/tracer/MOM_CFC_cap.F90 | 20 ++++--- 14 files changed, 187 insertions(+), 145 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index 26ab6269ef..f70cd34012 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -27,6 +27,7 @@ module MOM_surface_forcing_gfdl use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : read_netCDF_data use MOM_io, only : stdout_if_root @@ -153,8 +154,10 @@ module MOM_surface_forcing_gfdl !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' real, pointer, dimension(:,:) :: trestore_mask => NULL() !< Mask for SST restoring [nondim] - integer :: id_srestore = -1 !< An id number for time_interp_external. - integer :: id_trestore = -1 !< An id number for time_interp_external. + type(external_field) :: srestore_handle + !< Handle for time-interpolated salt restoration field + type(external_field) :: trestore_handle + !< Handle for time-interpolated temperature restoration field type(forcing_diags), public :: handles !< Diagnostics handles @@ -345,7 +348,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (CS%restore_salt) then - call time_interp_external(CS%id_srestore, Time, data_restore, scale=US%ppt_to_S) + call time_interp_external(CS%srestore_handle, Time, data_restore, scale=US%ppt_to_S) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice @@ -403,7 +406,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (CS%restore_temp) then - call time_interp_external(CS%id_trestore, Time, data_restore, scale=US%degC_to_C) + call time_interp_external(CS%trestore_handle, Time, data_restore, scale=US%degC_to_C) if ( CS%trestore_SPEAR_ECDA ) then do j=js,je ; do i=is,ie if (abs(data_restore(i,j)+1.8*US%degC_to_C) < 0.0001*US%degC_to_C) then @@ -1610,7 +1613,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if (CS%restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) - CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, MOM_domain=G%Domain) + CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, MOM_domain=G%Domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 if (CS%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' @@ -1620,7 +1623,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if (CS%restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) - CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, MOM_domain=G%Domain) + CS%trestore_handle = init_external_field(temp_file, CS%temp_restore_var_name, MOM_domain=G%Domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 if (CS%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index 0364d46ddc..9b858af94e 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 @@ -25,6 +25,7 @@ module MOM_surface_forcing_mct use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : stdout use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS @@ -134,8 +135,10 @@ module MOM_surface_forcing_mct !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. + type(external_field) :: srestore_handle + !< Handle for time-interpolated salt restoration field + type(external_field) :: trestore_handle + !< Handle for time-interpolated temperature restoration field type(forcing_diags), public :: handles !< diagnostics handles type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< restart pointer @@ -348,7 +351,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (restore_salinity) then - call time_interp_external(CS%id_srestore, Time, data_restore, scale=US%ppt_to_S) + call time_interp_external(CS%srestore_handle, Time, data_restore, scale=US%ppt_to_S) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice @@ -405,7 +408,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (restore_sst) then - call time_interp_external(CS%id_trestore, Time, data_restore, scale=US%degC_to_C) + call time_interp_external(CS%trestore_handle, Time, data_restore, scale=US%degC_to_C) do j=js,je ; do i=is,ie delta_sst = data_restore(i,j) - sfc_state%SST(i,j) delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) @@ -1292,7 +1295,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_salt)) then ; if (restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) - CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) + CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 if (CS%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' @@ -1302,7 +1305,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_temp)) then ; if (restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) - CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) + CS%trestore_handle = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 if (CS%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 2c8e3db8bd..b8162b4b59 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -26,6 +26,7 @@ module MOM_surface_forcing_nuopc use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_CFC_cap, only : CFC_cap_fluxes use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : stdout @@ -146,10 +147,14 @@ module MOM_surface_forcing_nuopc character(len=30) :: cfc11_var_name !< name of cfc11 in CFC_BC_file character(len=30) :: cfc12_var_name !< name of cfc11 in CFC_BC_file real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. - integer :: id_cfc11_atm = -1 !< id number for time_interp_external. - integer :: id_cfc12_atm = -1 !< id number for time_interp_external. + type(external_field) :: srestore_handle + !< Handle for time-interpolated salt restoration field + type(external_field) :: trestore_handle + !< Handle for time-interpolated temperature restoration field + type(external_field) :: cfc11_atm_handle + !< Handle for time-interpolated CFC11 restoration field + type(external_field) :: cfc12_atm_handle + !< Handle for time-interpolated CFC12 restoration field ! Diagnostics handles type(forcing_diags), public :: handles @@ -377,7 +382,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (restore_salinity) then - call time_interp_external(CS%id_srestore, Time, data_restore, scale=US%ppt_to_S) + call time_interp_external(CS%srestore_handle, Time, data_restore, scale=US%ppt_to_S) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice @@ -434,7 +439,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (restore_sst) then - call time_interp_external(CS%id_trestore, Time, data_restore, scale=US%degC_to_C) + call time_interp_external(CS%trestore_handle, Time, data_restore, scale=US%degC_to_C) do j=js,je ; do i=is,ie delta_sst = data_restore(i,j) - sfc_state%SST(i,j) delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) @@ -596,7 +601,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! CFCs if (CS%use_CFC) then - call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) + call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, & + CS%cfc11_atm_handle, CS%cfc11_atm_handle) endif if (associated(IOB%salt_flux)) then @@ -1394,7 +1400,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_salt)) then ; if (restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) - CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) + CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 if (CS%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' @@ -1404,7 +1410,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_temp)) then ; if (restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) - CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) + CS%trestore_handle = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 if (CS%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' @@ -1430,8 +1436,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "The name of the variable representing CFC-12 in "//& "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.) - CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) - CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) + CS%cfc11_atm_handle = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) + CS%cfc12_atm_handle = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) endif endif diff --git a/config_src/infra/FMS1/MOM_interp_infra.F90 b/config_src/infra/FMS1/MOM_interp_infra.F90 index 224e26a051..e14233a64b 100644 --- a/config_src/infra/FMS1/MOM_interp_infra.F90 +++ b/config_src/infra/FMS1/MOM_interp_infra.F90 @@ -18,6 +18,18 @@ module MOM_interp_infra public :: time_interp_extern, init_extern_field, time_interp_extern_init public :: get_external_field_info, axistype, get_axis_data public :: run_horiz_interp, build_horiz_interp_weights +public :: external_field + +!< Handle of an external field for interpolation +type :: external_field + private + integer :: id + !< FMS ID for the interpolated field + character(len=:), allocatable :: filename + !< Filename containing the field values + character(len=:), allocatable :: label + !< Field name in the file +end type external_field !> Read a field based on model time, and rotate to the model domain. interface time_interp_extern @@ -167,8 +179,8 @@ end function get_extern_field_missing !> Get information about the external fields. -subroutine get_external_field_info(field_id, size, axes, missing) - integer, intent(in) :: field_id !< The integer index of the external +subroutine get_external_field_info(field, size, axes, missing) + type(external_field), intent(in) :: field !< Handle for time interpolated external !! field returned from a previous !! call to init_external_field() integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data @@ -176,37 +188,35 @@ subroutine get_external_field_info(field_id, size, axes, missing) real, optional, intent(inout) :: missing !< Missing value for the input data if (present(size)) then - size(1:4) = get_extern_field_size(field_id) + size(1:4) = get_extern_field_size(field%id) endif if (present(axes)) then - axes(1:4) = get_extern_field_axes(field_id) + axes(1:4) = get_extern_field_axes(field%id) endif if (present(missing)) then - missing = get_extern_field_missing(field_id) + missing = get_extern_field_missing(field%id) endif end subroutine get_external_field_info !> Read a scalar field based on model time. -subroutine time_interp_extern_0d(field_id, time, data_in, verbose) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_0d(field, time, data_in, verbose) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, intent(inout) :: data_in !< The interpolated value logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging - call time_interp_external(field_id, time, data_in, verbose=verbose) + call time_interp_external(field%id, time, data_in, verbose=verbose) end subroutine time_interp_extern_0d !> Read a 2d field from an external based on model time, potentially including horizontal !! interpolation and rotation of the data -subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_2d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated values integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method @@ -216,15 +226,14 @@ subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_2d !> Read a 3d field based on model time, and rotate to the model grid -subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_3d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated values integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method @@ -234,14 +243,15 @@ subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_3d !> initialize an external field -integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) +function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency) & + result(field) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -261,17 +271,17 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, !! is in use, and (2) the modulo time period of the !! data is an integer number of years, then map !! a model date of Feb 29. onto a common year on Feb. 28. + type(external_field) :: field !< Handle to external field if (present(MOM_Domain)) then - init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & + field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) else - init_extern_field = init_external_field(file, fieldname, domain=domain, & + field%id = init_external_field(file, fieldname, domain=domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif - end function init_extern_field end module MOM_interp_infra diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index c29459aad1..7964b3537f 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -18,6 +18,18 @@ module MOM_interp_infra public :: time_interp_extern, init_extern_field, time_interp_extern_init public :: get_external_field_info, axistype, get_axis_data public :: run_horiz_interp, build_horiz_interp_weights +public :: external_field + +!< Handle of an external field for interpolation +type :: external_field + private + integer :: id + !< FMS ID for the interpolated field + character(len=:), allocatable :: filename + !< Filename containing the field values + character(len=:), allocatable :: label + !< Field name in the file +end type external_field !> Read a field based on model time, and rotate to the model domain. interface time_interp_extern @@ -166,8 +178,8 @@ end function get_extern_field_missing !> Get information about the external fields. -subroutine get_external_field_info(field_id, size, axes, missing) - integer, intent(in) :: field_id !< The integer index of the external +subroutine get_external_field_info(field, size, axes, missing) + type(external_field), intent(in) :: field !< Handle for time interpolated external !! field returned from a previous !! call to init_external_field() integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data @@ -175,37 +187,35 @@ subroutine get_external_field_info(field_id, size, axes, missing) real, optional, intent(inout) :: missing !< Missing value for the input data if (present(size)) then - size(1:4) = get_extern_field_size(field_id) + size(1:4) = get_extern_field_size(field%id) endif if (present(axes)) then - axes(1:4) = get_extern_field_axes(field_id) + axes(1:4) = get_extern_field_axes(field%id) endif if (present(missing)) then - missing = get_extern_field_missing(field_id) + missing = get_extern_field_missing(field%id) endif end subroutine get_external_field_info !> Read a scalar field based on model time. -subroutine time_interp_extern_0d(field_id, time, data_in, verbose) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_0d(field, time, data_in, verbose) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, intent(inout) :: data_in !< The interpolated value logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging - call time_interp_external(field_id, time, data_in, verbose=verbose) + call time_interp_external(field%id, time, data_in, verbose=verbose) end subroutine time_interp_extern_0d !> Read a 2d field from an external based on model time, potentially including horizontal !! interpolation and rotation of the data -subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_2d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated values integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method @@ -215,15 +225,14 @@ subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_2d !> Read a 3d field based on model time, and rotate to the model grid -subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_3d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated values integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method @@ -233,14 +242,15 @@ subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_3d !> initialize an external field -integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) +function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency) & + result(field) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -260,19 +270,20 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, !! is in use, and (2) the modulo time period of the !! data is an integer number of years, then map !! a model date of Feb 29. onto a common year on Feb. 28. + type(external_field) :: field !< Handle to external field - + field%filename = file + field%label = fieldname if (present(MOM_Domain)) then - init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & + field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) else - init_extern_field = init_external_field(file, fieldname, domain=domain, & + field%id = init_external_field(file, fieldname, domain=domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif - end function init_extern_field end module MOM_interp_infra diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9bd292e796..ba8b8ce818 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -24,6 +24,7 @@ module MOM_open_boundary use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_regridding, only : regridding_CS @@ -81,8 +82,9 @@ module MOM_open_boundary !> Open boundary segment data from files (mostly). type, public :: OBC_segment_data_type - integer :: fid !< handle from FMS associated with segment data on disk - integer :: fid_dz !< handle from FMS associated with segment thicknesses on disk + type(external_field) :: handle !< handle from FMS associated with segment data on disk + type(external_field) :: dz_handle !< handle from FMS associated with segment thicknesses on disk + logical :: use_IO = .false. !< True if segment data is based on file input character(len=32) :: name !< a name identifier for the segment data character(len=8) :: genre !< an identifier for the segment data real :: scale !< A scaling factor for converting input data to @@ -96,7 +98,7 @@ module MOM_open_boundary real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid. !! The values for tracers should have the same units as the field !! they are being applied to? - real :: value !< constant value if fid is equal to -1 + real :: value !< constant value if not read from file real :: resrv_lfac_in = 1. !< reservoir inverse length scale factor for IN direction per field !< the general 1/Lscale_IN is multiplied by this factor for each tracer real :: resrv_lfac_out= 1. !< reservoir inverse length scale factor for OUT direction per field @@ -842,6 +844,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) ! The scale factor for tracers may also be set in register_segment_tracer, and a constant input ! value is rescaled there. segment%field(m)%scale = scale_factor_from_name(fields(m), GV, US, segment%tr_Reg) + segment%field(m)%use_IO = .true. if (segment%field(m)%name == 'TEMP') then segment%temp_segment_data_exists = .true. segment%t_values_needed = .false. @@ -957,7 +960,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) endif endif segment%field(m)%buffer_src(:,:,:) = 0.0 - segment%field(m)%fid = init_external_field(trim(filename), trim(fieldname), & + segment%field(m)%handle = init_external_field(trim(filename), trim(fieldname), & ignore_axis_atts=.true., threading=SINGLE_FILE) if (siz(3) > 1) then if ((index(segment%field(m)%name, 'phase') > 0) .or. (index(segment%field(m)%name, 'amp') > 0)) then @@ -988,7 +991,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) endif segment%field(m)%dz_src(:,:,:)=0.0 segment%field(m)%nk_src=siz(3) - segment%field(m)%fid_dz = init_external_field(trim(filename), trim(fieldname), & + segment%field(m)%dz_handle = init_external_field(trim(filename), trim(fieldname), & ignore_axis_atts=.true., threading=SINGLE_FILE) endif else @@ -996,12 +999,12 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) endif endif else - segment%field(m)%fid = -1 segment%field(m)%name = trim(fields(m)) ! The scale factor for tracers may also be set in register_segment_tracer, and a constant input ! value is rescaled there. segment%field(m)%scale = scale_factor_from_name(fields(m), GV, US, segment%tr_Reg) segment%field(m)%value = segment%field(m)%scale * value + segment%field(m)%use_IO = .false. ! Check if this is a tidal field. If so, the number ! of expected constituents must be 1. @@ -3892,7 +3895,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) !a less frequent update as set by the parameter update_OBC_period_max in MOM.F90. !Cycle if it is not the time to update OBC segment data for this field. if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle - if (segment%field(m)%fid > 0) then + if (segment%field(m)%use_IO) then siz(1)=size(segment%field(m)%buffer_src,1) siz(2)=size(segment%field(m)%buffer_src,2) siz(3)=size(segment%field(m)%buffer_src,3) @@ -3972,7 +3975,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif ! This is where the data values are actually read in. - call time_interp_external(segment%field(m)%fid, Time, tmp_buffer_in, scale=segment%field(m)%scale) + call time_interp_external(segment%field(m)%handle, Time, tmp_buffer_in, scale=segment%field(m)%scale) ! NOTE: Rotation of face-points require that we skip the final value if (turns /= 0) then @@ -4045,7 +4048,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (segment%field(m)%nk_src > 1 .and.& (index(segment%field(m)%name, 'phase') <= 0 .and. index(segment%field(m)%name, 'amp') <= 0)) then ! This is where the 2-d tidal data values are actually read in. - call time_interp_external(segment%field(m)%fid_dz, Time, tmp_buffer_in, scale=US%m_to_Z) + call time_interp_external(segment%field(m)%dz_handle, Time, tmp_buffer_in, scale=US%m_to_Z) if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. if (segment%is_E_or_W & @@ -4211,7 +4214,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) deallocate(tmp_buffer) if (turns /= 0) & deallocate(tmp_buffer_in) - else ! fid <= 0 (Uniform value) + else ! use_IO = .false. (Uniform value) if (.not. allocated(segment%field(m)%buffer_dst)) then if (segment%is_E_or_W) then if (segment%field(m)%name == 'V') then @@ -4257,7 +4260,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) do m = 1,segment%num_fields !cycle if it is not the time to update OBGC tracers from source if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle - ! if (segment%field(m)%fid>0) then + ! if (segment%field(m)%use_IO) then ! calculate external BT velocity and transport if needed if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then if (trim(segment%field(m)%name) == 'U' .and. segment%is_E_or_W) then @@ -4684,7 +4687,7 @@ subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, & ! rescale the previously stored input values. Note that calls to register_segment_tracer ! can come before or after calls to initialize_segment_data. if (uppercase(segment%field(m)%name) == uppercase(segment%tr_Reg%Tr(ntseg)%name)) then - if (segment%field(m)%fid == -1) then + if (.not. segment%field(m)%use_IO) then rescale = scale if ((segment%field(m)%scale /= 0.0) .and. (segment%field(m)%scale /= 1.0)) & rescale = scale / segment%field(m)%scale @@ -5948,8 +5951,8 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%num_fields = segment_in%num_fields do n = 1, num_fields - segment%field(n)%fid = segment_in%field(n)%fid - segment%field(n)%fid_dz = segment_in%field(n)%fid_dz + segment%field(n)%handle = segment_in%field(n)%handle + segment%field(n)%dz_handle = segment_in%field(n)%dz_handle if (modulo(turns, 2) /= 0) then select case (segment_in%field(n)%name) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 83e7718311..bedf710582 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -17,6 +17,7 @@ module MOM_horizontal_regridding use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights use MOM_interp_infra, only : horiz_interp_type, horizontal_interp_init use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data +use MOM_interp_infra, only : external_field use MOM_time_manager, only : time_type use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data use MOM_io, only : read_attribute, read_variable @@ -598,12 +599,12 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr end subroutine horiz_interp_and_extrap_tracer_record !> Extrapolate and interpolate using a FMS time interpolation handle -subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, & +subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & z_in, z_edges_in, missing_value, scale, & homogenize, spongeOngrid, m_to_Z, & answers_2018, tr_iter_tol, answer_date) - integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator + type(external_field), intent(in) :: field !< Handle for the time interpolated field type(time_type), intent(in) :: Time !< A FMS time type type(ocean_grid_type), intent(inout) :: G !< Grid object real, allocatable, dimension(:,:,:), intent(out) :: tr_z @@ -716,7 +717,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, call cpu_clock_begin(id_clock_read) - call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_val_in) + call get_external_field_info(field, size=fld_sz, axes=axes_data, missing=missing_val_in) missing_value = scale*missing_val_in verbosity = MOM_get_verbosity() @@ -790,7 +791,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, if (.not.is_ongrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>5), turns=turns) + call time_interp_external(field, Time, data_in, verbose=(verbosity>5), turns=turns) ! Loop through each data level and interpolate to model grid. ! After interpolating, fill in points which will be needed to define the layers. @@ -897,7 +898,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>5), turns=turns) + call time_interp_external(field, Time, data_in, verbose=(verbosity>5), turns=turns) do k=1,kd do j=js,je do i=is,ie diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 index 38a786e593..e131e8db9d 100644 --- a/src/framework/MOM_interpolate.F90 +++ b/src/framework/MOM_interpolate.F90 @@ -9,12 +9,14 @@ module MOM_interpolate use MOM_interp_infra, only : time_interp_external_init=>time_interp_extern_init use MOM_interp_infra, only : horiz_interp_type, get_external_field_info use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights +use MOM_interp_infra, only : external_field use MOM_time_manager, only : time_type implicit none ; private public :: time_interp_external, init_external_field, time_interp_external_init, get_external_field_info public :: horiz_interp_type, run_horiz_interp, build_horiz_interp_weights +public :: external_field !> Read a field based on model time, and rotate to the model domain. interface time_interp_external @@ -26,9 +28,8 @@ module MOM_interpolate contains !> Read a scalar field based on model time. -subroutine time_interp_external_0d(field_id, time, data_in, verbose, scale) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_external_0d(field, time, data_in, verbose, scale) + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, intent(inout) :: data_in !< The interpolated value logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging @@ -48,7 +49,7 @@ subroutine time_interp_external_0d(field_id, time, data_in, verbose, scale) data_in = data_in * I_scale endif ; endif - call time_interp_extern(field_id, time, data_in, verbose=verbose) + call time_interp_extern(field, time, data_in, verbose=verbose) if (present(scale)) then ; if (scale /= 1.0) then ! Rescale data that has been newly set and restore the scaling of unset data. @@ -63,10 +64,9 @@ end subroutine time_interp_external_0d !> Read a 2d field from an external based on model time, potentially including horizontal !! interpolation and rotation of the data -subroutine time_interp_external_2d(field_id, time, data_in, interp, & +subroutine time_interp_external_2d(field, time, data_in, interp, & verbose, horz_interp, mask_out, turns, scale) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated values integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method @@ -105,11 +105,11 @@ subroutine time_interp_external_2d(field_id, time, data_in, interp, & qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) if (qturns == 0) then - call time_interp_extern(field_id, time, data_in, interp=interp, & + call time_interp_extern(field, time, data_in, interp=interp, & verbose=verbose, horz_interp=horz_interp) else call allocate_rotated_array(data_in, [1,1], -qturns, data_pre_rot) - call time_interp_extern(field_id, time, data_pre_rot, interp=interp, & + call time_interp_extern(field, time, data_pre_rot, interp=interp, & verbose=verbose, horz_interp=horz_interp) call rotate_array(data_pre_rot, turns, data_in) deallocate(data_pre_rot) @@ -136,10 +136,9 @@ end subroutine time_interp_external_2d !> Read a 3d field based on model time, and rotate to the model grid -subroutine time_interp_external_3d(field_id, time, data_in, interp, & +subroutine time_interp_external_3d(field, time, data_in, interp, & verbose, horz_interp, mask_out, turns, scale) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() + type(external_field), intent(in) :: field !< Handle for time interpolated field type(time_type), intent(in) :: time !< The target time for the data real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated values integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method @@ -178,11 +177,11 @@ subroutine time_interp_external_3d(field_id, time, data_in, interp, & qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) if (qturns == 0) then - call time_interp_extern(field_id, time, data_in, interp=interp, & + call time_interp_extern(field, time, data_in, interp=interp, & verbose=verbose, horz_interp=horz_interp) else call allocate_rotated_array(data_in, [1,1,1], -qturns, data_pre_rot) - call time_interp_extern(field_id, time, data_pre_rot, interp=interp, & + call time_interp_extern(field, time, data_pre_rot, interp=interp, & verbose=verbose, horz_interp=horz_interp) call rotate_array(data_pre_rot, turns, data_in) deallocate(data_pre_rot) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 113b6c045b..8e0e58c1b6 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -61,6 +61,7 @@ module MOM_ice_shelf use MOM_spatial_means, only : global_area_integral use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field implicit none ; private @@ -196,10 +197,10 @@ module MOM_ice_shelf id_shelf_sfc_mass_flux = -1 !>@} - integer :: id_read_mass !< An integer handle used in time interpolation of - !! the ice shelf mass read from a file - integer :: id_read_area !< An integer handle used in time interpolation of - !! the ice shelf mass read from a file + type(external_field) :: mass_handle + !< Handle for reading the time interpolated ice shelf mass from a file + type(external_field) :: area_handle + !< Handle for reading the time interpolated ice shelf area from a file type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to control diagnostic output. type(user_ice_shelf_CS), pointer :: user_CS => NULL() !< A pointer to the control structure for @@ -1118,7 +1119,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) do j=js,je ; do i=is,ie last_hmask(i,j) = ISS%hmask(i,j) ; last_area_shelf_h(i,j) = ISS%area_shelf_h(i,j) enddo ; enddo - call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) + call time_interp_external(CS%mass_handle, Time0, last_mass_shelf) do j=js,je ; do i=is,ie ! This should only be done if time_interp_extern did an update. last_mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp @@ -1937,7 +1938,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) filename = trim(slasher(inputdir))//trim(shelf_file) call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename) - CS%id_read_mass = init_external_field(filename, shelf_mass_var, & + CS%mass_handle = init_external_field(filename, shelf_mass_var, & MOM_domain=CS%Grid_in%Domain, verbose=CS%debug) if (read_shelf_area) then @@ -1945,7 +1946,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) "The variable in SHELF_FILE with the shelf area.", & default="shelf_area") - CS%id_read_area = init_external_field(filename, shelf_area_var, & + CS%area_handle = init_external_field(filename, shelf_area_var, & MOM_domain=CS%Grid_in%Domain) endif @@ -2040,7 +2041,7 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) allocate(tmp2d(is:ie,js:je), source=0.0) endif - call time_interp_external(CS%id_read_mass, Time, tmp2d) + call time_interp_external(CS%mass_handle, Time, tmp2d) call rotate_array(tmp2d, CS%turns, ISS%mass_shelf) deallocate(tmp2d) diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 8a1aab3328..53615b0063 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -17,6 +17,7 @@ module MOM_oda_driver_mod use MOM_io, only : SINGLE_FILE use MOM_interp_infra, only : init_extern_field, get_external_field_info use MOM_interp_infra, only : time_interp_extern +use MOM_interpolate, only : external_field use MOM_remapping, only : remappingSchemesDoc use MOM_time_manager, only : time_type, real_to_time, get_date use MOM_time_manager, only : operator(+), operator(>=), operator(/=) @@ -80,8 +81,8 @@ module MOM_oda_driver_mod !> A structure containing integer handles for bias adjustment of tracers type :: INC_CS integer :: fldno = 0 !< The number of tracers - integer :: T_id !< The integer handle for the temperature file - integer :: S_id !< The integer handle for the salinity file + type(external_field) :: T !< The handle for the temperature file + type(external_field) :: S !< The handle for the salinity file end type INC_CS !> Control structure that contains a transpose of the ocean state across ensemble members. @@ -391,11 +392,11 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) "tendency adjustments", default='temp_salt_adjustment.nc') inc_file = trim(inputdir) // trim(bias_correction_file) - CS%INC_CS%T_id = init_extern_field(inc_file, "temp_increment", & + CS%INC_CS%T = init_extern_field(inc_file, "temp_increment", & correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) - CS%INC_CS%S_id = init_extern_field(inc_file, "salt_increment", & + CS%INC_CS%S = init_extern_field(inc_file, "salt_increment", & correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) - call get_external_field_info(CS%INC_CS%T_id,size=fld_sz) + call get_external_field_info(CS%INC_CS%T, size=fld_sz) CS%INC_CS%fldno = 2 if (CS%nk /= fld_sz(3)) call MOM_error(FATAL,'Increment levels /= ODA levels') @@ -578,9 +579,9 @@ subroutine get_bias_correction_tracer(Time, US, CS) call cpu_clock_begin(id_clock_bias_adjustment) - call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id, Time, CS%G, T_bias, & + call horiz_interp_and_extrap_tracer(CS%INC_CS%T, Time, CS%G, T_bias, & valid_flag, z_in, z_edges_in, missing_value, scale=US%degC_to_C*US%s_to_T, spongeOngrid=.true.) - call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id, Time, CS%G, S_bias, & + call horiz_interp_and_extrap_tracer(CS%INC_CS%S, Time, CS%G, S_bias, & valid_flag, z_in, z_edges_in, missing_value, scale=US%ppt_to_S*US%s_to_T, spongeOngrid=.true.) ! This should be replaced to use mask_z instead of the following lines diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 2a5cef5974..6a439dfd22 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -21,6 +21,7 @@ module MOM_MEKE use MOM_interface_heights, only : find_eta use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : vardesc, var_desc, slasher use MOM_isopycnal_slopes, only : calc_isoneutral_slopes use MOM_restart, only : MOM_restart_CS, register_restart_field, query_initialized @@ -129,7 +130,7 @@ module MOM_MEKE integer :: id_Lrhines = -1, id_Leady = -1 integer :: id_MEKE_equilibrium = -1 !>@} - integer :: id_eke = -1 !< Handle for reading in EKE from a file + type(external_field) :: eke_handle !< Handle for reading in EKE from a file ! Infrastructure integer :: id_clock_pass !< Clock for group pass calls type(group_pass_type) :: pass_MEKE !< Group halo pass handle for MEKE%MEKE and maybe MEKE%Kh_diff @@ -627,7 +628,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif case(EKE_FILE) - call time_interp_external(CS%id_eke, Time, data_eke, scale=US%m_s_to_L_T**2) + call time_interp_external(CS%eke_handle, Time, data_eke, scale=US%m_s_to_L_T**2) do j=js,je ; do i=is,ie MEKE%MEKE(i,j) = data_eke(i,j) * G%mask2dT(i,j) enddo; enddo @@ -1153,7 +1154,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, dbcomms_CS, CS, MEKE, inputdir = slasher(inputdir) eke_filename = trim(inputdir) // trim(eke_filename) - CS%id_eke = init_external_field(eke_filename, eke_varname, domain=G%Domain%mpp_domain) + CS%eke_handle = init_external_field(eke_filename, eke_varname, domain=G%Domain%mpp_domain) case("prog") CS%eke_src = EKE_PROG ! Read all relevant parameters and write them to the model log. diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 584ccccc93..2a30f68b42 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -22,6 +22,7 @@ module MOM_ALE_sponge use MOM_grid, only : ocean_grid_type use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer use MOM_interpolate, only : init_external_field, get_external_field_info, time_interp_external_init +use MOM_interpolate, only : external_field use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping use MOM_spatial_means, only : global_i_mean use MOM_time_manager, only : time_type @@ -66,7 +67,7 @@ module MOM_ALE_sponge !> A structure for creating arrays of pointers to 3D arrays with extra gridding information type :: p3d - integer :: id !< id for FMS external time interpolator + !integer :: id !< id for FMS external time interpolator integer :: nz_data !< The number of vertical levels in the input field. integer :: num_tlevs !< The number of time records contained in the file real, dimension(:,:,:), pointer :: p => NULL() !< pointer to the data [various] @@ -75,7 +76,7 @@ module MOM_ALE_sponge !> A structure for creating arrays of pointers to 2D arrays with extra gridding information type :: p2d - integer :: id !< id for FMS external time interpolator + type(external_field) :: field !< Time interpolator field handle integer :: nz_data !< The number of vertical levels in the input field integer :: num_tlevs !< The number of time records contained in the file real :: scale = 1.0 !< A multiplicative factor by which to rescale input data [various] @@ -771,7 +772,6 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, !! if not given, use 'none' real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any !! contributions due to dimensional rescaling [various ~> 1]. - !! The default is 1. ! Local variables integer :: isd, ied, jsd, jed @@ -798,15 +798,15 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, ! get a unique time interp id for this field. If sponge data is on-grid, then setup ! to only read on the computational domain if (CS%spongeDataOngrid) then - CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname, MOM_domain=G%Domain) + CS%Ref_val(CS%fldno)%field = init_external_field(filename, fieldname, MOM_domain=G%Domain) else - CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname) + CS%Ref_val(CS%fldno)%field = init_external_field(filename, fieldname) endif CS%Ref_val(CS%fldno)%name = sp_name CS%Ref_val(CS%fldno)%long_name = long_name CS%Ref_val(CS%fldno)%unit = unit fld_sz(1:4) = -1 - call get_external_field_info(CS%Ref_val(CS%fldno)%id, size=fld_sz) + call get_external_field_info(CS%Ref_val(CS%fldno)%field, size=fld_sz) nz_data = fld_sz(3) CS%Ref_val(CS%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid CS%Ref_val(CS%fldno)%num_tlevs = fld_sz(4) @@ -899,23 +899,23 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename ! containing time-interpolated values from an external file corresponding ! to the current model date. if (CS%spongeDataOngrid) then - CS%Ref_val_u%id = init_external_field(filename_u, fieldname_u, domain=G%Domain%mpp_domain) + CS%Ref_val_u%field = init_external_field(filename_u, fieldname_u, domain=G%Domain%mpp_domain) else - CS%Ref_val_u%id = init_external_field(filename_u, fieldname_u) + CS%Ref_val_u%field = init_external_field(filename_u, fieldname_u) endif fld_sz(1:4) = -1 - call get_external_field_info(CS%Ref_val_u%id, size=fld_sz) + call get_external_field_info(CS%Ref_val_u%field, size=fld_sz) CS%Ref_val_u%nz_data = fld_sz(3) CS%Ref_val_u%num_tlevs = fld_sz(4) CS%Ref_val_u%scale = US%m_s_to_L_T ; if (present(scale)) CS%Ref_val_u%scale = scale if (CS%spongeDataOngrid) then - CS%Ref_val_v%id = init_external_field(filename_v, fieldname_v, domain=G%Domain%mpp_domain) + CS%Ref_val_v%field = init_external_field(filename_v, fieldname_v, domain=G%Domain%mpp_domain) else - CS%Ref_val_v%id = init_external_field(filename_v, fieldname_v) + CS%Ref_val_v%field = init_external_field(filename_v, fieldname_v) endif fld_sz(1:4) = -1 - call get_external_field_info(CS%Ref_val_v%id, size=fld_sz) + call get_external_field_info(CS%Ref_val_v%field, size=fld_sz) CS%Ref_val_v%nz_data = fld_sz(3) CS%Ref_val_v%num_tlevs = fld_sz(4) CS%Ref_val_v%scale = US%m_s_to_L_T ; if (present(scale)) CS%Ref_val_v%scale = scale @@ -989,7 +989,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%time_varying_sponges) then do m=1,CS%fldno nz_data = CS%Ref_val(m)%nz_data - call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, G, sp_val, & + call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%field, Time, G, sp_val, & mask_z, z_in, z_edges_in, missing_value, & scale=CS%Ref_val(m)%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & answer_date=CS%hor_regrid_answer_date) @@ -1073,7 +1073,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%time_varying_sponges) then nz_data = CS%Ref_val_u%nz_data ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, G, sp_val, & + call horiz_interp_and_extrap_tracer(CS%Ref_val_u%field, Time, G, sp_val, & mask_z, z_in, z_edges_in, missing_value, & scale=CS%Ref_val_u%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & answer_date=CS%hor_regrid_answer_date) @@ -1121,7 +1121,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) deallocate(sp_val, mask_u, mask_z, hsrc) nz_data = CS%Ref_val_v%nz_data ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, G, sp_val, & + call horiz_interp_and_extrap_tracer(CS%Ref_val_v%field, Time, G, sp_val, & mask_z, z_in, z_edges_in, missing_value, & scale=CS%Ref_val_v%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& answer_date=CS%hor_regrid_answer_date) @@ -1341,7 +1341,7 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, GV, turns, param_file) ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying() ! (time_interp_external_init, init_external_field, etc), so we manually ! do a portion of this function below. - sponge%Ref_val(n)%id = sponge_in%Ref_val(n)%id + sponge%Ref_val(n)%field = sponge_in%Ref_val(n)%field sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs nz_data = sponge_in%Ref_val(n)%nz_data diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 5f7acd982b..3096fe72cd 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -16,6 +16,7 @@ module MOM_diabatic_aux use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : slasher use MOM_opacity, only : set_opacity, opacity_CS, extract_optics_slice, extract_optics_fields use MOM_opacity, only : optics_type, optics_nbands, absorbRemainingSW, sumSWoverBands @@ -64,7 +65,7 @@ module MOM_diabatic_aux !! is added with a temperature of the local SST. logical :: var_pen_sw !< If true, use one of the CHL_A schemes to determine the !! e-folding depth of incoming shortwave radiation. - integer :: sbc_chl !< An integer handle used in time interpolation of + type(external_field) :: sbc_chl !< A handle used in time interpolation of !! chlorophyll read from a file. logical :: chl_from_file !< If true, chl_a is read from a file. diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 2a5e3f8854..ef8e712b7a 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -1,4 +1,4 @@ -!> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover + !> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover !! provided via cap (only NUOPC cap is implemented so far). module MOM_CFC_cap @@ -19,7 +19,8 @@ module MOM_CFC_cap use MOM_restart, only : query_initialized, set_initialized, MOM_restart_CS use MOM_spatial_means, only : global_mass_int_EFP use MOM_time_manager, only : time_type -use time_interp_external_mod, only : init_external_field, time_interp_external +use MOM_interpolate, only : time_interp_external +use MOM_interpolate, only : external_field use MOM_tracer_registry, only : register_tracer use MOM_tracer_types, only : tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut @@ -428,7 +429,8 @@ end subroutine CFC_cap_surface_state !> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM !! concentration, and calculating the solubility, Schmidt number, and gas exchange. -subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id_cfc12_atm) +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, & + cfc11_atm_handle, cfc12_atm_handle) type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type type(surface), intent(in ) :: sfc_state !< A structure containing fields @@ -439,8 +441,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id real, intent(in ) :: Rho0 !< The mean ocean density [R ~> kg m-3] type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the !! CFC's concentration in the atmosphere. - integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. - integer, optional, intent(inout):: id_cfc12_atm !< id number for time_interp_external. + type(external_field), optional, intent(inout) :: cfc11_atm_handle !< Handle for time-interpolated CFC11 + type(external_field), optional, intent(inout) :: cfc12_atm_handle !< Handle for time-interpolated CFC12 ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & @@ -463,8 +465,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ! CFC11 ATM concentration - if (present(id_cfc11_atm) .and. (id_cfc11_atm /= -1)) then - call time_interp_external(id_cfc11_atm, Time, cfc11_atm) + if (present(cfc11_atm_handle)) then + call time_interp_external(cfc11_atm_handle, Time, cfc11_atm) ! convert from ppt (pico mol/mol) to mol/mol cfc11_atm = cfc11_atm * 1.0e-12 else @@ -474,8 +476,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id endif ! CFC12 ATM concentration - if (present(id_cfc12_atm) .and. (id_cfc12_atm /= -1)) then - call time_interp_external(id_cfc12_atm, Time, cfc12_atm) + if (present(cfc12_atm_handle)) then + call time_interp_external(cfc12_atm_handle, Time, cfc12_atm) ! convert from ppt (pico mol/mol) to mol/mol cfc12_atm = cfc12_atm * 1.0e-12 else From 1c6dd7fef5a16585e467b59893d4ee9499150afb Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 6 Jun 2023 20:35:24 -0400 Subject: [PATCH 13/15] FMS2: Remove MPP-based axis data access With removal of axis-based operations in FMS2 I/O, this patch removes references to these calls and replaces them with MOM `axes_info` types. References to FMS1 read into an `axistype`, but the contents are transferred to an `axis_info`. FMS2 directly populates the `axis_info` content. The `get_external_field_info` calls are modified to return `axis_info` rather than `axistype`. The redundant `get_axis_data` function is also removed from `MOM_interp_infra`, since `get_axis_info` provides an equivalent operation. Generally speaking, this is not an improvement of the codebase. The FMS1 layer does a redundant copy of data from `axistype` to `axis_info`. The FMS2 layer is significantly worse, and re-opens the file to read the axis data for each field! But if the intention is to leverage the existing API, then I don't think we have any choice at the moment. Assuming this is a relatively infrequent operation, this should not cause any measureable issues, but it needs to be watched carefully. --- config_src/infra/FMS1/MOM_interp_infra.F90 | 44 +++++++++++++++------ config_src/infra/FMS2/MOM_interp_infra.F90 | 32 ++++++--------- src/framework/MOM_horizontal_regridding.F90 | 10 ++--- 3 files changed, 49 insertions(+), 37 deletions(-) diff --git a/config_src/infra/FMS1/MOM_interp_infra.F90 b/config_src/infra/FMS1/MOM_interp_infra.F90 index e14233a64b..70bc99827e 100644 --- a/config_src/infra/FMS1/MOM_interp_infra.F90 +++ b/config_src/infra/FMS1/MOM_interp_infra.F90 @@ -4,9 +4,11 @@ module MOM_interp_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domain_infra, only : MOM_domain_type, domain2d +use MOM_io, only : axis_info +use MOM_io, only : set_axis_info use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type -use mpp_io_mod, only : axistype, mpp_get_axis_data +use mpp_io_mod, only : axistype, mpp_get_axis_data, mpp_get_atts use time_interp_external_mod, only : time_interp_external use time_interp_external_mod, only : init_external_field, time_interp_external_init use time_interp_external_mod, only : get_external_field_size @@ -157,13 +159,33 @@ end function get_extern_field_size !> get axes of an external field from field index -function get_extern_field_axes(index) +function get_extern_field_axes(index) result(axes) - integer, intent(in) :: index !< field index - type(axistype), dimension(4) :: get_extern_field_axes !< field axes + integer, intent(in) :: index !< FMS interpolation field index + type(axis_info) :: axes(4) !< MOM IO field axes handle - get_extern_field_axes = get_external_field_axes(index) + type(axistype), dimension(4) :: fms_axes(4) + ! FMS axis handles + character(len=32) :: name + ! Axis name + real, allocatable :: points(:) + ! Axis line points + integer :: length + ! Axis line point length + integer :: i + ! Loop index + fms_axes = get_external_field_axes(index) + + do i = 1, 4 + call mpp_get_atts(fms_axes(i), name=name, len=length) + + allocate(points(length)) + call mpp_get_axis_data(fms_axes(i), points) + call set_axis_info(axes(i), name=name, ax_data=points) + + deallocate(points) + enddo end function get_extern_field_axes @@ -180,12 +202,12 @@ end function get_extern_field_missing !> Get information about the external fields. subroutine get_external_field_info(field, size, axes, missing) - type(external_field), intent(in) :: field !< Handle for time interpolated external - !! field returned from a previous - !! call to init_external_field() - integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data - type(axistype), dimension(4), optional, intent(inout) :: axes !< Axis types for the input data - real, optional, intent(inout) :: missing !< Missing value for the input data + type(external_field), intent(in) :: field !< Handle for time interpolated external + !! field returned from a previous + !! call to init_external_field() + integer, optional, intent(inout) :: size(4) !< Dimension sizes for the input data + type(axis_info), optional, intent(inout) :: axes(4) !< Axis types for the input data + real, optional, intent(inout) :: missing !< Missing value for the input data if (present(size)) then size(1:4) = get_extern_field_size(field%id) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index 7964b3537f..db471d0fc1 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -4,9 +4,10 @@ module MOM_interp_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domain_infra, only : MOM_domain_type, domain2d +use MOM_io, only : axis_info +use MOM_io, only : get_var_axes_info use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type -use mpp_io_mod, only : axistype, mpp_get_axis_data use time_interp_external_mod, only : time_interp_external use time_interp_external_mod, only : init_external_field, time_interp_external_init use time_interp_external_mod, only : get_external_field_size @@ -16,7 +17,7 @@ module MOM_interp_infra public :: horiz_interp_type, horizontal_interp_init public :: time_interp_extern, init_extern_field, time_interp_extern_init -public :: get_external_field_info, axistype, get_axis_data +public :: get_external_field_info public :: run_horiz_interp, build_horiz_interp_weights public :: external_field @@ -135,15 +136,6 @@ subroutine build_horiz_interp_weights_2d_to_2d(Interp, lon_in, lat_in, lon_out, end subroutine build_horiz_interp_weights_2d_to_2d -!> Extracts and returns the axis data stored in an axistype. -subroutine get_axis_data( axis, dat ) - type(axistype), intent(in) :: axis !< An axis type - real, dimension(:), intent(out) :: dat !< The data in the axis variable - - call mpp_get_axis_data( axis, dat ) -end subroutine get_axis_data - - !> get size of an external field from field index function get_extern_field_size(index) @@ -156,13 +148,11 @@ end function get_extern_field_size !> get axes of an external field from field index -function get_extern_field_axes(index) - - integer, intent(in) :: index !< field index - type(axistype), dimension(4) :: get_extern_field_axes !< field axes - - get_extern_field_axes = get_external_field_axes(index) +function get_extern_field_axes(field) result(axes) + type(external_field), intent(in) :: field !< Field handle + type(axis_info), dimension(4) :: axes !< Field axes + call get_var_axes_info(field%filename, field%label, axes) end function get_extern_field_axes @@ -182,16 +172,16 @@ subroutine get_external_field_info(field, size, axes, missing) type(external_field), intent(in) :: field !< Handle for time interpolated external !! field returned from a previous !! call to init_external_field() - integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data - type(axistype), dimension(4), optional, intent(inout) :: axes !< Axis types for the input data - real, optional, intent(inout) :: missing !< Missing value for the input data + integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data + type(axis_info), dimension(4), optional, intent(inout) :: axes !< Axis types for the input data + real, optional, intent(inout) :: missing !< Missing value for the input data if (present(size)) then size(1:4) = get_extern_field_size(field%id) endif if (present(axes)) then - axes(1:4) = get_extern_field_axes(field%id) + axes(1:4) = get_extern_field_axes(field) endif if (present(missing)) then diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index bedf710582..883653d715 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -16,7 +16,7 @@ module MOM_horizontal_regridding use MOM_interpolate, only : time_interp_external use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights use MOM_interp_infra, only : horiz_interp_type, horizontal_interp_init -use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data +use MOM_interp_infra, only : get_external_field_info use MOM_interp_infra, only : external_field use MOM_time_manager, only : time_type use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data @@ -668,7 +668,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim] logical :: add_np type(horiz_interp_type) :: Interp - type(axistype), dimension(4) :: axes_data + type(axis_info), dimension(4) :: axes_data integer :: is, ie, js, je ! compute domain indices integer :: isg, ieg, jsg, jeg ! global extent integer :: isd, ied, jsd, jed ! data domain indices @@ -728,8 +728,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & if (PRESENT(spongeOngrid)) is_ongrid = spongeOngrid if (.not. is_ongrid) then allocate(lon_in(id), lat_in(jd)) - call get_axis_data(axes_data(1), lon_in) - call get_axis_data(axes_data(2), lat_in) + call get_axis_info(axes_data(1), ax_data=lon_in) + call get_axis_info(axes_data(2), ax_data=lat_in) endif allocate(z_in(kd), z_edges_in(kd+1)) @@ -737,7 +737,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & allocate(tr_z(isd:ied,jsd:jed,kd), source=0.0) allocate(mask_z(isd:ied,jsd:jed,kd), source=0.0) - call get_axis_data(axes_data(3), z_in) + call get_axis_info(axes_data(3), ax_data=z_in) if (present(m_to_Z)) then ; do k=1,kd ; z_in(k) = m_to_Z * z_in(k) ; enddo ; endif From e3c82d24668c7329e4699048dbb77a87de5fab52 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 7 Jun 2023 00:24:30 -0400 Subject: [PATCH 14/15] FMS2: Update time_interp_external functions This patch shifts all remaining time_interp_external functions from time_interp_external to equivalent ones in time_interp_external2. Internally, time-interpolated fields are initialized with `ongrid` set to `.true.`, and such fields are assumed to be on-grid. This seems to hold for all existing instances of `time_interp_external`, but needs to be monitored in the future somehow. --- config_src/infra/FMS2/MOM_interp_infra.F90 | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index db471d0fc1..09a9eb0421 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -8,10 +8,10 @@ module MOM_interp_infra use MOM_io, only : get_var_axes_info use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type -use time_interp_external_mod, only : time_interp_external -use time_interp_external_mod, only : init_external_field, time_interp_external_init -use time_interp_external_mod, only : get_external_field_size -use time_interp_external_mod, only : get_external_field_axes, get_external_field_missing +use time_interp_external2_mod, only : time_interp_external +use time_interp_external2_mod, only : init_external_field, time_interp_external_init +use time_interp_external2_mod, only : get_external_field_size +use time_interp_external2_mod, only : get_external_field_missing implicit none ; private @@ -267,12 +267,14 @@ function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & if (present(MOM_Domain)) then field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & - correct_leap_year_inconsistency=correct_leap_year_inconsistency) + verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency, & + ongrid=.true.) else field%id = init_external_field(file, fieldname, domain=domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & - correct_leap_year_inconsistency=correct_leap_year_inconsistency) + verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency, & + ongrid=.true.) endif end function init_extern_field From dd1ee3422252e1bdd8613f3fb58d7807bc59d3c4 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 16 Jun 2023 10:40:57 -0400 Subject: [PATCH 15/15] FMS2: Case-insensitive init_external_field The FMS1 implementation of init_external_field is case-insensitive, but the FMS2 implementation is case-sensitive, which can cause errors in older established input files. This patch sweeps through the fields of the input files and checks for a case-insensitive match (using lowercase()). This requires an additional open/close of the file. --- config_src/infra/FMS2/MOM_interp_infra.F90 | 60 ++++++++++++++++++++-- 1 file changed, 56 insertions(+), 4 deletions(-) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index 09a9eb0421..0b45b752ae 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -7,7 +7,11 @@ module MOM_interp_infra use MOM_io, only : axis_info use MOM_io, only : get_var_axes_info use MOM_time_manager, only : time_type -use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use MOM_error_handler, only : MOM_error, FATAL +use MOM_string_functions, only : lowercase +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use netcdf_io_mod, only : FmsNetcdfFile_t, netcdf_file_open, netcdf_file_close +use netcdf_io_mod, only : get_num_variables, get_variable_names use time_interp_external2_mod, only : time_interp_external use time_interp_external2_mod, only : init_external_field, time_interp_external_init use time_interp_external2_mod, only : get_external_field_size @@ -262,16 +266,64 @@ function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & !! a model date of Feb 29. onto a common year on Feb. 28. type(external_field) :: field !< Handle to external field + type(FmsNetcdfFile_t) :: extern_file + ! Local instance of netCDF file used to locate case-insensitive field name + integer :: num_fields + ! Number of fields in external file + character(len=256), allocatable :: extern_fieldnames(:) + ! List of field names in file + ! NOTE: length should NF90_MAX_NAME, but I don't know how to read it + character(len=:), allocatable :: label + ! Case-insensitive match to fieldname in file + logical :: rc + ! Return status + integer :: i + ! Loop index + field%filename = file - field%label = fieldname + + ! FMS2's init_external_field is case sensitive, so we must replicate the + ! case-insensitivity of FMS1. This requires opening the file twice. + + rc = netcdf_file_open(extern_file, file, 'read') + if (.not. rc) then + call MOM_error(FATAL, 'init_extern_file: file ' // trim(file) & + // ' could not be opened.') + endif + + ! TODO: broadcast = .false.? + num_fields = get_num_variables(extern_file) + + allocate(extern_fieldnames(num_fields)) + call get_variable_names(extern_file, extern_fieldnames) + + do i = 1, num_fields + if (lowercase(extern_fieldnames(i)) == lowercase(fieldname)) then + field%label = extern_fieldnames(i) + exit + endif + enddo + + call netcdf_file_close(extern_file) + + if (.not. allocated(field%label)) then + call MOM_error(FATAL, 'init_extern_field: field ' // trim(fieldname) & + // ' not found in ' // trim(file) // '.') + endif + + ! Pass to FMS2 implementation of init_external_field + + ! NOTE: external fields are currently assumed to be on-grid, which holds + ! across the current codebase. In the future, we may need to either enforce + ! this or somehow relax this requirement. if (present(MOM_Domain)) then - field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & + field%id = init_external_field(file, field%label, domain=MOM_domain%mpp_domain, & verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency, & ongrid=.true.) else - field%id = init_external_field(file, fieldname, domain=domain, & + field%id = init_external_field(file, field%label, domain=domain, & verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency, & ongrid=.true.)