MOM_interface_heights.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!> Functions for calculating interface heights, including free surface height.
7
9use mom_debugging, only : hchksum
10use mom_error_handler, only : mom_error, fatal
11use mom_eos, only : calculate_density, average_specific_vol, eos_type, eos_domain
12use mom_file_parser, only : log_version
13use mom_grid, only : ocean_grid_type
17
18implicit none ; private
19
20#include <MOM_memory.h>
21
22public find_eta, find_dz_for_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
26
27!> Calculates the heights of the free surface or all interfaces from layer thicknesses.
28interface find_eta
29 module procedure find_eta_2d, find_eta_3d
30end interface find_eta
31
32!> Calculates layer thickness in thickness units from geometric distance between the
33!! interfaces around that layer in height units.
34interface dz_to_thickness
36end interface dz_to_thickness
37
38!> Converts layer thickness in thickness units into the vertical distance between the
39!! interfaces around a layer in height units.
40interface thickness_to_dz
42end interface thickness_to_dz
43
44contains
45
46!> Calculates the change in height across layers, using the appropriate form for
47!! consistency with the calculation of the pressure gradient forces.
48subroutine find_dz_for_eta(h, tv, G, GV, US, dz_lay, halo_size)
49 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
50 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
51 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
52 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
53 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
54 !! thermodynamic variables.
55 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: dz_lay !< Height change across layers [Z ~> m]
56 integer, optional, intent(in) :: halo_size !< width of halo points on
57 !! which to calculate eta.
58
59 ! Local variables
60 real :: p(szi_(g),szj_(g),szk_(gv)+1) ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
61 real :: dz_geo(szi_(g),szj_(g)) ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2]
62 real :: spv_lay_conv(szk_(gv)) ! The prescribed layer specific volume times a conversion factor from
63 ! the units of thickness to layer mass [Z H-1 ~> nondim or m3 kg-1]
64 real :: i_gearth ! The inverse of the gravitational acceleration times the
65 ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
66 integer :: i, j, k, isv, iev, jsv, jev, nz, halo
67
68 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
69
70 isv = g%isc-halo ; iev = g%iec+halo ; jsv = g%jsc-halo ; jev = g%jec+halo
71 nz = gv%ke
72
73 if ((isv<g%isd) .or. (iev>g%ied) .or. (jsv<g%jsd) .or. (jev>g%jed)) &
74 call mom_error(fatal,"find_dz_for_eta called with an overly large halo_size.")
75
76 if (gv%Boussinesq) then
77 do k=1,nz ; do j=jsv,jev ; do i=isv,iev
78 dz_lay(i,j,k) = h(i,j,k)*gv%H_to_Z
79 enddo ; enddo ; enddo
80 elseif (associated(tv%eqn_of_state)) then
81 i_gearth = 1.0 / gv%g_Earth
82 !$OMP parallel do default(shared)
83 do j=jsv,jev
84 if (associated(tv%p_surf)) then
85 do i=isv,iev ; p(i,j,1) = tv%p_surf(i,j) ; enddo
86 else
87 do i=isv,iev ; p(i,j,1) = 0.0 ; enddo
88 endif
89 do k=1,nz ; do i=isv,iev
90 p(i,j,k+1) = p(i,j,k) + gv%g_Earth*gv%H_to_RZ*h(i,j,k)
91 enddo ; enddo
92 enddo
93 !$OMP parallel do default(shared) private(dz_geo)
94 do k=1,nz
95 call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
96 0.0, g%HI, tv%eqn_of_state, us, dz_geo, halo_size=halo)
97 do j=jsv,jev ; do i=isv,iev
98 dz_lay(i,j,k) = i_gearth * dz_geo(i,j)
99 enddo ; enddo
100 enddo
101 else ! non-Boussinesq but with no equation of state
102 do k=1,nz ; do j=jsv,jev ; do i=isv,iev
103 dz_lay(i,j,k) = gv%H_to_RZ*h(i,j,k) / gv%Rlay(k)
104 enddo ; enddo ; enddo
105 ! This would be faster but could change answers.
106 ! do k=1,nz ; SpV_lay_conv(k) = GV%H_to_RZ / GV%Rlay(k) ; enddo
107 ! do k=1,nz ; do j=jsv,jev ; do i=isv,iev
108 ! dz_lay(i,j,K) = h(i,j,k) * SpV_lay_conv(k)
109 ! enddo ; enddo ; enddo
110 endif
111
112 ! To find eta, do the following:
113 ! do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) ; enddo ; enddo
114 ! do k=nz,1,-1 ; do j=jsv,jev ; do i=isv,iev
115 ! eta(i,j,K) = eta(i,j,K+1) + dz_lay(i,j,K)
116 ! enddo ; enddo ; enddo
117
118end subroutine find_dz_for_eta
119
120!> Calculates the heights of all interfaces between layers, using the appropriate
121!! form for consistency with the calculation of the pressure gradient forces.
122!! Additionally, these height may be dilated for consistency with the
123!! corresponding time-average quantity from the barotropic calculation.
124subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
125 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
126 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
127 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
128 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
129 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
130 !! thermodynamic variables.
131 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta !< layer interface heights [Z ~> m]
132 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable
133 !! that gives the "correct" free surface height (Boussinesq) or total water
134 !! column mass per unit area (non-Boussinesq). This is used to dilate the layer
135 !! thicknesses when calculating interface heights [H ~> m or kg m-2].
136 !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height.
137 integer, optional, intent(in) :: halo_size !< width of halo points on
138 !! which to calculate eta.
139 real, optional, intent(in) :: dZref !< The difference in the
140 !! reference height between G%bathyT and eta [Z ~> m]. The default is 0.
141
142 ! Local variables
143 real :: dz_lay(SZI_(G),SZJ_(G),SZK_(GV)) ! The change in height across a layer [Z ~> m]
144 real :: dilate(SZI_(G)) ! A non-dimensional dilation factor [nondim]
145 real :: htot(SZI_(G)) ! total thickness [H ~> m or kg m-2]
146 real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m].
147 ! dZ_ref is 0 unless the optional argument dZref is present.
148 integer :: i, j, k, isv, iev, jsv, jev, nz, halo
149
150 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
151
152 isv = g%isc-halo ; iev = g%iec+halo ; jsv = g%jsc-halo ; jev = g%jec+halo
153 nz = gv%ke
154
155 if ((isv<g%isd) .or. (iev>g%ied) .or. (jsv<g%jsd) .or. (jev>g%jed)) &
156 call mom_error(fatal,"find_eta called with an overly large halo_size.")
157
158 dz_ref = 0.0 ; if (present(dzref)) dz_ref = dzref
159
160 if (gv%Boussinesq) then
161 !$OMP parallel default(shared) private(dilate,htot)
162 !$OMP do
163 do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -(g%bathyT(i,j) + dz_ref) ; enddo ; enddo
164 !$OMP do
165 do j=jsv,jev ; do k=nz,1,-1 ; do i=isv,iev
166 eta(i,j,k) = eta(i,j,k+1) + h(i,j,k)*gv%H_to_Z
167 enddo ; enddo ; enddo
168 if (present(eta_bt)) then
169 ! Dilate the water column to agree with the free surface height
170 ! that is used for the dynamics.
171 !$OMP do
172 do j=jsv,jev !$OMP parallel do default(shared)
173
174 do i=isv,iev
175 dilate(i) = (eta_bt(i,j)*gv%H_to_Z + g%bathyT(i,j)) / &
176 (eta(i,j,1) + (g%bathyT(i,j) + dz_ref))
177 enddo
178 do k=1,nz ; do i=isv,iev
179 eta(i,j,k) = dilate(i) * (eta(i,j,k) + (g%bathyT(i,j) + dz_ref)) - &
180 (g%bathyT(i,j) + dz_ref)
181 enddo ; enddo
182 enddo
183 endif
184 !$OMP end parallel
185 else
186 call find_dz_for_eta(h, tv, g, gv, us, dz_lay, halo_size)
187 !$OMP parallel default(shared) private(dilate,htot)
188 !$OMP do
189 do j=jsv,jev
190 do i=isv,iev ; eta(i,j,nz+1) = -(g%bathyT(i,j) + dz_ref) ; enddo
191 do k=nz,1,-1 ; do i=isv,iev
192 eta(i,j,k) = eta(i,j,k+1) + dz_lay(i,j,k)
193 enddo ; enddo
194 enddo
195
196 if (present(eta_bt)) then
197 ! Dilate the water column to agree with the free surface height
198 ! from the time-averaged barotropic solution.
199 !$OMP do
200 do j=jsv,jev
201 do i=isv,iev ; htot(i) = gv%H_subroundoff ; enddo
202 do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
203 do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ; enddo
204 do k=1,nz ; do i=isv,iev
205 eta(i,j,k) = dilate(i) * (eta(i,j,k) + (g%bathyT(i,j) + dz_ref)) - &
206 (g%bathyT(i,j) + dz_ref)
207 enddo ; enddo
208 enddo
209 endif
210 !$OMP end parallel
211 endif
212
213end subroutine find_eta_3d
214
215!> Calculates the free surface height, using the appropriate form for consistency
216!! with the calculation of the pressure gradient forces. Additionally, the sea
217!! surface height may be adjusted for consistency with the corresponding
218!! time-average quantity from the barotropic calculation.
219subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
220 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
221 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
222 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
223 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
224 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
225 !! thermodynamic variables.
226 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta !< free surface height relative to
227 !! mean sea level (z=0) often [Z ~> m].
228 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic
229 !! variable that gives the "correct" free surface height (Boussinesq) or total
230 !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2].
231 !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height.
232 integer, optional, intent(in) :: halo_size !< width of halo points on
233 !! which to calculate eta.
234 real, optional, intent(in) :: dZref !< The difference in the
235 !! reference height between G%bathyT and eta [Z ~> m]. The default is 0.
236
237 ! Local variables
238 real :: dz_lay(SZI_(G),SZJ_(G),SZK_(GV)) ! The change in height across a layer [Z ~> m]
239 real :: htot(SZI_(G)) ! The sum of all layers' thicknesses [H ~> m or kg m-2].
240 real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m].
241 ! dZ_ref is 0 unless the optional argument dZref is present.
242 integer :: i, j, k, is, ie, js, je, nz, halo
243
244 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
245 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
246 nz = gv%ke
247
248 dz_ref = 0.0 ; if (present(dzref)) dz_ref = dzref
249
250 if (gv%Boussinesq) then
251 if (present(eta_bt)) then
252 !$OMP parallel do default(shared)
253 do j=js,je ; do i=is,ie
254 eta(i,j) = gv%H_to_Z*eta_bt(i,j) - dz_ref
255 enddo ; enddo
256 else
257 !$OMP parallel do default(shared)
258 do j=js,je
259 do i=is,ie ; eta(i,j) = -(g%bathyT(i,j) + dz_ref) ; enddo
260 do k=1,nz ; do i=is,ie
261 eta(i,j) = eta(i,j) + h(i,j,k)*gv%H_to_Z
262 enddo ; enddo
263 enddo
264 endif
265 else
266 call find_dz_for_eta(h, tv, g, gv, us, dz_lay, halo_size)
267 !$OMP parallel default(shared) private(htot)
268 !$OMP do
269 do j=js,je
270 do i=is,ie ; eta(i,j) = -(g%bathyT(i,j) + dz_ref) ; enddo
271 do k=1,nz ; do i=is,ie
272 eta(i,j) = eta(i,j) + dz_lay(i,j,k)
273 enddo ; enddo
274 enddo
275 if (present(eta_bt)) then
276 ! Dilate the water column to agree with the time-averaged column
277 ! mass from the barotropic solution.
278 !$OMP do
279 do j=js,je
280 do i=is,ie ; htot(i) = gv%H_subroundoff ; enddo
281 do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
282 do i=is,ie
283 eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + (g%bathyT(i,j) + dz_ref)) - &
284 (g%bathyT(i,j) + dz_ref)
285 enddo
286 enddo
287 endif
288 !$OMP end parallel
289 endif
290
291end subroutine find_eta_2d
292
293
294!> Calculate derived thermodynamic quantities for re-use later.
295subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
296 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
297 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
298 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
299 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various
300 !! thermodynamic variables, some of
301 !! which will be set here.
302 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
303 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
304 integer, optional, intent(in) :: halo !< Width of halo within which to
305 !! calculate thicknesses
306 logical, optional, intent(in) :: debug !< If present and true, write debugging checksums
307 ! Local variables
308 real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa]
309 real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa]
310 real, dimension(SZK_(GV)) :: spv_lay ! The specific volume of each layer when no equation of
311 ! state is used [R-1 ~> m3 kg-1]
312 logical :: do_debug ! If true, write checksums for debugging.
313 integer :: i, j, k, is, ie, js, je, halos, nz
314
315 do_debug = .false. ; if (present(debug)) do_debug = debug
316 halos = 0 ; if (present(halo)) halos = max(0,halo)
317 is = g%isc-halos ; ie = g%iec+halos ; js = g%jsc-halos ; je = g%jec+halos ; nz = gv%ke
318
319 if (allocated(tv%Spv_avg) .and. associated(tv%eqn_of_state)) then
320 if (associated(tv%p_surf)) then
321 do j=js,je ; do i=is,ie ; p_t(i,j) = tv%p_surf(i,j) ; enddo ; enddo
322 else
323 do j=js,je ; do i=is,ie ; p_t(i,j) = 0.0 ; enddo ; enddo
324 endif
325 do k=1,nz
326 do j=js,je ; do i=is,ie
327 dp(i,j) = gv%g_Earth*gv%H_to_RZ*h(i,j,k)
328 enddo ; enddo
329 call avg_specific_vol(tv%T(:,:,k), tv%S(:,:,k), p_t, dp, g%HI, tv%eqn_of_state, tv%SpV_avg(:,:,k), halo)
330 if (k<nz) then ; do j=js,je ; do i=is,ie
331 p_t(i,j) = p_t(i,j) + dp(i,j)
332 enddo ; enddo ; endif
333 enddo
334 tv%valid_SpV_halo = halos
335
336 if (do_debug) then
337 call hchksum(h, "derived_thermo h", g%HI, haloshift=halos, unscale=gv%H_to_MKS)
338 if (associated(tv%p_surf)) call hchksum(tv%p_surf, "derived_thermo p_surf", g%HI, &
339 haloshift=halos, unscale=us%RL2_T2_to_Pa)
340 call hchksum(tv%T, "derived_thermo T", g%HI, haloshift=halos, unscale=us%C_to_degC)
341 call hchksum(tv%S, "derived_thermo S", g%HI, haloshift=halos, unscale=us%S_to_ppt)
342 endif
343 elseif (allocated(tv%Spv_avg)) then
344 do k=1,nz ; spv_lay(k) = 1.0 / gv%Rlay(k) ; enddo
345 do k=1,nz ; do j=js,je ; do i=is,ie
346 tv%SpV_avg(i,j,k) = spv_lay(k)
347 enddo ; enddo ; enddo
348 tv%valid_SpV_halo = halos
349 endif
350
351end subroutine calc_derived_thermo
352
353
354!> Determine the column average specific volumes.
355subroutine find_col_avg_spv(h, SpV_avg, tv, G, GV, US, halo_size)
356 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
357 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
358 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
359 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
360 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
361 real, dimension(SZI_(G),SZJ_(G)), &
362 intent(inout) :: spv_avg !< Column average specific volume [R-1 ~> m3 kg-1]
363 ! SpV_avg is intent inout to retain excess halo values.
364 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
365 !! thermodynamic fields.
366 integer, optional, intent(in) :: halo_size !< width of halo points on which to work
367
368 ! Local variables
369 real :: h_tot(szi_(g)) ! Sum of the layer thicknesses [H ~> m or kg m-2]
370 real :: spv_x_h_tot(szi_(g)) ! Vertical sum of the layer average specific volume times
371 ! the layer thicknesses [H R-1 ~> m4 kg-1 or m]
372 real :: i_rho ! The inverse of the Boussiensq reference density [R-1 ~> m3 kg-1]
373 real :: spv_lay(szk_(gv)) ! The inverse of the layer target potential densities [R-1 ~> m3 kg-1]
374 character(len=128) :: mesg ! A string for error messages
375 integer :: i, j, k, is, ie, js, je, nz, halo
376
377 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
378
379 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
380 nz = gv%ke
381
382 if (gv%Boussinesq) then
383 i_rho = 1.0 / gv%Rho0
384 do j=js,je ; do i=is,ie
385 spv_avg(i,j) = i_rho
386 enddo ; enddo
387 elseif (.not.allocated(tv%SpV_avg)) then
388 do k=1,nz ; spv_lay(k) = 1.0 / gv%Rlay(k) ; enddo
389 do j=js,je
390 do i=is,ie ; spv_x_h_tot(i) = 0.0 ; h_tot(i) = 0.0 ; enddo
391 do k=1,nz ; do i=is,ie
392 h_tot(i) = h_tot(i) + max(h(i,j,k), gv%H_subroundoff)
393 spv_x_h_tot(i) = spv_x_h_tot(i) + spv_lay(k)*max(h(i,j,k), gv%H_subroundoff)
394 enddo ; enddo
395 do i=is,ie ; spv_avg(i,j) = spv_x_h_tot(i) / h_tot(i) ; enddo
396 enddo
397 else
398 ! Check that SpV_avg has been set.
399 if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then
400 if (tv%valid_SpV_halo < 0) then
401 mesg = "invalid values of SpV_avg."
402 else
403 write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
404 tv%valid_SpV_halo, halo
405 endif
406 call mom_error(fatal, "find_col_avg_SpV called in fully non-Boussinesq mode with "//trim(mesg))
407 endif
408
409 do j=js,je
410 do i=is,ie ; spv_x_h_tot(i) = 0.0 ; h_tot(i) = 0.0 ; enddo
411 do k=1,nz ; do i=is,ie
412 h_tot(i) = h_tot(i) + max(h(i,j,k), gv%H_subroundoff)
413 spv_x_h_tot(i) = spv_x_h_tot(i) + tv%SpV_avg(i,j,k)*max(h(i,j,k), gv%H_subroundoff)
414 enddo ; enddo
415 do i=is,ie ; spv_avg(i,j) = spv_x_h_tot(i) / h_tot(i) ; enddo
416 enddo
417 endif
418
419end subroutine find_col_avg_spv
420
421!> Calculate the integrated mass of the water column.
422subroutine find_col_mass(h, tv, G, GV, US, mass, p_bot, p_surf)
423 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
424 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
425 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
426 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
427 !! thermodynamic variables.
428 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
429 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass !< Integrated mass of the water column
430 !! [R Z ~> kg m-2]
431 real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: p_bot !< Bottom pressure = g * mass + psurf
432 !! [R L2 T-2 ~> Pa]
433 real, dimension(:,:), optional, pointer :: p_surf !< A pointer to surface pressure
434 !! [R L2 T-2 ~> Pa]
435
436 ! Local variables
437 real :: i_gearth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1]
438 real, dimension(SZI_(G),SZJ_(G)) :: &
439 z_top, & ! Height of the top of a layer [Z ~> m].
440 z_bot, & ! Height of the bottom of a layer [Z ~> m].
441 dp ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa].
442 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
443
444 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
445 isq = g%iscB ; ieq = g%iecB ; jsq = g%jscB ; jeq = g%jecB
446 nz = gv%ke
447
448 do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
449 if (gv%Boussinesq) then
450 if (associated(tv%eqn_of_state)) then
451 i_gearth = 1.0 / gv%g_Earth
452 do j=jsq,jeq+1 ; do i=isq,ieq+1 ; z_bot(i,j) = 0.0 ; enddo ; enddo
453 do k=1,nz
454 ! NOTE: int_density_z expects z_top and z_bot values from [ij]sq to [ij]eq+1
455 do j=jsq,jeq+1 ; do i=isq,ieq+1
456 z_top(i,j) = z_bot(i,j)
457 z_bot(i,j) = z_top(i,j) - gv%H_to_Z * h(i,j,k)
458 enddo ; enddo
459 call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, gv%Rho0, gv%g_Earth, &
460 g%HI, tv%eqn_of_state, us, dp)
461 do j=js,je ; do i=is,ie
462 mass(i,j) = mass(i,j) + dp(i,j) * i_gearth
463 enddo ; enddo
464 enddo
465 else
466 do k=1,nz ; do j=js,je ; do i=is,ie
467 mass(i,j) = mass(i,j) + (gv%H_to_Z * gv%Rlay(k)) * h(i,j,k)
468 enddo ; enddo ; enddo
469 endif
470 else
471 do k=1,nz ; do j=js,je ; do i=is,ie
472 mass(i,j) = mass(i,j) + gv%H_to_RZ * h(i,j,k)
473 enddo ; enddo ; enddo
474 endif
475
476 if (present(p_bot)) then
477 do j=js,je ; do i=is,ie
478 p_bot(i,j) = gv%g_Earth * mass(i,j)
479 enddo ; enddo
480 if (present(p_surf) .and. associated(p_surf)) then ; do j=js,je ; do i=is,ie
481 p_bot(i,j) = p_bot(i,j) + p_surf(i,j)
482 enddo ; enddo ; endif
483 endif
484
485end subroutine find_col_mass
486
487!> Determine the in situ density averaged over a specified distance from the bottom,
488!! calculating it as the inverse of the mass-weighted average specific volume.
489subroutine find_rho_bottom(G, GV, US, tv, h, dz, pres_int, dz_avg, j, Rho_bot, h_bot, k_bot)
490 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
491 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
492 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
493 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
494 !! thermodynamic fields.
495 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
496 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
497 real, dimension(SZI_(G),SZK_(GV)), &
498 intent(in) :: dz !< Height change across layers [Z ~> m]
499 real, dimension(SZI_(G),SZK_(GV)+1), &
500 intent(in) :: pres_int !< Pressure at each interface [R L2 T-2 ~> Pa]
501 real, dimension(SZI_(G)), intent(in) :: dz_avg !< The vertical distance over which to average [Z ~> m]
502 integer, intent(in) :: j !< j-index of row to work on
503 real, dimension(SZI_(G)), intent(out) :: rho_bot !< Near-bottom density [R ~> kg m-3].
504 real, dimension(SZI_(G)), intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
505 integer, dimension(SZI_(G)), intent(out) :: k_bot !< Bottom boundary layer top layer index
506
507 ! Local variables
508 real :: hb(szi_(g)) ! Running sum of the thickness in the bottom boundary layer [H ~> m or kg m-2]
509 real :: spv_h_bot(szi_(g)) ! Running sum of the specific volume times thickness in the bottom
510 ! boundary layer [H R-1 ~> m4 kg-1 or m]
511 real :: dz_bbl_rem(szi_(g)) ! Vertical extent of the boundary layer that has yet to be accounted
512 ! for [Z ~> m]
513 real :: h_bbl_frac(szi_(g)) ! Thickness of the fractional layer that makes up the top of the
514 ! boundary layer [H ~> m or kg m-2]
515 real :: t_bbl(szi_(g)) ! Temperature of the fractional layer that makes up the top of the
516 ! boundary layer [C ~> degC]
517 real :: s_bbl(szi_(g)) ! Salinity of the fractional layer that makes up the top of the
518 ! boundary layer [S ~> ppt]
519 real :: p_bbl(szi_(g)) ! Pressure the top of the boundary layer [R L2 T-2 ~> Pa]
520 real :: dp(szi_(g)) ! Pressure change across the fractional layer that makes up the top
521 ! of the boundary layer [R L2 T-2 ~> Pa]
522 real :: spv_bbl(szi_(g)) ! In situ specific volume of the fractional layer that makes up the
523 ! top of the boundary layer [R-1 ~> m3 kg-1]
524 real :: frac_in ! The fraction of a layer that is within the bottom boundary layer [nondim]
525 logical :: do_i(szi_(g)), do_any
526 logical :: use_eos
527 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
528 integer :: i, k, is, ie, nz
529
530 is = g%isc ; ie = g%iec ; nz = gv%ke
531
532 use_eos = associated(tv%T) .and. associated(tv%S) .and. associated(tv%eqn_of_state)
533
534 if (gv%Boussinesq .or. gv%semi_Boussinesq .or. .not.allocated(tv%SpV_avg)) then
535 do i=is,ie
536 rho_bot(i) = gv%Rho0
537 enddo
538
539 ! Obtain bottom boundary layer thickness and index of top layer
540 do i=is,ie
541 hb(i) = 0.0 ; h_bot(i) = 0.0 ; k_bot(i) = nz
542 dz_bbl_rem(i) = g%mask2dT(i,j) * max(0.0, dz_avg(i))
543 do_i(i) = .true.
544 if (g%mask2dT(i,j) <= 0.0) then
545 h_bbl_frac(i) = 0.0
546 do_i(i) = .false.
547 endif
548 enddo
549
550 do k=nz,1,-1
551 do_any = .false.
552 do i=is,ie ; if (do_i(i)) then
553 if (dz(i,k) < dz_bbl_rem(i)) then
554 ! This layer is fully within the averaging depth.
555 dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k)
556 hb(i) = hb(i) + h(i,j,k)
557 k_bot(i) = k
558 do_any = .true.
559 else
560 if (dz(i,k) > 0.0) then
561 frac_in = dz_bbl_rem(i) / dz(i,k)
562 if (frac_in >= 0.5) k_bot(i) = k ! update bbl top index if >= 50% of layer
563 else
564 frac_in = 0.0
565 endif
566 h_bbl_frac(i) = frac_in * h(i,j,k)
567 dz_bbl_rem(i) = 0.0
568 do_i(i) = .false.
569 endif
570 endif ; enddo
571 if (.not.do_any) exit
572 enddo
573 do i=is,ie ; if (do_i(i)) then
574 ! The nominal bottom boundary layer is thicker than the water column, but layer 1 is
575 ! already included in the averages. These values are set so that the call to find
576 ! the layer-average specific volume will behave sensibly.
577 h_bbl_frac(i) = 0.0
578 endif ; enddo
579
580 do i=is,ie
581 if (hb(i) + h_bbl_frac(i) < gv%H_subroundoff) h_bbl_frac(i) = gv%H_subroundoff
582 h_bot(i) = hb(i) + h_bbl_frac(i)
583 enddo
584
585 else
586 ! Check that SpV_avg has been set.
587 if (tv%valid_SpV_halo < 0) call mom_error(fatal, &
588 "find_rho_bottom called in fully non-Boussinesq mode with invalid values of SpV_avg.")
589
590 ! Set the bottom density to the inverse of the in situ specific volume averaged over the
591 ! specified distance, with care taken to avoid having compressibility lead to an imprint
592 ! of the layer thicknesses on this density.
593 do i=is,ie
594 hb(i) = 0.0 ; spv_h_bot(i) = 0.0 ; h_bot(i) = 0.0 ; k_bot(i) = nz
595 dz_bbl_rem(i) = g%mask2dT(i,j) * max(0.0, dz_avg(i))
596 do_i(i) = .true.
597 if (g%mask2dT(i,j) <= 0.0) then
598 ! Set acceptable values for calling the equation of state over land.
599 t_bbl(i) = 0.0 ; s_bbl(i) = 0.0 ; dp(i) = 0.0 ; p_bbl(i) = 0.0
600 spv_bbl(i) = 1.0 ! This value is arbitrary, provided it is non-zero.
601 h_bbl_frac(i) = 0.0
602 do_i(i) = .false.
603 endif
604 enddo
605
606 do k=nz,1,-1
607 do_any = .false.
608 do i=is,ie ; if (do_i(i)) then
609 if (dz(i,k) < dz_bbl_rem(i)) then
610 ! This layer is fully within the averaging depth.
611 spv_h_bot(i) = spv_h_bot(i) + h(i,j,k) * tv%SpV_avg(i,j,k)
612 dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k)
613 hb(i) = hb(i) + h(i,j,k)
614 k_bot(i) = k
615 do_any = .true.
616 else
617 if (dz(i,k) > 0.0) then
618 frac_in = dz_bbl_rem(i) / dz(i,k)
619 if (frac_in >= 0.5) k_bot(i) = k ! update bbl top index if >= 50% of layer
620 else
621 frac_in = 0.0
622 endif
623 if (use_eos) then
624 ! Store the properties of this layer to determine the average
625 ! specific volume of the portion that is within the BBL.
626 t_bbl(i) = tv%T(i,j,k) ; s_bbl(i) = tv%S(i,j,k)
627 dp(i) = frac_in * (gv%g_Earth*gv%H_to_RZ * h(i,j,k))
628 p_bbl(i) = pres_int(i,k) + (1.0-frac_in) * (gv%g_Earth*gv%H_to_RZ * h(i,j,k))
629 else
630 spv_bbl(i) = tv%SpV_avg(i,j,k)
631 endif
632 h_bbl_frac(i) = frac_in * h(i,j,k)
633 dz_bbl_rem(i) = 0.0
634 do_i(i) = .false.
635 endif
636 endif ; enddo
637 if (.not.do_any) exit
638 enddo
639 do i=is,ie ; if (do_i(i)) then
640 ! The nominal bottom boundary layer is thicker than the water column, but layer 1 is
641 ! already included in the averages. These values are set so that the call to find
642 ! the layer-average specific volume will behave sensibly.
643 if (use_eos) then
644 t_bbl(i) = tv%T(i,j,1) ; s_bbl(i) = tv%S(i,j,1)
645 dp(i) = 0.0
646 p_bbl(i) = pres_int(i,1)
647 else
648 spv_bbl(i) = tv%SpV_avg(i,j,1)
649 endif
650 h_bbl_frac(i) = 0.0
651 endif ; enddo
652
653 if (use_eos) then
654 ! Find the average specific volume of the fractional layer atop the BBL.
655 eosdom(:) = eos_domain(g%HI)
656 call average_specific_vol(t_bbl, s_bbl, p_bbl, dp, spv_bbl, tv%eqn_of_state, eosdom)
657 endif
658
659 do i=is,ie
660 if (hb(i) + h_bbl_frac(i) < gv%H_subroundoff) h_bbl_frac(i) = gv%H_subroundoff
661 rho_bot(i) = g%mask2dT(i,j) * (hb(i) + h_bbl_frac(i)) / (spv_h_bot(i) + h_bbl_frac(i)*spv_bbl(i))
662 h_bot(i) = hb(i) + h_bbl_frac(i)
663 enddo
664 endif
665
666end subroutine find_rho_bottom
667
668
669!> Converts thickness from geometric height units to thickness units, perhaps via an
670!! inversion of the integral of the density in pressure using variables stored in
671!! the thermo_var_ptrs type when in non-Boussinesq mode.
672subroutine dz_to_thickness_tv(dz, tv, h, G, GV, US, halo_size)
673 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
674 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
675 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
676 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
677 intent(in) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
678 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
679 !! thermodynamic variables
680 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
681 intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2].
682 !! This is essentially intent out, but declared as intent
683 !! inout to preserve any initialized values in halo points.
684 integer, optional, intent(in) :: halo_size !< Width of halo within which to
685 !! calculate thicknesses
686 ! Local variables
687 integer :: i, j, k, is, ie, js, je, halo, nz
688
689 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
690 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
691
692 if (gv%Boussinesq) then
693 do k=1,nz ; do j=js,je ; do i=is,ie
694 h(i,j,k) = gv%Z_to_H * dz(i,j,k)
695 enddo ; enddo ; enddo
696 else
697 if (associated(tv%eqn_of_state)) then
698 if (associated(tv%p_surf)) then
699 call dz_to_thickness_eos(dz, tv%T, tv%S, tv%eqn_of_state, h, g, gv, us, halo, tv%p_surf)
700 else
701 call dz_to_thickness_eos(dz, tv%T, tv%S, tv%eqn_of_state, h, g, gv, us, halo)
702 endif
703 else
704 do k=1,nz ; do j=js,je ; do i=is,ie
705 h(i,j,k) = (gv%RZ_to_H * gv%Rlay(k)) * dz(i,j,k)
706 enddo ; enddo ; enddo
707 endif
708 endif
709
710end subroutine dz_to_thickness_tv
711
712!> Converts thickness from geometric height units to thickness units, working via an
713!! inversion of the integral of the density in pressure when in non-Boussinesq mode.
714subroutine dz_to_thickness_eos(dz, Temp, Saln, EoS, h, G, GV, US, halo_size, p_surf)
715 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
716 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
717 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
718 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
719 intent(in) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
720 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
721 intent(in) :: Temp !< Input layer temperatures [C ~> degC]
722 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
723 intent(in) :: Saln !< Input layer salinities [S ~> ppt]
724 type(eos_type), intent(in) :: EoS !< Equation of state structure
725 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
726 intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2].
727 !! This is essentially intent out, but declared as intent
728 !! inout to preserve any initialized values in halo points.
729 integer, optional, intent(in) :: halo_size !< Width of halo within which to
730 !! calculate thicknesses
731 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: p_surf !< Surface pressures [R L2 T-2 ~> Pa]
732 ! Local variables
733 real, dimension(SZI_(G),SZJ_(G)) :: &
734 p_top, p_bot ! Pressure at the interfaces above and below a layer [R L2 T-2 ~> Pa]
735 real :: dp(SZI_(G),SZJ_(G)) ! Pressure change across a layer [R L2 T-2 ~> Pa]
736 real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2]
737 real :: rho(SZI_(G)) ! The in situ density [R ~> kg m-3]
738 real :: dp_adj ! The amount by which to change the bottom pressure in an
739 ! iteration [R L2 T-2 ~> Pa]
740 real :: I_gEarth ! Unit conversion factors divided by the gravitational
741 ! acceleration [H T2 R-1 L-2 ~> s2 m2 kg-1 or s2 m-1]
742 logical :: do_more(SZI_(G),SZJ_(G)) ! If true, additional iterations would be beneficial.
743 logical :: do_any ! True if there are points in this layer that need more itertions.
744 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
745 integer :: i, j, k, is, ie, js, je, halo, nz
746 integer :: itt, max_itt
747
748 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
749 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
750 max_itt = 10
751
752 if (gv%Boussinesq) then
753 do k=1,nz ; do j=js,je ; do i=is,ie
754 h(i,j,k) = gv%Z_to_H * dz(i,j,k)
755 enddo ; enddo ; enddo
756 else
757 i_gearth = gv%RZ_to_H / gv%g_Earth
758
759 if (present(p_surf)) then
760 do j=js,je ; do i=is,ie
761 p_bot(i,j) = 0.0 ; p_top(i,j) = p_surf(i,j)
762 enddo ; enddo
763 else
764 do j=js,je ; do i=is,ie
765 p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
766 enddo ; enddo
767 endif
768 eosdom(:) = eos_domain(g%HI)
769
770 ! The iterative approach here is inherited from very old code that was in the
771 ! MOM_state_initialization module. It does converge, but it is very inefficient and
772 ! should be revised, although doing so would change answers in non-Boussinesq mode.
773 do k=1,nz
774 do j=js,je
775 do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
776 call calculate_density(temp(:,j,k), saln(:,j,k), p_top(:,j), rho, &
777 eos, eosdom)
778 ! The following two expressions are mathematically equivalent.
779 if (gv%semi_Boussinesq) then
780 do i=is,ie
781 p_bot(i,j) = p_top(i,j) + (gv%g_Earth*gv%H_to_Z) * ((gv%Z_to_H*dz(i,j,k)) * rho(i))
782 dp(i,j) = (gv%g_Earth*gv%H_to_Z) * ((gv%Z_to_H*dz(i,j,k)) * rho(i))
783 enddo
784 else
785 do i=is,ie
786 p_bot(i,j) = p_top(i,j) + rho(i) * (gv%g_Earth * dz(i,j,k))
787 dp(i,j) = rho(i) * (gv%g_Earth * dz(i,j,k))
788 enddo
789 endif
790 enddo
791
792 do_more(:,:) = .true.
793 do itt=1,max_itt
794 do_any = .false.
795 call int_specific_vol_dp(temp(:,:,k), saln(:,:,k), p_top, p_bot, 0.0, g%HI, eos, us, dz_geo)
796 if (itt < max_itt) then ; do j=js,je
797 call calculate_density(temp(:,j,k), saln(:,j,k), p_bot(:,j), rho, eos, eosdom)
798 ! Use Newton's method to correct the bottom value.
799 ! The hydrostatic equation is sufficiently linear that no bounds-checking is needed.
800 if (gv%semi_Boussinesq) then
801 do i=is,ie
802 dp_adj = rho(i) * ((gv%g_Earth*gv%H_to_Z)*(gv%Z_to_H*dz(i,j,k)) - dz_geo(i,j))
803 p_bot(i,j) = p_bot(i,j) + dp_adj
804 dp(i,j) = dp(i,j) + dp_adj
805 enddo
806 do_any = .true. ! To avoid changing answers, always use the maximum number of itertions.
807 else
808 do i=is,ie ; if (do_more(i,j)) then
809 dp_adj = rho(i) * (gv%g_Earth*dz(i,j,k) - dz_geo(i,j))
810 p_bot(i,j) = p_bot(i,j) + dp_adj
811 dp(i,j) = dp(i,j) + dp_adj
812 ! Check for convergence to roundoff.
813 do_more(i,j) = (abs(dp_adj) > 1.0e-15*dp(i,j))
814 if (do_more(i,j)) do_any = .true.
815 endif ; enddo
816 endif
817 enddo ; endif
818 if (.not.do_any) exit
819 enddo
820
821 if (gv%semi_Boussinesq) then
822 do j=js,je ; do i=is,ie
823 h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * i_gearth
824 enddo ; enddo
825 else
826 do j=js,je ; do i=is,ie
827 h(i,j,k) = dp(i,j) * i_gearth
828 enddo ; enddo
829 endif
830 enddo
831 endif
832
833end subroutine dz_to_thickness_eos
834
835!> Converts thickness from geometric height units to thickness units, perhaps using
836!! a simple conversion factor that may be problematic in non-Boussinesq mode.
837subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode)
838 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
839 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
840 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
841 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
842 intent(in) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
843 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
844 intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2].
845 !! This is essentially intent out, but declared as intent
846 !! inout to preserve any initialized values in halo points.
847 integer, optional, intent(in) :: halo_size !< Width of halo within which to
848 !! calculate thicknesses
849 logical, optional, intent(in) :: layer_mode !< If present and true, do the conversion that
850 !! is appropriate in pure isopycnal layer mode with
851 !! no state variables or equation of state. Otherwise
852 !! use a simple constant rescaling factor and avoid the
853 !! use of GV%Rlay.
854 ! Local variables
855 logical :: layered ! If true and the model is non-Boussinesq, do calculations appropriate for use
856 ! in pure isopycnal layered mode with no state variables or equation of state.
857 integer :: i, j, k, is, ie, js, je, halo, nz
858
859 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
860 layered = .false. ; if (present(layer_mode)) layered = layer_mode
861 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
862
863 if (gv%Boussinesq) then
864 do k=1,nz ; do j=js,je ; do i=is,ie
865 h(i,j,k) = gv%Z_to_H * dz(i,j,k)
866 enddo ; enddo ; enddo
867 elseif (layered) then
868 do k=1,nz ; do j=js,je ; do i=is,ie
869 h(i,j,k) = (gv%RZ_to_H * gv%Rlay(k)) * dz(i,j,k)
870 enddo ; enddo ; enddo
871 else
872 do k=1,nz ; do j=js,je ; do i=is,ie
873 h(i,j,k) = (us%Z_to_m * gv%m_to_H) * dz(i,j,k)
874 enddo ; enddo ; enddo
875 endif
876
877end subroutine dz_to_thickness_simple
878
879!> Converts layer thicknesses in thickness units to the vertical distance between edges in height
880!! units, perhaps by multiplication by the precomputed layer-mean specific volume stored in an
881!! array in the thermo_var_ptrs type when in non-Boussinesq mode.
882subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size)
883 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
884 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
885 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
886 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
887 intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2].
888 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
889 !! thermodynamic variables
890 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
891 intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
892 !! This is essentially intent out, but declared as intent
893 !! inout to preserve any initialized values in halo points.
894 integer, optional, intent(in) :: halo_size !< Width of halo within which to
895 !! calculate thicknesses
896 ! Local variables
897 character(len=128) :: mesg ! A string for error messages
898 integer :: i, j, k, is, ie, js, je, halo, nz
899
900 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
901 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
902
903 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
904 if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then
905 if (tv%valid_SpV_halo < 0) then
906 mesg = "invalid values of SpV_avg."
907 else
908 write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
909 tv%valid_SpV_halo, halo
910 endif
911 call mom_error(fatal, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg))
912 endif
913
914 do k=1,nz ; do j=js,je ; do i=is,ie
915 dz(i,j,k) = gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
916 enddo ; enddo ; enddo
917 else
918 do k=1,nz ; do j=js,je ; do i=is,ie
919 dz(i,j,k) = gv%H_to_Z * h(i,j,k)
920 enddo ; enddo ; enddo
921 endif
922
923end subroutine thickness_to_dz_3d
924
925
926!> Converts a vertical i- / k- slice of layer thicknesses in thickness units to the vertical
927!! distance between edges in height units, perhaps by multiplication by the precomputed layer-mean
928!! specific volume stored in an array in the thermo_var_ptrs type when in non-Boussinesq mode.
929subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size)
930 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
931 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
932 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
933 intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2].
934 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
935 !! thermodynamic variables
936 real, dimension(SZI_(G),SZK_(GV)), &
937 intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
938 !! This is essentially intent out, but declared as intent
939 !! inout to preserve any initialized values in halo points.
940 integer, intent(in) :: j !< The second (j-) index of the input thicknesses to work with
941 integer, optional, intent(in) :: halo_size !< Width of halo within which to
942 !! calculate thicknesses
943 ! Local variables
944 character(len=128) :: mesg ! A string for error messages
945 integer :: i, k, is, ie, halo, nz
946
947 halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
948 is = g%isc-halo ; ie = g%iec+halo ; nz = gv%ke
949
950 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
951 if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then
952 if (tv%valid_SpV_halo < 0) then
953 mesg = "invalid values of SpV_avg."
954 else
955 write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
956 tv%valid_SpV_halo, halo
957 endif
958 call mom_error(fatal, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg))
959 endif
960
961 do k=1,nz ; do i=is,ie
962 dz(i,k) = gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
963 enddo ; enddo
964 else
965 do k=1,nz ; do i=is,ie
966 dz(i,k) = gv%H_to_Z * h(i,j,k)
967 enddo ; enddo
968 endif
969
970end subroutine thickness_to_dz_jslice
971
972
973!> Convert mixed layer depths in height units into the thickness of water in the mixed
974!! in thickness units.
975subroutine convert_mld_to_ml_thickness(MLD_in, h, h_MLD, tv, G, GV, halo)
976 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
977 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
978 real, dimension(SZI_(G),SZJ_(G)), &
979 intent(in) :: mld_in !< Input mixed layer depth [Z ~> m].
980 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
981 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
982 real, dimension(SZI_(G),SZJ_(G)), &
983 intent(out) :: h_mld !< Thickness of water in the mixed layer [H ~> m or kg m-2]
984 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
985 !! thermodynamic fields.
986 integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil
987
988 ! Local variables
989 real :: mld_rem(szi_(g)) ! The vertical extent of the MLD_in that has not yet been accounted for [Z ~> m]
990 character(len=128) :: mesg ! A string for error messages
991 logical :: keep_going
992 integer :: i, j, k, is, ie, js, je, nz, halos
993
994 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
995
996 halos = 0 ; if (present(halo)) halos = halo
997 if (present(halo)) then
998 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
999 endif
1000
1001 if (gv%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
1002 do j=js,je ; do i=is,ie
1003 h_mld(i,j) = gv%Z_to_H * mld_in(i,j)
1004 enddo ; enddo
1005 else ! The fully non-Boussinesq conversion between height in MLD_in and thickness.
1006 if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halos)) then
1007 if (tv%valid_SpV_halo < 0) then
1008 mesg = "invalid values of SpV_avg."
1009 else
1010 write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
1011 tv%valid_SpV_halo, halos
1012 endif
1013 call mom_error(fatal, "convert_MLD_to_ML_thickness called in fully non-Boussinesq mode with "//trim(mesg))
1014 endif
1015
1016 do j=js,je
1017 do i=is,ie ; mld_rem(i) = mld_in(i,j) ; h_mld(i,j) = 0.0 ; enddo
1018 do k=1,nz
1019 keep_going = .false.
1020 do i=is,ie ; if (mld_rem(i) > 0.0) then
1021 if (mld_rem(i) > gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)) then
1022 h_mld(i,j) = h_mld(i,j) + h(i,j,k)
1023 mld_rem(i) = mld_rem(i) - gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
1024 keep_going = .true.
1025 else
1026 h_mld(i,j) = h_mld(i,j) + gv%RZ_to_H * mld_rem(i) / tv%SpV_avg(i,j,k)
1027 mld_rem(i) = 0.0
1028 endif
1029 endif ; enddo
1030 if (.not.keep_going) exit
1031 enddo
1032 enddo
1033 endif
1034
1035end subroutine convert_mld_to_ml_thickness
1036
1037end module mom_interface_heights