MOM_self_attr_load.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
6
7use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_module
8use mom_domains, only : pass_var
9use mom_error_handler, only : mom_error, fatal, warning
10use mom_file_parser, only : get_param, log_param, log_version, param_file_type
11use mom_grid, only : ocean_grid_type
13use mom_io, only : mom_infra_file, mom_field, vardesc, slasher
14use mom_io, only : create_mom_file, mom_read_data, mom_write_field, var_desc
15use mom_load_love_numbers, only : love_data
24
25implicit none ; private
26
28
29#include <MOM_memory.h>
30
31!> The control structure for the MOM_self_attr_load module
32type, public :: sal_cs ; private
33 logical :: use_sal_scalar = .false.
34 !< If true, use the scalar approximation to calculate SAL.
35 logical :: use_sal_sht = .false.
36 !< If true, use online spherical harmonics to calculate SAL
37 logical :: use_tidal_sal_prev = .false.
38 !< If true, read the tidal SAL from the previous iteration of the tides to
39 !! facilitate convergence.
40 logical :: use_bpa = .false.
41 !< If true, use bottom pressure anomaly instead of SSH to calculate SAL.
42 real :: eta_prop
43 !< The partial derivative of eta_sal with the local value of eta [nondim].
44 real :: linear_scaling
45 !< Dimensional coefficients for scalar SAL [nondim] or [Z T2 L-2 R-1 ~> m Pa-1]
46 type(sht_cs), allocatable :: sht
47 !< Spherical harmonic transforms (SHT) control structure
48 integer :: sal_sht_nd
49 !< Maximum degree for spherical harmonic transforms [nondim]
50 real, allocatable :: pbot_ref(:,:)
51 !< Reference bottom pressure [R L2 T-2 ~> Pa]
52 real, allocatable :: love_scaling(:)
53 !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers
54 !! [nondim] or [Z T2 L-2 R-1 ~> m Pa-1], depending on the value of use_ppa.
55 real, allocatable :: snm_re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m]
56 snm_im(:) !< Imaginary SHT coefficient for SHT SAL [Z ~> m]
57end type sal_cs
58
59integer :: id_clock_sal !< CPU clock for self-attraction and loading
60
61contains
62
63!> This subroutine calculates seawater self-attraction and loading based on either sea surface height (SSH) or bottom
64!! pressure anomaly. Note that the SAL calculation applies to all motions across the spectrum. Tidal-specific methods
65!! that assume periodicity, i.e. iterative and read-in SAL, are stored in MOM_tidal_forcing module.
66!! The input field can be either SSH [Z ~> m] or total bottom pressure [R L2 T-2 ~> Pa]. If total bottom pressure is
67!! used, bottom pressure anomaly is first calculated by subtracting a reference bottom pressure from an input file.
68!! The output field is expressed as geopotential height anomaly, and therefore has the unit of [Z ~> m].
69subroutine calc_sal(eta, eta_sal, G, CS, tmp_scale)
70 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
71 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta !< The sea surface height anomaly from
72 !! a time-mean geoid or total bottom pressure [Z ~> m] or [R L2 T-2 ~> Pa].
73 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The geopotential height anomaly from
74 !! self-attraction and loading [Z ~> m].
75 type(sal_cs), intent(inout) :: cs !< The control structure returned by a previous call to SAL_init.
76 real, optional, intent(in) :: tmp_scale !< A rescaling factor to temporarily convert eta
77 !! to MKS units in reproducing sumes [m Z-1 ~> 1]
78
79 ! Local variables
80 real, dimension(SZI_(G),SZJ_(G)) :: bpa ! SSH or bottom pressure anomaly [Z ~> m] or [R L2 T-2 ~> Pa]
81 integer :: n, m, l
82 integer :: isq, ieq, jsq, jeq
83 integer :: i, j
84
85 call cpu_clock_begin(id_clock_sal)
86
87 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
88
89 if (cs%use_bpa) then ; do j=jsq,jeq+1 ; do i=isq,ieq+1
90 bpa(i,j) = eta(i,j) - cs%pbot_ref(i,j)
91 enddo ; enddo ; else ; do j=jsq,jeq+1 ; do i=isq,ieq+1
92 bpa(i,j) = eta(i,j)
93 enddo ; enddo ; endif
94
95 ! use the scalar approximation and/or iterative tidal SAL
96 if (cs%use_sal_scalar .or. cs%use_tidal_sal_prev) then
97 do j=jsq,jeq+1 ; do i=isq,ieq+1
98 eta_sal(i,j) = cs%linear_scaling * bpa(i,j)
99 enddo ; enddo
100
101 ! use the spherical harmonics method
102 elseif (cs%use_sal_sht) then
103 call spherical_harmonics_forward(g, cs%sht, bpa, cs%Snm_Re, cs%Snm_Im, cs%sal_sht_Nd, tmp_scale=tmp_scale)
104
105 ! Multiply scaling factors to each mode
106 do m = 0,cs%sal_sht_Nd
107 l = order2index(m, cs%sal_sht_Nd)
108 do n = m,cs%sal_sht_Nd
109 cs%Snm_Re(l+n-m) = cs%Snm_Re(l+n-m) * cs%Love_scaling(l+n-m)
110 cs%Snm_Im(l+n-m) = cs%Snm_Im(l+n-m) * cs%Love_scaling(l+n-m)
111 enddo
112 enddo
113
114 call spherical_harmonics_inverse(g, cs%sht, cs%Snm_Re, cs%Snm_Im, eta_sal, cs%sal_sht_Nd)
115 ! Halo was not calculated in spherical harmonic transforms.
116 call pass_var(eta_sal, g%domain)
117
118 else
119 do j=jsq,jeq+1 ; do i=isq,ieq+1
120 eta_sal(i,j) = 0.0
121 enddo ; enddo
122 endif
123
124 call cpu_clock_end(id_clock_sal)
125end subroutine calc_sal
126
127!> This subroutine returns eta_prop member of SAL_CS type, which is the non-dimensional partial
128!! derivative of the local geopotential height with the input sea surface height due to the scalar
129!! approximation of self-attraction and loading.
130subroutine scalar_sal_sensitivity(CS, deta_sal_deta)
131 type(sal_cs), intent(in) :: cs !< The control structure returned by a previous call to SAL_init.
132 real, intent(out) :: deta_sal_deta !< The partial derivative of eta_sal with
133 !! the local value of eta [nondim].
134 deta_sal_deta = cs%eta_prop
135end subroutine scalar_sal_sensitivity
136
137!> This subroutine calculates coefficients of the spherical harmonic modes for self-attraction and loading.
138!! The algorithm is based on the SAL implementation in MPAS-ocean, which was modified by Kristin Barton from
139!! routine written by K. Quinn (March 2010) and modified by M. Schindelegger (May 2017).
140subroutine calc_love_scaling(rhoW, rhoE, grav, CS)
141 real, intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3]
142 real, intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3]
143 real, intent(in) :: grav !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
144 type(sal_cs), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init.
145
146 ! Local variables
147 real :: coef_rhoE ! A scaling coefficient of solid Earth density. coef_rhoE = rhoW / rhoE with USE_BPA=False
148 ! and coef_rhoE = 1.0 / (rhoE * grav) with USE_BPA=True. [nondim] or [Z T2 L-2 R-1 ~> m Pa-1]
149 real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames [nondim]
150 real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers [nondim]
151 integer :: n_tot ! Size of the stored Love numbers [nondim]
152 integer :: nlm ! Maximum spherical harmonics degree [nondim]
153 integer :: n, m, l
154
155 n_tot = size(love_data, dim=2)
156 nlm = cs%sal_sht_Nd
157
158 if (nlm+1 > n_tot) call mom_error(fatal, "MOM_tidal_forcing " // &
159 "calc_love_scaling: maximum spherical harmonics degree is larger than " // &
160 "the size of the stored Love numbers in MOM_load_love_number.")
161
162 allocate(hdat(nlm+1), ldat(nlm+1), kdat(nlm+1))
163 hdat(:) = love_data(2,1:nlm+1) ; ldat(:) = love_data(3,1:nlm+1) ; kdat(:) = love_data(4,1:nlm+1)
164
165 ! Convert reference frames from CM to CF
166 if (nlm > 0) then
167 h1 = hdat(2) ; l1 = ldat(2) ; k1 = kdat(2)
168 hdat(2) = ( 2.0 / 3.0) * (h1 - l1)
169 ldat(2) = (-1.0 / 3.0) * (h1 - l1)
170 kdat(2) = (-1.0 / 3.0) * h1 - (2.0 / 3.0) * l1 - 1.0
171 endif
172
173 if (cs%use_bpa) then
174 coef_rhoe = 1.0 / (rhoe * grav) ! [Z T2 L-2 R-1 ~> m Pa-1]
175 else
176 coef_rhoe = rhow / rhoe ! [nondim]
177 endif
178
179 do m=0,nlm ; do n=m,nlm
180 l = order2index(m, nlm)
181 ! Love_scaling has the same as coef_rhoE.
182 cs%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * coef_rhoe * (1.0 + kdat(n+1) - hdat(n+1))
183 enddo ; enddo
184end subroutine calc_love_scaling
185
186!> This subroutine initializes the self-attraction and loading control structure.
187subroutine sal_init(h, tv, G, GV, US, param_file, CS, restart_CS)
188 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
189 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
190 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
191 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
192 type(sal_cs), intent(inout) :: cs !< Self-attraction and loading control structure
193 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables.
194 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
195 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
196 type(mom_restart_cs), optional, intent(in) :: restart_cs !< MOM restart control structure
197 ! Local variables
198# include "version_variable.h"
199 character(len=40) :: mdl = "MOM_self_attr_load" ! This module's name.
200 integer :: lmax ! Total modes of the real spherical harmonics [nondim]
201 real :: rhoe ! The average density of Earth [R ~> kg m-3].
202 character(len=20) :: bpa_config ! String for reference bottom pressure config option
203 real :: tmp(g%isd:g%ied, g%jsd:g%jed) ! Temporary field storing mass returned by find_col_mass
204 ! [R Z ~> kg m-2]
205 logical :: restart_sim ! If true, this is a restart run
206 character(len=200) :: filename, ref_pbot_file, inputdir ! Strings for file/path
207 character(len=200) :: ref_pbot_varname ! Variable name in file
208 type(mom_infra_file) :: io_handle ! used to write ref_pbot file
209 type(vardesc) :: vars(1) ! used to write ref_pbot file
210 type(mom_field) :: fields(1) ! used to write ref_pbot file
211 logical :: calculate_sal, tides, use_tidal_sal_file
212 integer :: default_answer_date, tides_answer_date ! Recover old answers with tides
213 real :: sal_scalar_value ! Scaling SAL factors [nondim]
214 integer :: isd, ied, jsd, jed
215
216 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
217
218 ! Read all relevant parameters and write them to the model log.
219 call log_version(param_file, mdl, version, "")
220
221 call get_param(param_file, '', "TIDES", tides, default=.false., do_not_log=.true.)
222 call get_param(param_file, '', "CALCULATE_SAL", calculate_sal, default=tides, do_not_log=.true.)
223 if (.not. calculate_sal) return
224
225 call get_param(param_file, mdl, "SAL_USE_BPA", cs%use_bpa, &
226 "If true, use bottom pressure anomaly to calculate self-attraction and "// &
227 "loading (SAL). Otherwise sea surface height anomaly is used, which is "// &
228 "only accurate for uniform density fluid.", default=.false.)
229 if (cs%use_bpa) then
230 allocate(cs%pbot_ref(isd:ied, jsd:jed), source=0.0)
231 call get_param(param_file, mdl, "SAL_REF_PBOT_CONFIG", bpa_config, default="file", &
232 do_not_log=.true.)
233 restart_sim = .false. ; if (present(restart_cs)) restart_sim = (.not. is_new_run(restart_cs))
234 if (restart_sim .and. (trim(lowercase(bpa_config))/='file')) then
235 call mom_error(warning, "SAL_init: 'file' is not used by SAL_PBOT_REF_CONFIG for a restart "//&
236 "run, SAL_PBOT_REF_CONFIG is reset to 'file'.")
237 bpa_config = 'file'
238 endif
239 call get_param(param_file, mdl, "SAL_REF_PBOT_CONFIG", bpa_config, &
240 "A string that determines how the reference bottom pressure for SAL "//&
241 "is specified:\n"//&
242 "\t init - calculated by thickness, temperature and salinity from \n"//&
243 "\t initialization and assuming surface pressure is zero.\n"//&
244 "\t This option can only be used by new simulations.\n"//&
245 "\t file - read from the file specified by REF_PBOT_FILE.", &
246 default="file", do_not_read=.true.)
247 call get_param(param_file, '', "INPUTDIR", inputdir, default=".", do_not_log=.true.)
248 call get_param(param_file, mdl, "REF_PBOT_FILE", ref_pbot_file, &
249 "Reference bottom pressure file used by self-attraction and loading (SAL).", &
250 default="pbot.nc")
251 call get_param(param_file, mdl, "REF_PBOT_VARNAME", ref_pbot_varname, &
252 "The name of the variable in REF_PBOT_FILE with reference bottom "//&
253 "pressure. The variable should have the unit of Pa.", default="pbot")
254 filename = trim(slasher(inputdir))//trim(ref_pbot_file)
255 call log_param(param_file, mdl, "INPUTDIR/REF_PBOT_FILE", filename)
256 select case (trim(lowercase(bpa_config)))
257 case ("file")
258 call mom_read_data(filename, trim(ref_pbot_varname), cs%pbot_ref, g%Domain,&
259 scale=us%Pa_to_RL2_T2)
260 case ("init")
261 call find_col_mass(h, tv, g, gv, us, tmp, cs%pbot_ref)
262 ! Write reference bottom pressure file
263 vars(1) = var_desc(trim(ref_pbot_varname), units="Pa", &
264 longname="Reference bottom pressure", &
265 hor_grid='h', z_grid='1', t_grid='1')
266 call create_mom_file(io_handle, trim(filename), vars, 1, fields, g=g)
267 call mom_write_field(io_handle, fields(1), g%Domain, cs%pbot_ref, unscale=us%RL2_T2_to_Pa)
268 call io_handle%close()
269 case default
270 call mom_error(fatal, "SAL_init: Unsupported SAL_PBOT_REF_CONFIG option "//trim(bpa_config))
271 end select
272 call pass_var(cs%pbot_ref, g%Domain)
273 endif
274
275 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
276 "This sets the default value for the various _ANSWER_DATE parameters.", &
277 default=99991231, do_not_log=.true.) ! used to check SAL_USE_BPA
278 call get_param(param_file, '', "TIDES_ANSWER_DATE", tides_answer_date, &
279 default=default_answer_date, do_not_log=.true.) ! used to check SAL_USE_BPA
280 if (tides_answer_date<=20250131 .and. cs%use_bpa) &
281 call mom_error(fatal, trim(mdl) // ", SAL_init: SAL_USE_BPA needs to be false to recover "//&
282 "tide answers before 20250131.")
283 call get_param(param_file, '', "TIDAL_SAL_FROM_FILE", use_tidal_sal_file, default=.false., &
284 do_not_log=.true.) ! used to set default of SAL_SCALAR_APPROX
285 call get_param(param_file, mdl, "SAL_SCALAR_APPROX", cs%use_sal_scalar, &
286 "If true, use the scalar approximation to calculate self-attraction and "//&
287 "loading.", default=tides .and. (.not. use_tidal_sal_file))
288 if (cs%use_sal_scalar .and. cs%use_bpa) &
289 call mom_error(warning, trim(mdl) // ", SAL_init: Using bottom pressure anomaly for scalar "//&
290 "approximation SAL is unsubstantiated.")
291 call get_param(param_file, mdl, "SAL_SCALAR_VALUE", sal_scalar_value, "The constant of "//&
292 "proportionality between self-attraction and loading (SAL) geopotential "//&
293 "anomaly and barotropic geopotential anomaly. This is only used if "//&
294 "SAL_SCALAR_APPROX is true or USE_PREVIOUS_TIDES is true.", default=0.0, &
295 units="m m-1", do_not_log=.not.(cs%use_sal_scalar .or. cs%use_tidal_sal_prev), &
296 old_name='TIDE_SAL_SCALAR_VALUE')
297 call get_param(param_file, '', "USE_PREVIOUS_TIDES", cs%use_tidal_sal_prev, &
298 default=.false., do_not_log=.true.)
299 call get_param(param_file, mdl, "SAL_HARMONICS", cs%use_sal_sht, &
300 "If true, use the online spherical harmonics method to calculate "//&
301 "self-attraction and loading.", default=.false.)
302 call get_param(param_file, mdl, "SAL_HARMONICS_DEGREE", cs%sal_sht_Nd, &
303 "The maximum degree of the spherical harmonics transformation used for "// &
304 "calculating the self-attraction and loading term.", &
305 default=0, do_not_log=(.not. cs%use_sal_sht))
306 call get_param(param_file, mdl, "RHO_SOLID_EARTH", rhoe, &
307 "The mean solid earth density. This is used for calculating the "// &
308 "self-attraction and loading term.", units="kg m-3", &
309 default=5517.0, scale=us%kg_m3_to_R, do_not_log=(.not. cs%use_sal_sht))
310
311 ! Set scaling coefficients for scalar approximation
312 if (cs%use_sal_scalar .or. cs%use_tidal_sal_prev) then
313 if (cs%use_sal_scalar .and. cs%use_tidal_sal_prev) then
314 cs%eta_prop = 2.0 * sal_scalar_value
315 else
316 cs%eta_prop = sal_scalar_value
317 endif
318 if (cs%use_bpa) then
319 cs%linear_scaling = cs%eta_prop / (gv%Rho0 * gv%g_Earth)
320 else
321 cs%linear_scaling = cs%eta_prop
322 endif
323 else
324 cs%eta_prop = 0.0 ; cs%linear_scaling = 0.0
325 endif
326
327 ! Set scaling coefficients for spherical harmonics
328 if (cs%use_sal_sht) then
329 lmax = calc_lmax(cs%sal_sht_Nd)
330 allocate(cs%Snm_Re(lmax), source=0.0)
331 allocate(cs%Snm_Im(lmax), source=0.0)
332
333 allocate(cs%Love_scaling(lmax), source=0.0)
334 call calc_love_scaling(gv%Rho0, rhoe, gv%g_Earth, cs)
335
336 allocate(cs%sht)
337 call spherical_harmonics_init(g, param_file, cs%sht)
338 endif
339
340 id_clock_sal = cpu_clock_id('(Ocean SAL)', grain=clock_module)
341
342end subroutine sal_init
343
344!> This subroutine deallocates memory associated with the SAL module.
345subroutine sal_end(CS)
346 type(sal_cs), intent(inout) :: cs !< The control structure returned by a previous call
347 !! to SAL_init; it is deallocated here.
348
349 if (allocated(cs%pbot_ref)) deallocate(cs%pbot_ref)
350
351 if (cs%use_sal_sht) then
352 if (allocated(cs%Love_scaling)) deallocate(cs%Love_scaling)
353 if (allocated(cs%Snm_Re)) deallocate(cs%Snm_Re)
354 if (allocated(cs%Snm_Im)) deallocate(cs%Snm_Im)
355 call spherical_harmonics_end(cs%sht)
356 deallocate(cs%sht)
357 endif
358end subroutine sal_end
359
360!> \namespace self_attr_load
361!!
362!! \section section_SAL Self attraction and loading
363!!
364!! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height or
365!! bottom pressure anomaly. SAL is primarily used for fast evolving processes like tides or storm surges, but the
366!! effect applies to all motions.
367!!
368!! If <code>SAL_SCALAR_APPROX</code> is true, a scalar approximation is applied (\cite Accad1978) and the SAL is simply
369!! a fraction (set by <code>SAL_SCALAR_VALUE</code>, usually around 10% for global tides) of local SSH. For tides, the
370!! scalar approximation can also be used to iterate the SAL to convergence [see <code>USE_PREVIOUS_TIDES</code> in
371!! MOM_tidal_forcing, \cite Arbic2004].
372!!
373!! If <code>SAL_HARMONICS</code> is true, a more accurate online spherical harmonic transforms are used to calculate
374!! SAL. Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is set
375!! by <code>SAL_HARMONICS_DEGREE</code>. The algorithm is based on SAL calculation in Model for Prediction Across
376!! Scales (MPAS)-Ocean developed by Los Alamos National Laboratory and University of Michigan
377!! [\cite Barton2022 and \cite Brus2023].
378!!
379!! References:
380!!
381!! Accad, Y. and Pekeris, C.L., 1978. Solution of the tidal equations for the M2 and S2 tides in the world oceans from a
382!! knowledge of the tidal potential alone. Philosophical Transactions of the Royal Society of London. Series A,
383!! Mathematical and Physical Sciences, 290(1368), pp.235-266.
384!! https://doi.org/10.1098/rsta.1978.0083
385!!
386!! Arbic, B.K., Garner, S.T., Hallberg, R.W. and Simmons, H.L., 2004. The accuracy of surface elevations in forward
387!! global barotropic and baroclinic tide models. Deep Sea Research Part II: Topical Studies in Oceanography, 51(25-26),
388!! pp.3069-3101.
389!! https://doi.org/10.1016/j.dsr2.2004.09.014
390!!
391!! Barton, K.N., Pal, N., Brus, S.R., Petersen, M.R., Arbic, B.K., Engwirda, D., Roberts, A.F., Westerink, J.J.,
392!! Wirasaet, D. and Schindelegger, M., 2022. Global Barotropic Tide Modeling Using Inline Self‐Attraction and Loading in
393!! MPAS‐Ocean. Journal of Advances in Modeling Earth Systems, 14(11), p.e2022MS003207.
394!! https://doi.org/10.1029/2022MS003207
395!!
396!! Brus, S.R., Barton, K.N., Pal, N., Roberts, A.F., Engwirda, D., Petersen, M.R., Arbic, B.K., Wirasaet, D.,
397!! Westerink, J.J. and Schindelegger, M., 2023. Scalable self attraction and loading calculations for unstructured ocean
398!! tide models. Ocean Modelling, p.102160.
399!! https://doi.org/10.1016/j.ocemod.2023.102160
400end module mom_self_attr_load