@@ -29,8 +29,6 @@ module MOM_spherical_harmonics
29
29
sin_lonT_wtd(:,:,:) ! < Precomputed area-weighted sine factors at the t-cells [nondim]
30
30
real , allocatable :: a_recur(:,:), & ! < Precomputed recurrence coefficients a [nondim].
31
31
b_recur(:,:) ! < Precomputed recurrence coefficients b [nondim].
32
- real , allocatable :: Snm_Re_raw(:,:,:), & ! < Array to store un-summed SHT coefficients
33
- Snm_Im_raw(:,:,:) ! < at the t-cells for reproducing sums [same as input variable]
34
32
logical :: reprod_sum ! < True if use reproducible global sums
35
33
end type sht_CS
36
34
@@ -46,9 +44,13 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale
46
44
type (ocean_grid_type), intent (in ) :: G ! < The ocean's grid structure.
47
45
type (sht_CS), intent (inout ) :: CS ! < Control structure for SHT
48
46
real , dimension (SZI_(G),SZJ_(G)), &
49
- intent (in ) :: var ! < Input 2-D variable [A]
50
- real , intent (out ) :: Snm_Re(:) ! < SHT coefficients for the real modes (cosine) [A]
51
- real , intent (out ) :: Snm_Im(:) ! < SHT coefficients for the imaginary modes (sine) [A]
47
+ intent (in ) :: var ! < Input 2-D variable in arbitrary mks units [a]
48
+ ! ! or in arbitrary rescaled units [A ~> a] if
49
+ ! ! tmp_scale is present
50
+ real , intent (out ) :: Snm_Re(:) ! < SHT coefficients for the real modes (cosine) in
51
+ ! ! the same arbitrary units as var [a] or [A ~> a]
52
+ real , intent (out ) :: Snm_Im(:) ! < SHT coefficients for the imaginary modes (sine) in
53
+ ! ! the same arbitrary units as var [a] or [A ~> a]
52
54
integer , optional , intent (in ) :: Nd ! < Maximum degree of the spherical harmonics
53
55
! ! overriding ndegree in the CS [nondim]
54
56
real , optional , intent (in ) :: tmp_scale ! < A temporary rescaling factor to convert
@@ -61,10 +63,13 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale
61
63
pmn, & ! Current associated Legendre polynomials of degree n and order m [nondim]
62
64
pmnm1, & ! Associated Legendre polynomials of degree n-1 and order m [nondim]
63
65
pmnm2 ! Associated Legendre polynomials of degree n-2 and order m [nondim]
64
- real :: scale ! A rescaling factor to temporarily convert var to MKS units during the
65
- ! reproducing sums [a A-1 ~> 1]
66
- real :: I_scale ! The inverse of scale [A a-1 ~> 1]
67
- real :: sum_tot ! The total of all components output by the reproducing sum in arbitrary units [a]
66
+ real , allocatable , dimension (:,:,:) :: &
67
+ Snm_Re_raw, & ! Array of un-summed real spherical harmonics transform coefficients for
68
+ ! reproducing sums in the same arbitrary units as var, [a] or [A ~> a]
69
+ Snm_Im_raw ! Array of un-summed imaginary spherical harmonics transform coefficients for
70
+ ! reproducing sums in the same arbitrary units as var, [a] or [A ~> a]
71
+ real :: sum_tot ! The total of all components output by the reproducing sum in the same
72
+ ! arbitrary units as var, [a] or [A ~> a]
68
73
integer :: i, j, k
69
74
integer :: is, ie, js, je, isd, ied, jsd, jed
70
75
integer :: m, n, l
@@ -75,35 +80,36 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale
75
80
if (id_clock_sht> 0 ) call cpu_clock_begin(id_clock_sht)
76
81
if (id_clock_sht_forward> 0 ) call cpu_clock_begin(id_clock_sht_forward)
77
82
78
- Nmax = CS% ndegree; if (present (Nd)) Nmax = Nd
83
+ Nmax = CS% ndegree ; if (present (Nd)) Nmax = Nd
79
84
Ltot = calc_lmax(Nmax)
80
85
81
86
is = G% isc ; ie = G% iec ; js = G% jsc ; je = G% jec
82
87
isd = G% isd ; ied = G% ied ; jsd = G% jsd ; jed = G% jed
83
88
84
89
do j= jsd,jed ; do i= isd,ied
85
- pmn(i,j) = 0.0 ; pmnm1(i,j) = 0.0 ; pmnm2(i,j) = 0.0
90
+ pmn(i,j) = 0.0 ; pmnm1(i,j) = 0.0 ; pmnm2(i,j) = 0.0
86
91
enddo ; enddo
87
92
88
- do l= 1 ,Ltot ; Snm_Re(l) = 0.0 ; Snm_Im(l) = 0.0 ; enddo
93
+ do l= 1 ,Ltot ; Snm_Re(l) = 0.0 ; Snm_Im(l) = 0.0 ; enddo
89
94
90
95
if (CS% reprod_sum) then
91
- scale = 1.0 ; if (present (tmp_scale)) scale = tmp_scale
96
+ allocate (Snm_Re_raw(is:ie, js:je, Ltot), source= 0.0 )
97
+ allocate (Snm_Im_raw(is:ie, js:je, Ltot), source= 0.0 )
92
98
do m= 0 ,Nmax
93
99
l = order2index(m, Nmax)
94
100
95
101
do j= js,je ; do i= is,ie
96
- CS % Snm_Re_raw(i,j,l) = (scale * var(i,j) ) * CS% Pmm(i,j,m+1 ) * CS% cos_lonT_wtd(i,j,m+1 )
97
- CS % Snm_Im_raw(i,j,l) = (scale * var(i,j) ) * CS% Pmm(i,j,m+1 ) * CS% sin_lonT_wtd(i,j,m+1 )
102
+ Snm_Re_raw(i,j,l) = var(i,j) * CS% Pmm(i,j,m+1 ) * CS% cos_lonT_wtd(i,j,m+1 )
103
+ Snm_Im_raw(i,j,l) = var(i,j) * CS% Pmm(i,j,m+1 ) * CS% sin_lonT_wtd(i,j,m+1 )
98
104
pmnm2(i,j) = 0.0
99
105
pmnm1(i,j) = CS% Pmm(i,j,m+1 )
100
106
enddo ; enddo
101
107
102
108
do n = m+1 , Nmax ; do j= js,je ; do i= is,ie
103
109
pmn(i,j) = &
104
110
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)
105
- CS % Snm_Re_raw(i,j,l+ n- m) = (scale * var(i,j) ) * pmn(i,j) * CS% cos_lonT_wtd(i,j,m+1 )
106
- CS % Snm_Im_raw(i,j,l+ n- m) = (scale * var(i,j) ) * pmn(i,j) * CS% sin_lonT_wtd(i,j,m+1 )
111
+ Snm_Re_raw(i,j,l+ n- m) = var(i,j) * pmn(i,j) * CS% cos_lonT_wtd(i,j,m+1 )
112
+ Snm_Im_raw(i,j,l+ n- m) = var(i,j) * pmn(i,j) * CS% sin_lonT_wtd(i,j,m+1 )
107
113
pmnm2(i,j) = pmnm1(i,j)
108
114
pmnm1(i,j) = pmn(i,j)
109
115
enddo ; enddo ; enddo
@@ -133,15 +139,9 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale
133
139
if (id_clock_sht_global_sum> 0 ) call cpu_clock_begin(id_clock_sht_global_sum)
134
140
135
141
if (CS% reprod_sum) then
136
- sum_tot = reproducing_sum(CS% Snm_Re_raw(:,:,1 :Ltot), sums= Snm_Re(1 :Ltot))
137
- sum_tot = reproducing_sum(CS% Snm_Im_raw(:,:,1 :Ltot), sums= Snm_Im(1 :Ltot))
138
- if (scale /= 1.0 ) then
139
- I_scale = 1.0 / scale
140
- do l= 1 ,Ltot
141
- Snm_Re(l) = I_scale * Snm_Re(l)
142
- Snm_Im(l) = I_scale * Snm_Im(l)
143
- enddo
144
- endif
142
+ sum_tot = reproducing_sum(Snm_Re_raw(:,:,1 :Ltot), sums= Snm_Re(1 :Ltot), unscale= tmp_scale)
143
+ sum_tot = reproducing_sum(Snm_Im_raw(:,:,1 :Ltot), sums= Snm_Im(1 :Ltot), unscale= tmp_scale)
144
+ deallocate (Snm_Re_raw, Snm_Im_raw)
145
145
else
146
146
call sum_across_PEs(Snm_Re, Ltot)
147
147
call sum_across_PEs(Snm_Im, Ltot)
@@ -156,10 +156,13 @@ end subroutine spherical_harmonics_forward
156
156
subroutine spherical_harmonics_inverse (G , CS , Snm_Re , Snm_Im , var , Nd )
157
157
type (ocean_grid_type), intent (in ) :: G ! < The ocean's grid structure.
158
158
type (sht_CS), intent (in ) :: CS ! < Control structure for SHT
159
- real , intent (in ) :: Snm_Re(:) ! < SHT coefficients for the real modes (cosine) [A]
160
- real , intent (in ) :: Snm_Im(:) ! < SHT coefficients for the imaginary modes (sine) [A]
159
+ real , intent (in ) :: Snm_Re(:) ! < SHT coefficients for the real modes (cosine)
160
+ ! ! in arbitrary units [a] or [A ~> a]
161
+ real , intent (in ) :: Snm_Im(:) ! < SHT coefficients for the imaginary modes (sine) in
162
+ ! ! the same arbitrary units as Snm_Re [a] or [A ~> a]
161
163
real , dimension (SZI_(G),SZJ_(G)), &
162
- intent (out ) :: var ! < Output 2-D variable [A]
164
+ intent (out ) :: var ! < Output 2-D variable in the same arbitrary units
165
+ ! ! as Snm_Re and Snm_Im [a] or [A ~> a]
163
166
integer , optional , intent (in ) :: Nd ! < Maximum degree of the spherical harmonics
164
167
! ! overriding ndegree in the CS [nondim]
165
168
! local variables
@@ -179,13 +182,13 @@ subroutine spherical_harmonics_inverse(G, CS, Snm_Re, Snm_Im, var, Nd)
179
182
if (id_clock_sht> 0 ) call cpu_clock_begin(id_clock_sht)
180
183
if (id_clock_sht_inverse> 0 ) call cpu_clock_begin(id_clock_sht_inverse)
181
184
182
- Nmax = CS% ndegree; if (present (Nd)) Nmax = Nd
185
+ Nmax = CS% ndegree ; if (present (Nd)) Nmax = Nd
183
186
184
187
is = G% isc ; ie = G% iec ; js = G% jsc ; je = G% jec
185
188
isd = G% isd ; ied = G% ied ; jsd = G% jsd ; jed = G% jed
186
189
187
190
do j= jsd,jed ; do i= isd,ied
188
- pmn(i,j) = 0.0 ; pmnm1(i,j) = 0.0 ; pmnm2(i,j) = 0.0
191
+ pmn(i,j) = 0.0 ; pmnm1(i,j) = 0.0 ; pmnm2(i,j) = 0.0
189
192
var(i,j) = 0.0
190
193
enddo ; enddo
191
194
@@ -250,19 +253,19 @@ subroutine spherical_harmonics_init(G, param_file, CS)
250
253
default= .False. )
251
254
252
255
! Calculate recurrence relationship coefficients
253
- allocate (CS% a_recur(CS% ndegree+1 , CS% ndegree+1 )); CS % a_recur(:,:) = 0.0
254
- allocate (CS% b_recur(CS% ndegree+1 , CS% ndegree+1 )); CS % b_recur(:,:) = 0.0
256
+ allocate (CS% a_recur(CS% ndegree+1 , CS% ndegree+1 ), source = 0.0 )
257
+ allocate (CS% b_recur(CS% ndegree+1 , CS% ndegree+1 ), source = 0.0 )
255
258
do m= 0 ,CS% ndegree ; do n= m+1 ,CS% ndegree
256
259
! These expressione will give NaNs with 32-bit integers for n > 23170, but this is trapped elsewhere.
257
260
CS% a_recur(n+1 ,m+1 ) = sqrt (real ((2 * n-1 ) * (2 * n+1 )) / real ((n- m) * (n+ m)))
258
261
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 )))
259
262
enddo ; enddo
260
263
261
264
! Calculate complex exponential factors
262
- allocate (CS% cos_lonT_wtd(is:ie, js:je, CS% ndegree+1 )); CS % cos_lonT_wtd(:,:,:) = 0.0
263
- allocate (CS% sin_lonT_wtd(is:ie, js:je, CS% ndegree+1 )); CS % sin_lonT_wtd(:,:,:) = 0.0
264
- allocate (CS% cos_lonT(is:ie, js:je, CS% ndegree+1 )); CS % cos_lonT(:,:,:) = 0.0
265
- allocate (CS% sin_lonT(is:ie, js:je, CS% ndegree+1 )); CS % sin_lonT(:,:,:) = 0.0
265
+ allocate (CS% cos_lonT_wtd(is:ie, js:je, CS% ndegree+1 ), source = 0.0 )
266
+ allocate (CS% sin_lonT_wtd(is:ie, js:je, CS% ndegree+1 ), source = 0.0 )
267
+ allocate (CS% cos_lonT(is:ie, js:je, CS% ndegree+1 ), source = 0.0 )
268
+ allocate (CS% sin_lonT(is:ie, js:je, CS% ndegree+1 ), source = 0.0 )
266
269
do m= 0 ,CS% ndegree
267
270
do j= js,je ; do i= is,ie
268
271
CS% cos_lonT(i,j,m+1 ) = cos (real (m) * (G% geolonT(i,j)* RADIAN))
@@ -273,28 +276,23 @@ subroutine spherical_harmonics_init(G, param_file, CS)
273
276
enddo
274
277
275
278
! Calculate sine and cosine of colatitude
276
- allocate (CS% cos_clatT(is:ie, js:je)); CS % cos_clatT(:,:) = 0.0
279
+ allocate (CS% cos_clatT(is:ie, js:je), source = 0.0 )
277
280
do j= js,je ; do i= is,ie
278
281
CS% cos_clatT(i,j) = cos (0.5 * PI - G% geolatT(i,j)* RADIAN)
279
282
sin_clatT(i,j) = sin (0.5 * PI - G% geolatT(i,j)* RADIAN)
280
283
enddo ; enddo
281
284
282
285
! Calculate the diagonal elements of the associated Legendre polynomials (n=m)
283
- allocate (CS% Pmm(is:ie,js:je,m+1 )); CS % Pmm(:,:,:) = 0.0
286
+ allocate (CS% Pmm(is:ie,js:je,m+1 ), source = 0.0 )
284
287
do m= 0 ,CS% ndegree
285
288
Pmm_coef = 1.0 / (4.0 * PI)
286
- do k= 1 ,m ; Pmm_coef = Pmm_coef * (real (2 * k+1 ) / real (2 * k)); enddo
289
+ do k= 1 ,m ; Pmm_coef = Pmm_coef * (real (2 * k+1 ) / real (2 * k)) ; enddo
287
290
Pmm_coef = sqrt (Pmm_coef)
288
291
do j= js,je ; do i= is,ie
289
292
CS% Pmm(i,j,m+1 ) = Pmm_coef * (sin_clatT(i,j)** m)
290
293
enddo ; enddo
291
294
enddo
292
295
293
- if (CS% reprod_sum) then
294
- allocate (CS% Snm_Re_raw(is:ie, js:je, CS% lmax)); CS% Snm_Re_raw = 0.0
295
- allocate (CS% Snm_Im_raw(is:ie, js:je, CS% lmax)); CS% Snm_Im_raw = 0.0
296
- endif
297
-
298
296
id_clock_sht = cpu_clock_id(' (Ocean spherical harmonics)' , grain= CLOCK_MODULE)
299
297
id_clock_sht_forward = cpu_clock_id(' (Ocean SHT forward)' , grain= CLOCK_ROUTINE)
300
298
id_clock_sht_inverse = cpu_clock_id(' (Ocean SHT inverse)' , grain= CLOCK_ROUTINE)
@@ -310,8 +308,6 @@ subroutine spherical_harmonics_end(CS)
310
308
deallocate (CS% Pmm)
311
309
deallocate (CS% cos_lonT_wtd, CS% sin_lonT_wtd, CS% cos_lonT, CS% sin_lonT)
312
310
deallocate (CS% a_recur, CS% b_recur)
313
- if (CS% reprod_sum) &
314
- deallocate (CS% Snm_Re_raw, CS% Snm_Im_raw)
315
311
end subroutine spherical_harmonics_end
316
312
317
313
! > Calculates the number of real elements (cosine) of spherical harmonics given maximum degree Nd.
0 commit comments