MOM_spherical_harmonics.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
5!> Laplace's spherical harmonic transforms (SHT)
7use mom_coms_infra, only : sum_across_pes
8use mom_coms, only : reproducing_sum
9use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, &
10 clock_module, clock_routine, clock_loop
11use mom_error_handler, only : mom_error, fatal
12use mom_file_parser, only : get_param, log_version, param_file_type
13use mom_grid, only : ocean_grid_type
14
15implicit none ; private
16
19
20#include <MOM_memory.h>
21
22!> Control structure for spherical harmonic transforms
23type, public :: sht_cs ; private
24 logical :: initialized = .false. !< True if this control structure has been initialized.
25 integer :: ndegree !< Maximum degree of the spherical harmonics [nondim].
26 integer :: lmax !< Number of associated Legendre polynomials of nonnegative m
27 !! [lmax=(ndegree+1)*(ndegree+2)/2] [nondim].
28 real, allocatable :: cos_clatt(:,:) !< Precomputed cosine of colatitude at the t-cells [nondim].
29 real, allocatable :: pmm(:,:,:) !< Precomputed associated Legendre polynomials (m=n) at the t-cells [nondim].
30 real, allocatable :: cos_lont(:,:,:), & !< Precomputed cosine factors at the t-cells [nondim].
31 sin_lont(:,:,:) !< Precomputed sine factors at the t-cells [nondim].
32 real, allocatable :: cos_lont_wtd(:,:,:), & !< Precomputed area-weighted cosine factors at the t-cells [nondim]
33 sin_lont_wtd(:,:,:) !< Precomputed area-weighted sine factors at the t-cells [nondim]
34 real, allocatable :: a_recur(:,:), & !< Precomputed recurrence coefficients a [nondim].
35 b_recur(:,:) !< Precomputed recurrence coefficients b [nondim].
36 logical :: reprod_sum !< True if use reproducible global sums
37end type sht_cs
38
39integer :: id_clock_sht=-1 !< CPU clock for SHT [MODULE]
40integer :: id_clock_sht_forward=-1 !< CPU clock for forward transforms [ROUTINE]
41integer :: id_clock_sht_inverse=-1 !< CPU clock for inverse transforms [ROUTINE]
42integer :: id_clock_sht_global_sum=-1 !< CPU clock for global summation in forward transforms [LOOP]
43
44contains
45
46!> Calculates forward spherical harmonics transforms
47subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale)
48 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
49 type(sht_cs), intent(inout) :: cs !< Control structure for SHT
50 real, dimension(SZI_(G),SZJ_(G)), &
51 intent(in) :: var !< Input 2-D variable in arbitrary mks units [a]
52 !! or in arbitrary rescaled units [A ~> a] if
53 !! tmp_scale is present
54 real, intent(out) :: snm_re(:) !< SHT coefficients for the real modes (cosine) in
55 !! the same arbitrary units as var [a] or [A ~> a]
56 real, intent(out) :: snm_im(:) !< SHT coefficients for the imaginary modes (sine) in
57 !! the same arbitrary units as var [a] or [A ~> a]
58 integer, optional, intent(in) :: nd !< Maximum degree of the spherical harmonics
59 !! overriding ndegree in the CS [nondim]
60 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor to convert
61 !! var to MKS units during the reproducing
62 !! sums [a A-1 ~> 1]
63 ! local variables
64 integer :: nmax ! Local copy of the maximum degree of the spherical harmonics
65 integer :: ltot ! Local copy of the number of spherical harmonics
66 real, dimension(SZI_(G),SZJ_(G)) :: &
67 pmn, & ! Current associated Legendre polynomials of degree n and order m [nondim]
68 pmnm1, & ! Associated Legendre polynomials of degree n-1 and order m [nondim]
69 pmnm2 ! Associated Legendre polynomials of degree n-2 and order m [nondim]
70 real, allocatable, dimension(:,:,:) :: &
71 snm_re_raw, & ! Array of un-summed real spherical harmonics transform coefficients for
72 ! reproducing sums in the same arbitrary units as var, [a] or [A ~> a]
73 snm_im_raw ! Array of un-summed imaginary spherical harmonics transform coefficients for
74 ! reproducing sums in the same arbitrary units as var, [a] or [A ~> a]
75 real :: sum_tot ! The total of all components output by the reproducing sum in the same
76 ! arbitrary units as var, [a] or [A ~> a]
77 integer :: i, j, k
78 integer :: is, ie, js, je, isd, ied, jsd, jed
79 integer :: m, n, l
80
81 if (.not.cs%initialized) call mom_error(fatal, "MOM_spherical_harmonics " // &
82 "spherical_harmonics_forward: Module must be initialized before it is used.")
83
84 if (id_clock_sht>0) call cpu_clock_begin(id_clock_sht)
85 if (id_clock_sht_forward>0) call cpu_clock_begin(id_clock_sht_forward)
86
87 nmax = cs%ndegree ; if (present(nd)) nmax = nd
88 ltot = calc_lmax(nmax)
89
90 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
91 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
92
93 do j=jsd,jed ; do i=isd,ied
94 pmn(i,j) = 0.0 ; pmnm1(i,j) = 0.0 ; pmnm2(i,j) = 0.0
95 enddo ; enddo
96
97 do l=1,ltot ; snm_re(l) = 0.0 ; snm_im(l) = 0.0 ; enddo
98
99 if (cs%reprod_sum) then
100 allocate(snm_re_raw(is:ie, js:je, ltot), source=0.0)
101 allocate(snm_im_raw(is:ie, js:je, ltot), source=0.0)
102 do m=0,nmax
103 l = order2index(m, nmax)
104
105 do j=js,je ; do i=is,ie
106 snm_re_raw(i,j,l) = var(i,j) * cs%Pmm(i,j,m+1) * cs%cos_lonT_wtd(i,j,m+1)
107 snm_im_raw(i,j,l) = var(i,j) * cs%Pmm(i,j,m+1) * cs%sin_lonT_wtd(i,j,m+1)
108 pmnm2(i,j) = 0.0
109 pmnm1(i,j) = cs%Pmm(i,j,m+1)
110 enddo ; enddo
111
112 do n = m+1, nmax ; do j=js,je ; do i=is,ie
113 pmn(i,j) = &
114 cs%a_recur(n+1,m+1) * cs%cos_clatT(i,j) * pmnm1(i,j) - cs%b_recur(n+1,m+1) * pmnm2(i,j)
115 snm_re_raw(i,j,l+n-m) = var(i,j) * pmn(i,j) * cs%cos_lonT_wtd(i,j,m+1)
116 snm_im_raw(i,j,l+n-m) = var(i,j) * pmn(i,j) * cs%sin_lonT_wtd(i,j,m+1)
117 pmnm2(i,j) = pmnm1(i,j)
118 pmnm1(i,j) = pmn(i,j)
119 enddo ; enddo ; enddo
120 enddo
121 else
122 do m=0,nmax
123 l = order2index(m, nmax)
124
125 do j=js,je ; do i=is,ie
126 snm_re(l) = snm_re(l) + var(i,j) * cs%Pmm(i,j,m+1) * cs%cos_lonT_wtd(i,j,m+1)
127 snm_im(l) = snm_im(l) + var(i,j) * cs%Pmm(i,j,m+1) * cs%sin_lonT_wtd(i,j,m+1)
128 pmnm2(i,j) = 0.0
129 pmnm1(i,j) = cs%Pmm(i,j,m+1)
130 enddo ; enddo
131
132 do n=m+1, nmax ; do j=js,je ; do i=is,ie
133 pmn(i,j) = &
134 cs%a_recur(n+1,m+1) * cs%cos_clatT(i,j) * pmnm1(i,j) - cs%b_recur(n+1,m+1) * pmnm2(i,j)
135 snm_re(l+n-m) = snm_re(l+n-m) + var(i,j) * pmn(i,j) * cs%cos_lonT_wtd(i,j,m+1)
136 snm_im(l+n-m) = snm_im(l+n-m) + var(i,j) * pmn(i,j) * cs%sin_lonT_wtd(i,j,m+1)
137 pmnm2(i,j) = pmnm1(i,j)
138 pmnm1(i,j) = pmn(i,j)
139 enddo ; enddo ; enddo
140 enddo
141 endif
142
143 if (id_clock_sht_global_sum>0) call cpu_clock_begin(id_clock_sht_global_sum)
144
145 if (cs%reprod_sum) then
146 sum_tot = reproducing_sum(snm_re_raw(:,:,1:ltot), sums=snm_re(1:ltot), unscale=tmp_scale)
147 sum_tot = reproducing_sum(snm_im_raw(:,:,1:ltot), sums=snm_im(1:ltot), unscale=tmp_scale)
148 deallocate(snm_re_raw, snm_im_raw)
149 else
150 call sum_across_pes(snm_re, ltot)
151 call sum_across_pes(snm_im, ltot)
152 endif
153
154 if (id_clock_sht_global_sum>0) call cpu_clock_end(id_clock_sht_global_sum)
155 if (id_clock_sht_forward>0) call cpu_clock_end(id_clock_sht_forward)
156 if (id_clock_sht>0) call cpu_clock_end(id_clock_sht)
157end subroutine spherical_harmonics_forward
158
159!> Calculates inverse spherical harmonics transforms
160subroutine spherical_harmonics_inverse(G, CS, Snm_Re, Snm_Im, var, Nd)
161 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
162 type(sht_cs), intent(in) :: cs !< Control structure for SHT
163 real, intent(in) :: snm_re(:) !< SHT coefficients for the real modes (cosine)
164 !! in arbitrary units [a] or [A ~> a]
165 real, intent(in) :: snm_im(:) !< SHT coefficients for the imaginary modes (sine) in
166 !! the same arbitrary units as Snm_Re [a] or [A ~> a]
167 real, dimension(SZI_(G),SZJ_(G)), &
168 intent(out) :: var !< Output 2-D variable in the same arbitrary units
169 !! as Snm_Re and Snm_Im [a] or [A ~> a]
170 integer, optional, intent(in) :: nd !< Maximum degree of the spherical harmonics
171 !! overriding ndegree in the CS [nondim]
172 ! local variables
173 integer :: nmax ! Local copy of the maximum degree of the spherical harmonics [nondim]
174 real :: mfac ! A constant multiplier. mFac = 1 (if m==0) or 2 (if m>0) [nondim]
175 real, dimension(SZI_(G),SZJ_(G)) :: &
176 pmn, & ! Current associated Legendre polynomials of degree n and order m [nondim]
177 pmnm1, & ! Associated Legendre polynomials of degree n-1 and order m [nondim]
178 pmnm2 ! Associated Legendre polynomials of degree n-2 and order m [nondim]
179 integer :: i, j, k
180 integer :: is, ie, js, je, isd, ied, jsd, jed
181 integer :: m, n, l
182
183 if (.not.cs%initialized) call mom_error(fatal, "MOM_spherical_harmonics " // &
184 "spherical_harmonics_inverse: Module must be initialized before it is used.")
185
186 if (id_clock_sht>0) call cpu_clock_begin(id_clock_sht)
187 if (id_clock_sht_inverse>0) call cpu_clock_begin(id_clock_sht_inverse)
188
189 nmax = cs%ndegree ; if (present(nd)) nmax = nd
190
191 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
192 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
193
194 do j=jsd,jed ; do i=isd,ied
195 pmn(i,j) = 0.0 ; pmnm1(i,j) = 0.0 ; pmnm2(i,j) = 0.0
196 var(i,j) = 0.0
197 enddo ; enddo
198
199 do m=0,nmax
200 mfac = sign(1.0, m-0.5)*0.5 + 1.5
201 l = order2index(m, nmax)
202
203 do j=js,je ; do i=is,ie
204 var(i,j) = var(i,j) &
205 + mfac * cs%Pmm(i,j,m+1) * ( snm_re(l) * cs%cos_lonT(i,j,m+1) &
206 + snm_im(l) * cs%sin_lonT(i,j,m+1))
207 pmnm2(i,j) = 0.0
208 pmnm1(i,j) = cs%Pmm(i,j,m+1)
209 enddo ; enddo
210
211 do n=m+1,nmax ; do j=js,je ; do i=is,ie
212 pmn(i,j) = &
213 cs%a_recur(n+1,m+1) * cs%cos_clatT(i,j) * pmnm1(i,j) - cs%b_recur(n+1,m+1) * pmnm2(i,j)
214 var(i,j) = var(i,j) &
215 + mfac * pmn(i,j) * ( snm_re(l+n-m) * cs%cos_lonT(i,j,m+1) &
216 + snm_im(l+n-m) * cs%sin_lonT(i,j,m+1))
217 pmnm2(i,j) = pmnm1(i,j)
218 pmnm1(i,j) = pmn(i,j)
219 enddo ; enddo ; enddo
220 enddo
221
222 if (id_clock_sht_inverse>0) call cpu_clock_end(id_clock_sht_inverse)
223 if (id_clock_sht>0) call cpu_clock_end(id_clock_sht)
224end subroutine spherical_harmonics_inverse
225
226!> Calculate precomputed coefficients
227subroutine spherical_harmonics_init(G, param_file, CS)
228 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
229 type(param_file_type), intent(in) :: param_file !< A structure indicating
230 type(sht_cs), intent(inout) :: cs !< Control structure for spherical harmonic transforms
231
232 ! local variables
233 real, parameter :: pi = 4.0*atan(1.0) ! 3.1415926... calculated as 4*atan(1) [nondim]
234 real, parameter :: radian = pi / 180.0 ! Degree to Radian constant [radian degree-1]
235 real, dimension(SZI_(G),SZJ_(G)) :: sin_clatt ! sine of colatitude at the t-cells [nondim].
236 real :: pmm_coef ! = sqrt{ 1.0/(4.0*PI) * prod[(2k+1)/2k)] } [nondim].
237 integer :: is, ie, js, je
238 integer :: i, j, k
239 integer :: m, n
240 integer :: nd_sal ! Maximum degree for SAL
241 ! This include declares and sets the variable "version".
242# include "version_variable.h"
243 character(len=40) :: mdl = "MOM_spherical_harmonics" ! This module's name.
244
245 if (cs%initialized) return
246 cs%initialized = .true.
247
248 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
249
250 call log_version(param_file, mdl, version, "")
251 call get_param(param_file, mdl, "SAL_HARMONICS_DEGREE", nd_sal, "", default=0, do_not_log=.true.)
252 cs%ndegree = nd_sal
253 cs%lmax = calc_lmax(cs%ndegree)
254 call get_param(param_file, mdl, "SHT_REPRODUCING_SUM", cs%reprod_sum, &
255 "If true, use reproducing sums (invariant to PE layout) in inverse transform "// &
256 "of spherical harmonics. Otherwise use a simple sum of floating point numbers. ", &
257 default=.false.)
258
259 ! Calculate recurrence relationship coefficients
260 allocate(cs%a_recur(cs%ndegree+1, cs%ndegree+1), source=0.0)
261 allocate(cs%b_recur(cs%ndegree+1, cs%ndegree+1), source=0.0)
262 do m=0,cs%ndegree ; do n=m+1,cs%ndegree
263 ! These expressione will give NaNs with 32-bit integers for n > 23170, but this is trapped elsewhere.
264 cs%a_recur(n+1,m+1) = sqrt(real((2*n-1) * (2*n+1)) / real((n-m) * (n+m)))
265 cs%b_recur(n+1,m+1) = sqrt((real(2*n+1) * real((n+m-1) * (n-m-1))) / (real((n-m) * (n+m)) * real(2*n-3)))
266 enddo ; enddo
267
268 ! Calculate complex exponential factors
269 allocate(cs%cos_lonT_wtd(is:ie, js:je, cs%ndegree+1), source=0.0)
270 allocate(cs%sin_lonT_wtd(is:ie, js:je, cs%ndegree+1), source=0.0)
271 allocate(cs%cos_lonT(is:ie, js:je, cs%ndegree+1), source=0.0)
272 allocate(cs%sin_lonT(is:ie, js:je, cs%ndegree+1), source=0.0)
273 do m=0,cs%ndegree
274 do j=js,je ; do i=is,ie
275 cs%cos_lonT(i,j,m+1) = cos(real(m) * (g%geolonT(i,j)*radian))
276 cs%sin_lonT(i,j,m+1) = sin(real(m) * (g%geolonT(i,j)*radian))
277 cs%cos_lonT_wtd(i,j,m+1) = cs%cos_lonT(i,j,m+1) * g%areaT(i,j) / g%Rad_Earth_L**2
278 cs%sin_lonT_wtd(i,j,m+1) = cs%sin_lonT(i,j,m+1) * g%areaT(i,j) / g%Rad_Earth_L**2
279 enddo ; enddo
280 enddo
281
282 ! Calculate sine and cosine of colatitude
283 allocate(cs%cos_clatT(is:ie, js:je), source=0.0)
284 do j=js,je ; do i=is,ie
285 cs%cos_clatT(i,j) = cos(0.5*pi - g%geolatT(i,j)*radian)
286 sin_clatt(i,j) = sin(0.5*pi - g%geolatT(i,j)*radian)
287 enddo ; enddo
288
289 ! Calculate the diagonal elements of the associated Legendre polynomials (n=m)
290 allocate(cs%Pmm(is:ie,js:je,m+1), source=0.0)
291 do m=0,cs%ndegree
292 pmm_coef = 1.0/(4.0*pi)
293 do k=1,m ; pmm_coef = pmm_coef * (real(2*k+1) / real(2*k)) ; enddo
294 pmm_coef = sqrt(pmm_coef)
295 do j=js,je ; do i=is,ie
296 cs%Pmm(i,j,m+1) = pmm_coef * (sin_clatt(i,j)**m)
297 enddo ; enddo
298 enddo
299
300 id_clock_sht = cpu_clock_id('(Ocean spherical harmonics)', grain=clock_module)
301 id_clock_sht_forward = cpu_clock_id('(Ocean SHT forward)', grain=clock_routine)
302 id_clock_sht_inverse = cpu_clock_id('(Ocean SHT inverse)', grain=clock_routine)
303 id_clock_sht_global_sum = cpu_clock_id('(Ocean SHT global sum)', grain=clock_loop)
304
305end subroutine spherical_harmonics_init
306
307!> Deallocate any variables allocated in spherical_harmonics_init
308subroutine spherical_harmonics_end(CS)
309 type(sht_cs), intent(inout) :: cs !< Control structure for spherical harmonic transforms
310
311 deallocate(cs%cos_clatT)
312 deallocate(cs%Pmm)
313 deallocate(cs%cos_lonT_wtd, cs%sin_lonT_wtd, cs%cos_lonT, cs%sin_lonT)
314 deallocate(cs%a_recur, cs%b_recur)
315end subroutine spherical_harmonics_end
316
317!> Calculates the number of real elements (cosine) of spherical harmonics given maximum degree Nd.
318function calc_lmax(Nd) result(lmax)
319 integer :: lmax !< Number of real spherical harmonic modes [nondim]
320 integer, intent(in) :: nd !< Maximum degree [nondim]
321
322 lmax = (nd+2) * (nd+1) / 2
323end function calc_lmax
324
325!> Calculates the one-dimensional index number at (n=0, m=m), given order m and maximum degree Nd.
326!! It is sequenced with degree (n) changing first and order (m) changing second.
327function order2index(m, Nd) result(l)
328 integer :: l !< One-dimensional index number [nondim]
329 integer, intent(in) :: m !< Current order number [nondim]
330 integer, intent(in) :: nd !< Maximum degree [nondim]
331
332 l = ((nd+1) + (nd+1-(m-1)))*m/2 + 1
333end function order2index
334
335!> \namespace mom_spherical_harmonics
336!!
337!! \section section_spherical_harmonics Spherical harmonics
338!!
339!! This module contains the subroutines to calculate spherical harmonic transforms (SHT), namely, forward transform
340!! of a two-dimensional field into a given number of spherical harmonic modes and its inverse transform. This module
341!! is primarily used to but not limited to calculate self-attraction and loading (SAL) term, which is mostly relevant to
342!! high frequency motions such as tides. Should other needs arise in the future, this API can be easily modified.
343!! Currently, the transforms are for t-cell fields only.
344!!
345!! This module is stemmed from SAL calculation in Model for Prediction Across Scales (MPAS)-Ocean developed by Los
346!! Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. The algorithm
347!! for forward and inverse transforms loosely follows \cite Schaeffer2013.
348!!
349!! In forward transform, a two-dimensional physical field can be projected into a series of spherical harmonics. The
350!! spherical harmonic coefficient of degree n and order m for a field \f$f(\theta, \phi)\f$ is calculated as follows:
351!! \f[
352!! f^m_n = \int^{2\pi}_{0}\int^{\pi}_{0}f(\theta,\phi)Y^m_n(\theta,\phi)\sin\theta d\theta d\phi
353!! \f]
354!! and
355!! \f[
356!! Y^m_n(\theta,\phi) = P^m_n(\cos\theta)\exp(im\phi)
357!! \f]
358!! where \f$P^m_n(\cos \theta)\f$ is the normalized associated Legendre polynomial of degree n and order m. \f$\phi\f$
359!! is the longitude and \f$\theta\f$ is the colatitude.
360!! Or, written in the discretized form:
361!! \f[
362!! f^m_n = \sum^{Nj}_{0}\sum^{Ni}_{0}f(i,j)Y^m_n(i,j)A(i,j)/r_e^2
363!! \f]
364!! where \f$A\f$ is the area of the cell and \f$r_e\f$ is the radius of the Earth.
365!!
366!! In inverse transform, the first N degree spherical harmonic coefficients are used to reconstruct a two-dimensional
367!! physical field:
368!! \f[
369!! f(\theta,\phi) = \sum^N_{n=0}\sum^{n}_{m=-n}f^m_nY^m_n(\theta,\phi)
370!! \f]
371!!
372!! The exponential coefficients are pre-computed and stored in the memory. The associated Legendre polynomials are
373!! computed "on-the-fly", using the recurrence relationships to avoid large memory usage and take the advantage of
374!! array vectorization.
375!!
376!! The maximum degree of the spherical harmonics is a runtime parameter and the maximum used by all SHT applications.
377!! At the moment, it is only decided by <code>SAL_HARMONICS_DEGREE</code>.
378!!
379!! The forward transforms involve a global summation. Runtime flag <code>SHT_REPRODUCING_SUM</code> controls
380!! whether this is done in a bit-wise reproducing way or not.
381!!
382!! References:
383!!
384!! Barton, K.N., Pal, N., Brus, S.R., Petersen, M.R., Arbic, B.K., Engwirda, D., Roberts, A.F., Westerink, J.J.,
385!! Wirasaet, D. and Schindelegger, M., 2022. Global Barotropic Tide Modeling Using Inline Self‐Attraction and Loading in
386!! MPAS‐Ocean. Journal of Advances in Modeling Earth Systems, 14(11), p.e2022MS003207.
387!! https://doi.org/10.1029/2022MS003207
388!!
389!! Brus, S.R., Barton, K.N., Pal, N., Roberts, A.F., Engwirda, D., Petersen, M.R., Arbic, B.K., Wirasaet, D.,
390!! Westerink, J.J. and Schindelegger, M., 2023. Scalable self attraction and loading calculations for unstructured ocean
391!! tide models. Ocean Modelling, p.102160.
392!! https://doi.org/10.1016/j.ocemod.2023.102160
393!!
394!! Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations.
395!! Geochemistry, Geophysics, Geosystems, 14(3), pp.751-758.
396!! https://doi.org/10.1002/ggge.20071