MOM_porous_barriers.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!> Module for calculating curve fit for porous topography.
6!written by sjd
8
9use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_module
10use mom_error_handler, only : mom_error, fatal
11use mom_grid, only : ocean_grid_type
15use mom_interface_heights, only : find_eta
16use mom_time_manager, only : time_type
18use mom_file_parser, only : param_file_type, get_param, log_version
20use mom_debugging, only : hchksum, uvchksum
21
22implicit none ; private
23
25
26#include <MOM_memory.h>
27
28!> The control structure for the MOM_porous_barriers module
29type, public :: porous_barrier_cs ; private
30 logical :: initialized = .false. !< True if this control structure has been initialized.
31 type(diag_ctrl), pointer :: &
32 diag => null() !< A structure to regulate diagnostic output timing
33 logical :: debug !< If true, write verbose checksums for debugging purposes.
34 real :: mask_depth !< The depth shallower than which porous barrier is not applied [Z ~> m]
35 integer :: eta_interp !< An integer indicating how the interface heights at the velocity
36 !! points are calculated. Valid values are given by the parameters
37 !! defined below: MAX, MIN, ARITHMETIC and HARMONIC.
38 integer :: answer_date !< The vintage of the porous barrier weight function calculations.
39 !! Values below 20220806 recover the old answers in which the layer
40 !! averaged weights are not strictly limited by an upper-bound of 1.0 .
41 !>@{ Diagnostic IDs
42 integer :: id_por_layer_widthu = -1, id_por_layer_widthv = -1, &
43 id_por_face_areau = -1, id_por_face_areav = -1
44 !>@}
45end type porous_barrier_cs
46
47integer :: id_clock_porous_barrier !< CPU clock for porous barrier
48
49!>@{ Enumeration values for eta interpolation schemes
50integer, parameter :: eta_interp_max = 1
51integer, parameter :: eta_interp_min = 2
52integer, parameter :: eta_interp_arith = 3
53integer, parameter :: eta_interp_harm = 4
54character(len=20), parameter :: eta_interp_max_string = "MAX"
55character(len=20), parameter :: eta_interp_min_string = "MIN"
56character(len=20), parameter :: eta_interp_arith_string = "ARITHMETIC"
57character(len=20), parameter :: eta_interp_harm_string = "HARMONIC"
58!>@}
59
60contains
61
62!> subroutine to assign porous barrier widths averaged over a layer
63subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt)
64 ! Note: eta_bt is not currently used
65 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
66 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
67 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
68 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
69 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
70 !! thermodynamic variables.
71 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable
72 !! used to dilate the layer thicknesses
73 !! [H ~> m or kg m-2].
74 type(porous_barrier_type), intent(inout) :: pbv !< porous barrier fractional cell metrics
75 type(porous_barrier_cs), intent(in) :: cs !< Control structure for porous barrier
76
77 !local variables
78 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_u ! Layer interface heights at u points [Z ~> m]
79 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_v ! Layer interface heights at v points [Z ~> m]
80 real, dimension(SZIB_(G),SZJB_(G)) :: a_layer_prev ! Integral of fractional open width from the bottom
81 ! to the previous layer at u or v points [Z ~> m]
82 logical, dimension(SZIB_(G),SZJB_(G)) :: do_i ! Booleans for calculation at u or v points
83 ! updated while moving up layers
84 real :: a_layer ! Integral of fractional open width from bottom to current layer [Z ~> m]
85 real :: dz_min ! The minimum layer thickness [Z ~> m]
86 real :: dmask ! The depth below which porous barrier is not applied [Z ~> m]
87 integer :: i, j, k, nk, is, ie, js, je, isq, ieq, jsq, jeq
88
89 if (.not.cs%initialized) call mom_error(fatal, &
90 "MOM_Porous_barrier: Module must be initialized before it is used.")
91
92 call cpu_clock_begin(id_clock_porous_barrier)
93
94 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nk = gv%ke
95 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
96
97 if (cs%answer_date < 20220806) then
98 dmask = 0.0
99 else
100 dmask = cs%mask_depth
101 endif
102
103 call calc_eta_at_uv(eta_u, eta_v, cs%eta_interp, dmask, h, tv, g, gv, us)
104
105 dz_min = gv%Angstrom_Z
106
107 ! u-points
108 do j=js,je ; do i=isq,ieq ; do_i(i,j) = .false. ; enddo ; enddo
109
110 do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
111 call calc_por_layer(g%porous_DminU(i,j), g%porous_DmaxU(i,j), g%porous_DavgU(i,j), &
112 eta_u(i,j,nk+1), a_layer_prev(i,j), do_i(i,j))
113 endif ; enddo ; enddo
114
115 if (cs%answer_date < 20220806) then
116 do k=nk,1,-1 ; do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
117 call calc_por_layer(g%porous_DminU(i,j), g%porous_DmaxU(i,j), g%porous_DavgU(i,j), &
118 eta_u(i,j,k), a_layer, do_i(i,j))
119 if (eta_u(i,j,k) - eta_u(i,j,k+1) > 0.0) then
120 pbv%por_face_areaU(i,j,k) = (a_layer - a_layer_prev(i,j)) / (eta_u(i,j,k) - eta_u(i,j,k+1))
121 else
122 pbv%por_face_areaU(i,j,k) = 0.0
123 endif
124 a_layer_prev(i,j) = a_layer
125 endif ; enddo ; enddo ; enddo
126 else
127 do k=nk,1,-1 ; do j=js,je ; do i=isq,ieq
128 if (do_i(i,j)) then
129 call calc_por_layer(g%porous_DminU(i,j), g%porous_DmaxU(i,j), g%porous_DavgU(i,j), &
130 eta_u(i,j,k), a_layer, do_i(i,j))
131 if (eta_u(i,j,k) - (eta_u(i,j,k+1)+dz_min) > 0.0) then
132 pbv%por_face_areaU(i,j,k) = min(1.0, (a_layer - a_layer_prev(i,j)) / (eta_u(i,j,k) - eta_u(i,j,k+1)))
133 else
134 pbv%por_face_areaU(i,j,k) = 0.0 ! use calc_por_interface() might be a better choice
135 endif
136 a_layer_prev(i,j) = a_layer
137 else
138 pbv%por_face_areaU(i,j,k) = 1.0
139 endif
140 enddo ; enddo ; enddo
141 endif
142
143 ! v-points
144 do j=jsq,jeq ; do i=is,ie ; do_i(i,j) = .false. ; enddo ; enddo
145
146 do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
147 call calc_por_layer(g%porous_DminV(i,j), g%porous_DmaxV(i,j), g%porous_DavgV(i,j), &
148 eta_v(i,j,nk+1), a_layer_prev(i,j), do_i(i,j))
149 endif ; enddo ; enddo
150
151 if (cs%answer_date < 20220806) then
152 do k=nk,1,-1 ; do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
153 call calc_por_layer(g%porous_DminV(i,j), g%porous_DmaxV(i,j), g%porous_DavgV(i,j), &
154 eta_v(i,j,k), a_layer, do_i(i,j))
155 if (eta_v(i,j,k) - eta_v(i,j,k+1) > 0.0) then
156 pbv%por_face_areaV(i,j,k) = (a_layer - a_layer_prev(i,j)) / (eta_v(i,j,k) - eta_v(i,j,k+1))
157 else
158 pbv%por_face_areaV(i,j,k) = 0.0
159 endif
160 a_layer_prev(i,j) = a_layer
161 endif ; enddo ; enddo ; enddo
162 else
163 do k=nk,1,-1 ; do j=jsq,jeq ; do i=is,ie
164 if (do_i(i,j)) then
165 call calc_por_layer(g%porous_DminV(i,j), g%porous_DmaxV(i,j), g%porous_DavgV(i,j), &
166 eta_v(i,j,k), a_layer, do_i(i,j))
167 if (eta_v(i,j,k) - (eta_v(i,j,k+1)+dz_min) > 0.0) then
168 pbv%por_face_areaV(i,j,k) = min(1.0, (a_layer - a_layer_prev(i,j)) / (eta_v(i,j,k) - eta_v(i,j,k+1)))
169 else
170 pbv%por_face_areaV(i,j,k) = 0.0 ! use calc_por_interface() might be a better choice
171 endif
172 a_layer_prev(i,j) = a_layer
173 else
174 pbv%por_face_areaV(i,j,k) = 1.0
175 endif
176 enddo ; enddo ; enddo
177 endif
178
179 if (cs%debug) then
180 call uvchksum("Interface height used by porous barrier for layer weights", &
181 eta_u, eta_v, g%HI, haloshift=0, scalar_pair=.true.)
182 call uvchksum("Porous barrier layer-averaged weights: por_face_area[UV]", &
183 pbv%por_face_areaU, pbv%por_face_areaV, g%HI, haloshift=0, &
184 scalar_pair=.true.)
185 endif
186
187 if (cs%id_por_face_areaU > 0) call post_data(cs%id_por_face_areaU, pbv%por_face_areaU, cs%diag)
188 if (cs%id_por_face_areaV > 0) call post_data(cs%id_por_face_areaV, pbv%por_face_areaV, cs%diag)
189
190 call cpu_clock_end(id_clock_porous_barrier)
191end subroutine porous_widths_layer
192
193!> subroutine to assign porous barrier widths at the layer interfaces
194subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt)
195 ! Note: eta_bt is not currently used
196 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
197 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
198 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
199 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
200 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
201 !! thermodynamic variables.
202 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable
203 !! used to dilate the layer thicknesses
204 !! [H ~> m or kg m-2].
205 type(porous_barrier_type), intent(inout) :: pbv !< porous barrier fractional cell metrics
206 type(porous_barrier_cs), intent(in) :: cs !< Control structure for porous barrier
207
208 !local variables
209 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: eta_u ! Layer interface height at u points [Z ~> m]
210 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: eta_v ! Layer interface height at v points [Z ~> m]
211 logical, dimension(SZIB_(G),SZJB_(G)) :: do_i ! Booleans for calculation at u or v points
212 ! updated while moving up layers
213 real :: dmask ! The depth below which porous barrier is not applied [Z ~> m]
214 integer :: i, j, k, nk, is, ie, js, je, isq, ieq, jsq, jeq
215
216 if (.not.cs%initialized) call mom_error(fatal, &
217 "MOM_Porous_barrier: Module must be initialized before it is used.")
218
219 call cpu_clock_begin(id_clock_porous_barrier)
220
221 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nk = gv%ke
222 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
223
224 if (cs%answer_date < 20220806) then
225 dmask = 0.0
226 else
227 dmask = cs%mask_depth
228 endif
229
230 call calc_eta_at_uv(eta_u, eta_v, cs%eta_interp, dmask, h, tv, g, gv, us)
231
232 ! u-points
233 do j=js,je ; do i=isq,ieq
234 do_i(i,j) = .false.
235 if (g%porous_DavgU(i,j) < dmask) do_i(i,j) = .true.
236 enddo ; enddo
237
238 if (cs%answer_date < 20220806) then
239 do k=1,nk+1 ; do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
240 call calc_por_interface(g%porous_DminU(i,j), g%porous_DmaxU(i,j), g%porous_DavgU(i,j), &
241 eta_u(i,j,k), pbv%por_layer_widthU(i,j,k), do_i(i,j))
242 endif ; enddo ; enddo ; enddo
243 else
244 do k=1,nk+1 ; do j=js,je ; do i=isq,ieq
245 if (do_i(i,j)) then
246 call calc_por_interface(g%porous_DminU(i,j), g%porous_DmaxU(i,j), g%porous_DavgU(i,j), &
247 eta_u(i,j,k), pbv%por_layer_widthU(i,j,k), do_i(i,j))
248 else
249 pbv%por_layer_widthU(i,j,k) = 1.0
250 endif
251 enddo ; enddo ; enddo
252 endif
253
254 ! v-points
255 do j=jsq,jeq ; do i=is,ie
256 do_i(i,j) = .false.
257 if (g%porous_DavgV(i,j) < dmask) do_i(i,j) = .true.
258 enddo ; enddo
259
260 if (cs%answer_date < 20220806) then
261 do k=1,nk+1 ; do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
262 call calc_por_interface(g%porous_DminV(i,j), g%porous_DmaxV(i,j), g%porous_DavgV(i,j), &
263 eta_v(i,j,k), pbv%por_layer_widthV(i,j,k), do_i(i,j))
264 endif ; enddo ; enddo ; enddo
265 else
266 do k=1,nk+1 ; do j=jsq,jeq ; do i=is,ie
267 if (do_i(i,j)) then
268 call calc_por_interface(g%porous_DminV(i,j), g%porous_DmaxV(i,j), g%porous_DavgV(i,j), &
269 eta_v(i,j,k), pbv%por_layer_widthV(i,j,k), do_i(i,j))
270 else
271 pbv%por_layer_widthV(i,j,k) = 1.0
272 endif
273 enddo ; enddo ; enddo
274 endif
275
276 if (cs%debug) then
277 call uvchksum("Interface height used by porous barrier for interface weights", &
278 eta_u, eta_v, g%HI, haloshift=0, scalar_pair=.true.)
279 call uvchksum("Porous barrier weights at the layer-interface: por_layer_width[UV]", &
280 pbv%por_layer_widthU, pbv%por_layer_widthV, g%HI, &
281 haloshift=0, scalar_pair=.true.)
282 endif
283
284 if (cs%id_por_layer_widthU > 0) call post_data(cs%id_por_layer_widthU, pbv%por_layer_widthU, cs%diag)
285 if (cs%id_por_layer_widthV > 0) call post_data(cs%id_por_layer_widthV, pbv%por_layer_widthV, cs%diag)
286
287 call cpu_clock_end(id_clock_porous_barrier)
288end subroutine porous_widths_interface
289
290subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt)
291 !variables needed to call find_eta
292 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
293 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
294 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
295 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
296 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
297 !! thermodynamic variables.
298 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable
299 !! used to dilate the layer thicknesses
300 !! [H ~> m or kg m-2].
301 real, intent(in) :: dmask !< The depth shallower than which
302 !! porous barrier is not applied [Z ~> m]
303 integer, intent(in) :: interp !< eta interpolation method
304 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta_u !< Layer interface heights at u points [Z ~> m]
305 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(out) :: eta_v !< Layer interface heights at v points [Z ~> m]
306
307 ! local variables
308 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Layer interface heights [Z ~> m].
309 real :: dz_neglect ! A negligible height difference [Z ~> m]
310 integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq
311
312 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nk = gv%ke
313 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
314
315 ! currently no treatment for using optional find_eta arguments if present
316 call find_eta(h, tv, g, gv, us, eta, halo_size=1)
317
318 dz_neglect = gv%dZ_subroundoff
319
320 do k=1,nk+1
321 do j=js,je ; do i=isq,ieq ; eta_u(i,j,k) = dmask ; enddo ; enddo
322 do j=jsq,jeq ; do i=is,ie ; eta_v(i,j,k) = dmask ; enddo ; enddo
323 enddo
324
325 select case (interp)
326 case (eta_interp_max) ! The shallower interface height
327 do k=1,nk+1
328 do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
329 eta_u(i,j,k) = max(eta(i,j,k), eta(i+1,j,k))
330 endif ; enddo ; enddo
331 do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
332 eta_v(i,j,k) = max(eta(i,j,k), eta(i,j+1,k))
333 endif ; enddo ; enddo
334 enddo
335 case (eta_interp_min) ! The deeper interface height
336 do k=1,nk+1
337 do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
338 eta_u(i,j,k) = min(eta(i,j,k), eta(i+1,j,k))
339 endif ; enddo ; enddo
340 do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
341 eta_v(i,j,k) = min(eta(i,j,k), eta(i,j+1,k))
342 endif ; enddo ; enddo
343 enddo
344 case (eta_interp_arith) ! Arithmetic mean
345 do k=1,nk+1
346 do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
347 eta_u(i,j,k) = 0.5 * (eta(i,j,k) + eta(i+1,j,k))
348 endif ; enddo ; enddo
349 do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
350 eta_v(i,j,k) = 0.5 * (eta(i,j,k) + eta(i,j+1,k))
351 endif ; enddo ; enddo
352 enddo
353 case (eta_interp_harm) ! Harmonic mean
354 do k=1,nk+1
355 do j=js,je ; do i=isq,ieq ; if (g%porous_DavgU(i,j) < dmask) then
356 eta_u(i,j,k) = 2.0 * (eta(i,j,k) * eta(i+1,j,k)) / (eta(i,j,k) + eta(i+1,j,k) + dz_neglect)
357 endif ; enddo ; enddo
358 do j=jsq,jeq ; do i=is,ie ; if (g%porous_DavgV(i,j) < dmask) then
359 eta_v(i,j,k) = 2.0 * (eta(i,j,k) * eta(i,j+1,k)) / (eta(i,j,k) + eta(i,j+1,k) + dz_neglect)
360 endif ; enddo ; enddo
361 enddo
362 case default
363 call mom_error(fatal, "porous_widths::calc_eta_at_uv: "//&
364 "invalid value for eta interpolation method.")
365 end select
366end subroutine calc_eta_at_uv
367
368!> subroutine to calculate the profile fit (the three parameter fit from Adcroft 2013)
369! of the open face area fraction below a certain depth (eta_layer) in a column
370subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, A_layer, do_next)
371 real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m]
372 real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m]
373 real, intent(in) :: D_avg !< mean topographic height [Z ~> m]
374 real, intent(in) :: eta_layer !< height of interface [Z ~> m]
375 real, intent(out) :: A_layer !< frac. open face area of below eta_layer [Z ~> m]
376 logical, intent(out) :: do_next !< False if eta_layer>D_max
377
378 ! local variables
379 real :: m ! convenience constant for fit [nondim]
380 real :: zeta ! normalized vertical coordinate [nondim]
381
382 do_next = .true.
383 if (eta_layer <= d_min) then
384 a_layer = 0.0
385 elseif (eta_layer > d_max) then
386 a_layer = eta_layer - d_avg
387 do_next = .false.
388 else
389 m = (d_avg - d_min) / (d_max - d_min)
390 zeta = (eta_layer - d_min) / (d_max - d_min)
391 if (m < 0.5) then
392 a_layer = (d_max - d_min) * ((1.0 - m) * zeta**(1.0 / (1.0 - m)))
393 elseif (m == 0.5) then
394 a_layer = (d_max - d_min) * (0.5 * zeta * zeta)
395 else
396 a_layer = (d_max - d_min) * (zeta - m + m * ((1.0 - zeta)**(1.0 / m)))
397 endif
398 endif
399end subroutine calc_por_layer
400
401!> subroutine to calculate the profile fit (the three parameter fit from Adcroft 2013)
402! of the open interface fraction at a certain depth (eta_layer) in a column
403subroutine calc_por_interface(D_min, D_max, D_avg, eta_layer, w_layer, do_next)
404 real, intent(in) :: D_min !< minimum topographic height (deepest) [Z ~> m]
405 real, intent(in) :: D_max !< maximum topographic height (shallowest) [Z ~> m]
406 real, intent(in) :: D_avg !< mean topographic height [Z ~> m]
407 real, intent(in) :: eta_layer !< height of interface [Z ~> m]
408 real, intent(out) :: w_layer !< frac. open interface width at eta_layer [nondim]
409 logical, intent(out) :: do_next !< False if eta_layer>D_max
410
411 ! local variables
412 real :: m, a ! convenience constants for fit [nondim]
413 real :: zeta ! normalized vertical coordinate [nondim]
414
415 do_next = .true.
416 if (eta_layer <= d_min) then
417 w_layer = 0.0
418 elseif (eta_layer > d_max) then
419 w_layer = 1.0
420 do_next = .false.
421 else ! The following option could be refactored for stability and efficiency (with fewer divisions)
422 m = (d_avg - d_min) / (d_max - d_min)
423 a = (1.0 - m) / m
424 zeta = (eta_layer - d_min) / (d_max - d_min)
425 if (m < 0.5) then
426 w_layer = zeta**(1.0 / a)
427 ! Note that this would be safer and more efficent if it were rewritten as:
428 ! w_layer = zeta**( (D_avg - D_min) / (D_max - D_avg) )
429 elseif (m == 0.5) then
430 w_layer = zeta
431 else
432 w_layer = 1.0 - (1.0 - zeta)**a
433 endif
434 endif
435end subroutine calc_por_interface
436
437subroutine porous_barriers_init(Time, GV, US, param_file, diag, CS)
438 type(time_type), intent(in) :: time !< Current model time
439 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
440 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
441 type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse
442 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
443 type(porous_barrier_cs), intent(inout) :: cs !< Module control structure
444
445 ! local variables
446 character(len=40) :: mdl = "MOM_porous_barriers" ! This module's name.
447 character(len=20) :: interp_method ! String storing eta interpolation method
448 integer :: default_answer_date ! Global answer date
449 !> This include declares and sets the variable "version".
450# include "version_variable.h"
451
452 cs%initialized = .true.
453 cs%diag => diag
454
455 call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.false., &
456 debugging=.false.)
457 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
458 "This sets the default value for the various _ANSWER_DATE parameters.", &
459 default=99991231)
460 call get_param(param_file, mdl, "PORBAR_ANSWER_DATE", cs%answer_date, &
461 "The vintage of the porous barrier weight function calculations. Values below "//&
462 "20220806 recover the old answers in which the layer averaged weights are not "//&
463 "strictly limited by an upper-bound of 1.0 .", &
464 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
465 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)
466 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
467 call get_param(param_file, mdl, "PORBAR_MASKING_DEPTH", cs%mask_depth, &
468 "If the effective average depth at the velocity cell is shallower than this "//&
469 "number, then porous barrier is not applied at that location. "//&
470 "PORBAR_MASKING_DEPTH is assumed to be positive below the sea surface.", &
471 units="m", default=0.0, scale=us%m_to_Z)
472 ! The sign needs to be inverted to be consistent with the sign convention of Davg_[UV]
473 cs%mask_depth = -cs%mask_depth
474 call get_param(param_file, mdl, "PORBAR_ETA_INTERP", interp_method, &
475 "A string describing the method that decides how the "//&
476 "interface heights at the velocity points are calculated. "//&
477 "Valid values are:\n"//&
478 "\t MAX (the default) - maximum of the adjacent cells \n"//&
479 "\t MIN - minimum of the adjacent cells \n"//&
480 "\t ARITHMETIC - arithmetic mean of the adjacent cells \n"//&
481 "\t HARMONIC - harmonic mean of the adjacent cells \n", &
482 default=eta_interp_max_string)
483 select case (interp_method)
484 case (eta_interp_max_string) ; cs%eta_interp = eta_interp_max
485 case (eta_interp_min_string) ; cs%eta_interp = eta_interp_min
486 case (eta_interp_arith_string) ; cs%eta_interp = eta_interp_arith
487 case (eta_interp_harm_string) ; cs%eta_interp = eta_interp_harm
488 case default
489 call mom_error(fatal, "porous_barriers_init: Unrecognized setting "// &
490 "#define PORBAR_ETA_INTERP "//trim(interp_method)//" found in input file.")
491 end select
492
493 cs%id_por_layer_widthU = register_diag_field('ocean_model', 'por_layer_widthU', diag%axesCui, time, &
494 'Porous barrier open width fraction (at the layer interfaces) of the u-faces', 'nondim')
495 cs%id_por_layer_widthV = register_diag_field('ocean_model', 'por_layer_widthV', diag%axesCvi, time, &
496 'Porous barrier open width fraction (at the layer interfaces) of the v-faces', 'nondim')
497 cs%id_por_face_areaU = register_diag_field('ocean_model', 'por_face_areaU', diag%axesCuL, time, &
498 'Porous barrier open area fraction (layer averaged) of U-faces', 'nondim')
499 cs%id_por_face_areaV = register_diag_field('ocean_model', 'por_face_areaV', diag%axesCvL, time, &
500 'Porous barrier open area fraction (layer averaged) of V-faces', 'nondim')
501
502 id_clock_porous_barrier = cpu_clock_id('(Ocean porous barrier)', grain=clock_module)
503end subroutine
504
505end module mom_porous_barriers