MOM_unique_scales.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!> This module provides tools that can be used to check the uniqueness of the dimensional
6!! scaling factors used by the MOM6 ocean model or other models
8
10
11implicit none ; private
12
13! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
14! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
15! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
16! vary with the Boussinesq approximation, the Boussinesq variant is given first.
17
19
20contains
21
22!> This subroutine does a checks whether the provided dimensional scaling factors give a unique
23!! overall scaling for each of the combinations of units in description, and suggests a better
24!! combination if it is not unique. However, this subroutine does nothing if the verbosity level
25!! for this run is below 3.
26subroutine check_scaling_uniqueness(component, descs, weights, key, scales, max_powers)
27 character(len=*), intent(in) :: component !< The name of the component (e.g., MOM6) to use in messages
28 character(len=*), intent(in) :: descs(:) !< The descriptions for each combination of units
29 integer, intent(in) :: weights(:) !< A list of the weights for each described combination
30 character(len=*), intent(in) :: key(:) !< The key for the unit scaling
31 integer, intent(in) :: scales(:) !< The powers of 2 that give the scaling for each unit in key
32 integer, optional, intent(in) :: max_powers !< The maximum range of powers of 2 to search for
33 !! suggestions of better scaling factors, or 0 to avoid
34 !! suggesting improved factors.
35
36 ! Local variables
37 integer, dimension(size(key)) :: next_scales, prev_scales, better_scales
38 character(len=512) :: mesg
39 character(len=64) :: msg_frag
40 integer, dimension(size(key), size(weights)) :: list
41 integer :: verbosity
42 logical :: same_key
43 integer :: orig_cost, test_cost, better_cost, prev_cost ! Various squared-weight mismatch costs.
44 integer :: better_dp ! The absolute change in powers with the better estimate.
45 integer :: ndims, ns, m, n, i, p, itt, max_itt, max_pow
46
47 call assert((size(scales) == size(key)), "check_scaling_factors: Mismatched scales and key sizes.")
48 call assert((size(descs) == size(weights)), "check_scaling_factors: Mismatched descs and weights.")
49
50 verbosity = mom_get_verbosity()
51 ! Skip the rest of this routine if it would not write anything out.
52 if (verbosity < 3) return
53
54 ndims = size(key)
55 ns = size(weights)
56 max_pow = 0 ; if (present(max_powers)) max_pow = max_powers
57
58 list(:,:) = 0
59 do n=1,ns
60 call encode_dim_powers(descs(n), key, list(:,n))
61 enddo
62
63 if (verbosity >= 7) then
64 write(mesg, '(I0)') ns
65 call mom_mesg(trim(component)//": Extracted "//trim(mesg)//" unit combinations from the list.")
66 mesg = "Dim Key: ["
67 do i=1,ndims ; mesg = trim(mesg)//" "//trim(key(i)) ; enddo
68 mesg = trim(mesg)//"]:"
69 call mom_mesg(mesg)
70 do n=1,ns
71 call mom_mesg(trim(component)//": Extracted ["//trim(int_array_msg(list(:,n)))//"] from "//trim(descs(n)))
72 enddo
73
74 do n=1,ns ; do m=1,n-1
75 same_key = .true.
76 do i=1,ndims ; if (list(i,n) /= list(i,m)) same_key = .false. ; enddo
77 if (same_key) then
78 call mom_mesg(trim(component)//": The same powers occur for "//&
79 trim(descs(n))//" and "//trim(descs(m))//"." )
80 endif
81 enddo ; enddo
82 endif
83
84 orig_cost = non_unique_scales(scales, list, descs, weights, silent=(verbosity<4))
85
86 max_itt = 3*ndims ! Do up to 3 iterations for each rescalable dimension.
87 if (orig_cost /= 0) then
88 call mom_mesg(trim(component)//": The dimensional scaling factors are not unique.")
89 prev_cost = orig_cost
90 prev_scales(:) = scales(:)
91 do itt=1,max_itt
92 ! Iterate to find a better solution.
93 better_scales(:) = prev_scales(:)
94 better_cost = prev_cost
95 better_dp = 0
96 do i=1,ndims
97 if (scales(i) == 0) cycle ! DO not optimize unscaled dimensions.
98 next_scales(:) = prev_scales(:)
99 do p=-max_pow,max_pow
100 if ((p==0) .or. (p==prev_scales(i))) cycle
101 next_scales(i) = p
102 test_cost = non_unique_scales(next_scales, list, descs, weights, silent=.true.)
103 if ((test_cost < better_cost) .or. &
104 ((test_cost == better_cost) .and. (abs(p-prev_scales(i)) < better_dp))) then
105 ! This is a better scaling or has the same weighted mismatches but smaller
106 ! changes in rescaling factors, so it could be the next guess.
107 better_scales(:) = next_scales(:)
108 better_cost = test_cost
109 better_dp = abs(p - prev_scales(i))
110 endif
111 enddo
112 enddo
113 if (better_cost < prev_cost) then
114 ! Store the new best guess and try again.
115 prev_scales(:) = better_scales(:)
116 prev_cost = better_cost
117 else ! No further optimization is possible.
118 exit
119 endif
120 if (better_cost == 0) exit
121 if (verbosity >= 7) then
122 write(mesg, '("Iteration ",I0," scaling cost reduced from ",I0," with original scales to ", I0)') &
123 itt, orig_cost, better_cost
124 call mom_mesg(trim(component)//": "//trim(mesg)//" with revised scaling factors.")
125 endif
126 enddo
127 if (prev_cost < orig_cost) then
128 test_cost = non_unique_scales(prev_scales, list, descs, weights, silent=(verbosity<4))
129 mesg = trim(component)//": Suggested improved scales: "
130 do i=1,ndims ; if ((prev_scales(i) /= scales(i)) .and. (scales(i) /= 0)) then
131 write(msg_frag, '(I0)') prev_scales(i)
132 mesg = trim(mesg)//" "//trim(key(i))//"_RESCALE_POWER = "//trim(msg_frag)
133 endif ; enddo
134 call mom_mesg(mesg)
135
136 write(mesg, '(I0)') orig_cost
137 write(msg_frag, '(I0)') test_cost
138 mesg = trim(component)//": Scaling overlaps reduced from "//trim(mesg)//&
139 " with original scales to "//trim(msg_frag)//" with suggested scales."
140 call mom_mesg(mesg)
141 endif
142
143 endif
144
145end subroutine check_scaling_uniqueness
146
147!> Convert a unit scaling descriptor into an array of the dimensions of powers given in the key
148subroutine encode_dim_powers(scaling, key, dim_powers)
149
150 character(len=*), intent(in) :: scaling !< The unit description that will be converted
151 character(len=*), dimension(:), intent(in) :: key(:) !< The key for the unit scaling
152 integer, dimension(size(key)), intent(out) :: dim_powers !< The dimensions in scaling of each
153 !! element of they key.
154
155 ! Local variables
156 character(len=:), allocatable :: actstr ! The full active remaining string to be parsed.
157 character(len=:), allocatable :: fragment ! The space-delimited fragment being parsed.
158 character(len=:), allocatable :: dimnm ! The probable dimension name
159 character(len=11) :: numbers ! The list of characters that could make up the exponent.
160 ! character(len=128) :: mesg
161 integer :: istart, iend, ieq, ief, ipow ! Positions in strings.
162 integer :: dp ! The power for this dimension.
163 integer :: ndim ! The number of dimensional scaling factors to consider.
164 integer :: n
165
166 dim_powers(:) = 0
167
168 iend = index(scaling, "~>") - 1
169 if (iend < 1) return
170
171 ! Parse the key.
172 ndim = size(key)
173 numbers = "-0123456789"
174
175 ! Strip away any leading square brace.
176 istart = index(scaling(:iend), "[") + 1
177 ! If there is an "=" in the string, start after this.
178 ieq = index(scaling(istart:iend), "=", back=.true.)
179 if (ieq > 0) istart = istart + ieq
180
181 ! Set up the active string to work on.
182 actstr = trim(adjustl(scaling(istart:iend)))
183 do ! Loop over each of the elements in the unit scaling descriptor.
184 if (len_trim(actstr) == 0) exit
185 ief = index(actstr, " ") - 1
186 if (ief <= 0) ief = len_trim(actstr)
187 fragment = actstr(:ief)
188 ipow = scan(fragment, "-")
189 if (ipow == 0) ipow = scan(fragment, numbers)
190
191 if (ipow == 0) then ! There is no exponent
192 dimnm = fragment
193 dp = 1
194 ! call MOM_mesg("Parsing powerless fragment "//trim(fragment)//" from "//trim(scaling))
195 else
196 if (verify(fragment(ipow:), numbers) == 0) then
197 read(fragment(ipow:),*) dp
198 dimnm = fragment(:ipow-1)
199 ! write(mesg, '(I0)') dp
200 ! call MOM_mesg("Parsed fragment "//trim(fragment)//" from "//trim(scaling)//&
201 ! " as "//trim(dimnm)//trim(mesg))
202 else
203 dimnm = fragment
204 dp = 1
205 ! call MOM_mesg("Unparsed fragment "//trim(fragment)//" from "//trim(scaling))
206 endif
207 endif
208
209 do n=1,ndim
210 if (trim(dimnm) == trim(key(n))) then
211 dim_powers(n) = dim_powers(n) + dp
212 exit
213 endif
214 enddo
215
216 ! Remove the leading fragment that has been parsed from actstr
217 actstr = trim(adjustl(actstr(ief+1:)))
218 enddo
219
220end subroutine encode_dim_powers
221
222!> Find the integer power of two that describe each of the scaling factors, or return 0 for
223!! scaling factors that are not exceptionally close to an integer power of 2.
224subroutine scales_to_powers(scale, pow2)
225 real, intent(in) :: scale(:) !< The scaling factor for each dimension
226 integer, intent(out) :: pow2(:) !< The exact powers of 2 for each scale, or 0 for non-exact powers of 2.
227
228 real :: log2_sc ! The log base 2 of an element of scale
229 integer :: n, ndim
230
231 ndim = size(scale)
232
233 ! Find the integer power of two for the scaling factors, but skip the analysis of any factors
234 ! that are not close enough to being integer powers of 2.
235 do n=1,ndim
236 if (abs(scale(n)) > 0.0) then
237 log2_sc = log(abs(scale(n))) / log(2.0)
238 else
239 log2_sc = 0.0
240 endif
241 if (abs(log2_sc - nint(log2_sc)) < 1.0e-6) then
242 ! This is close to an integer power of two.
243 pow2(n) = nint(log2_sc)
244 else
245 ! This is not being scaled by an integer power of 2, so return 0.
246 pow2(n) = 0
247 endif
248 enddo
249
250end subroutine scales_to_powers
251
252!> Determine from the list of scaling factors and the unit combinations that are in use whether
253!! all these combinations scale uniquely.
254integer function non_unique_scales(scales, list, descs, weights, silent)
255 integer, intent(in) :: scales(:) !< The power of 2 that gives the scaling factor for each dimension
256 integer, intent(in) :: list(:,:) !< A list of the integers for each scaling
257 character(len=*), intent(in) :: descs(:) !< The unit descriptions that have been converted
258 integer, intent(in) :: weights(:) !< A list of the weights for each scaling
259 logical, optional, intent(in) :: silent !< If present and true, do not write any output.
260
261 ! Local variables
262 integer, dimension(size(weights)) :: res_pow ! The net rescaling power for each combination.
263 integer, dimension(size(weights)) :: wt_merge ! The merged weights of scaling factors with common powers
264 ! for the dimensions being tested.
265 logical :: same_key, same_scales, verbose
266 character(len=256) :: mesg
267 integer :: nonzero_count ! The number of non-zero scaling factors
268 integer :: ndim ! The number of dimensional scaling factors to work with
269 integer :: i, n, m, ns
270
271 non_unique_scales = -9999 ! Set return value to a _dummy_ value
272
273 verbose = .true. ; if (present(silent)) verbose = .not.silent
274
275 ndim = size(scales)
276 ns = size(descs)
277 call assert((size(scales) == size(list, 1)), "non_unique_scales: Mismatched scales and list sizes.")
278 call assert((size(descs) == size(list, 2)), "non_unique_scales: Mismatched descs and list sizes.")
279 call assert((size(descs) == size(weights)), "non_unique_scales: Mismatched descs and weights.")
280
281 ! Return .true. if all scaling powers are 0, or there is only one scaling factor in use.
282 nonzero_count = 0 ; do n=1,ndim ; if (scales(n) /= 0) nonzero_count = nonzero_count + 1 ; enddo
283 if (nonzero_count <= 1) return
284
285 ! Figure out which unit combinations are unique for the set of dimensions and scaling factors
286 ! that are being tested, and combine the weights for scaling factors.
287 wt_merge(:) = weights(:)
288 do n=1,ns ; do m=1,n-1
289 same_key = .true.
290 same_scales = .true.
291 do i=1,ndim
292 if (list(i,n) /= list(i,m)) same_key = .false.
293 if ((scales(i) /= 0) .and. (list(i,n) /= list(i,m))) same_scales = .false.
294 enddo
295 if (same_key .or. same_scales) then
296 if (wt_merge(n) > wt_merge(m)) then
297 wt_merge(n) = wt_merge(n) + wt_merge(m)
298 wt_merge(m) = 0
299 else
300 wt_merge(m) = wt_merge(m) + wt_merge(n)
301 wt_merge(n) = 0
302 endif
303 endif
304 if (wt_merge(n) == 0) exit ! Go to the next value of n.
305 enddo ; enddo
306
307 do n=1,ns
308 res_pow(n) = 0
309 do i=1,ndim
310 res_pow(n) = res_pow(n) + scales(i) * list(i,n)
311 enddo
312 enddo
313
314 ! Determine the weighted cost of non-unique scaling factors.
316 do n=1,ns ; if (wt_merge(n) > 0) then ; do m=1,n-1 ; if (wt_merge(m) > 0) then
317 if (res_pow(n) == res_pow(m)) then
318 ! Use the product of the weights as the cost, as this should be vaguely proportional to
319 ! the likelihood that these factors would be combined in an expression.
320 non_unique_scales = min(non_unique_scales + wt_merge(n) * wt_merge(m), 99999999)
321 if (verbose) then
322 write(mesg, '(I0)') res_pow(n)
323 call mom_mesg("The factors "//trim(descs(n))//" and "//trim(descs(m))//" both scale to "//&
324 trim(mesg)//" for the given powers.")
325
326 ! call MOM_mesg("Powers ["//trim(int_array_msg(list(:,n)))//"] and ["//&
327 ! trim(int_array_msg(list(:,m)))//"] with rescaling by ["//&
328 ! trim(int_array_msg(scales))//"]")
329 endif
330 endif
331 endif ; enddo ; endif ; enddo
332
333end function non_unique_scales
334
335!> Return a string the elements of an array of integers
336function int_array_msg(array)
337 integer, intent(in) :: array(:) !< The array whose values are to be written.
338 character(len=16*size(array)) :: int_array_msg
339
340 character(len=12) :: msg_frag
341 integer :: i, ni
342 ni = size(array)
343
344 int_array_msg = ""
345 if (ni < 1) return
346
347 do i=1,ni
348 write(msg_frag, '(I0)') array(i)
349 if (i == 1) then
350 int_array_msg = trim(msg_frag)
351 else
352 int_array_msg = trim(int_array_msg)//" "//trim(msg_frag)
353 endif
354 enddo
355end function int_array_msg
356
357end module mom_unique_scales