coord_hycom.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!> Regrid columns for the HyCOM coordinate
6module coord_hycom
7
8use mom_error_handler, only : mom_error, is_root_pe, fatal, note
9use mom_variables, only : ocean_grid_type, thermo_var_ptrs
10use mom_eos, only : eos_type, calculate_density
13use regrid_interp, only : degree_max
14
15implicit none ; private
16
17#include <MOM_memory.h>
18
19!> Control structure containing required parameters for the HyCOM coordinate
20type, public :: hycom_cs ; private
21
22 !> Number of layers/levels in generated grid
23 integer :: nk
24
25 !> Nominal near-surface resolution [Z ~> m]
26 real, allocatable, dimension(:) :: coordinateresolution
27
28 !> Nominal density of interfaces [R ~> kg m-3]
29 real, allocatable, dimension(:) :: target_density
30
31 !> Maximum depths of interfaces [H ~> m or kg m-2]
32 real, allocatable, dimension(:) :: max_interface_depths
33
34 !> Maximum thicknesses of layers [H ~> m or kg m-2]
35 real, allocatable, dimension(:) :: max_layer_thickness
36
37 !> If true, an interface only moves if it improves the density fit
38 logical :: only_improves = .false.
39
40 !> If true, use 3-D control fields
41 logical :: use_3d = .false.
42
43 !> Nominal density of interfaces [R ~> kg m-3]
44 real, allocatable, dimension(:,:,:) :: target_density_3d
45
46 !> Nominal near-surface resolution [Z ~> m]
47 real, allocatable, dimension(:,:,:) :: coordinateresolution_3d
48
49 !> Interpolation control structure
50 type(interp_cs_type) :: interp_cs
51end type hycom_cs
52
54
55contains
56
57!> Initialise a hycom_CS with pointers to parameters
58subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS)
59 type(hycom_cs), pointer :: cs !< Unassociated pointer to hold the control structure
60 integer, intent(in) :: nk !< Number of layers in generated grid
61 real, dimension(nk), intent(in) :: coordinateresolution !< Nominal near-surface resolution [Z ~> m]
62 real, dimension(nk+1),intent(in) :: target_density !< Interface target densities [R ~> kg m-3]
63 type(interp_cs_type), intent(in) :: interp_cs !< Controls for interpolation
64
65 if (associated(cs)) call mom_error(fatal, "init_coord_hycom: CS already associated!")
66 allocate(cs)
67 allocate(cs%coordinateResolution(nk))
68 allocate(cs%target_density(nk+1))
69
70 cs%nk = nk
71 cs%coordinateResolution(:) = coordinateresolution(:)
72 cs%target_density(:) = target_density(:)
73 cs%use_3d = .false.
74 cs%interp_CS = interp_cs
75
76 if (is_root_pe()) call mom_error(note, "init_coord_hycom: use_3d = .false.")
77
78end subroutine init_coord_hycom
79
80!> Initialise a hycom_CS with pointers to parameters
81subroutine init_3d_coord_hycom(CS, G, nk, coordinateResolution, target_density, interp_CS)
82 type(hycom_cs), pointer :: cs !< Unassociated pointer to hold the control structure
83 type(ocean_grid_type),intent(in) :: g !< Ocean grid structure
84 integer, intent(in) :: nk !< Number of layers in generated grid
85 real, dimension(SZI_(G),SZJ_(G),nk), intent(in) :: coordinateresolution !< Nominal near-surface resolution [Z ~> m]
86 real, dimension(SZI_(G),SZJ_(G),nk+1), intent(in) :: target_density !< Interface target densities [R ~> kg m-3]
87 type(interp_cs_type), intent(in) :: interp_cs !< Controls for interpolation
88 ! Local variables
89 integer :: i,j,k
90
91 if (associated(cs)) call mom_error(fatal, "init_3d_coord_hycom: CS already associated!")
92
93 allocate(cs)
94 allocate(cs%coordinateResolution_3d(nk,szi_(g),szj_(g)), source=0.0)
95 allocate(cs%target_density_3d(nk+1,szi_(g),szj_(g)), source=0.0)
96
97 cs%nk = nk
98 cs%use_3d = .true.
99 cs%interp_CS = interp_cs
100
101 do i=g%isc-1,g%iec+1 ; do j=g%jsc-1,g%jec+1
102 if (g%mask2dT(i,j)>0.) then
103 do k= 1,nk
104 cs%coordinateResolution_3d(k,i,j) = coordinateresolution(i,j,k)
105 cs%target_density_3d(k,i,j) = target_density(i,j,k)
106 enddo
107 cs%target_density_3d(nk+1,i,j) = target_density(i,j,nk+1)
108 endif !mask2dT
109 enddo ; enddo
110
111 if (is_root_pe()) call mom_error(note, "init_3d_coord_hycom: use_3d = .true.")
112
113end subroutine init_3d_coord_hycom
114
115!> This subroutine deallocates memory in the control structure for the coord_hycom module
116subroutine end_coord_hycom(CS)
117 type(hycom_cs), pointer :: cs !< Coordinate control structure
118
119 ! nothing to do
120 if (.not. associated(cs)) return
121
122 if (allocated(cs%coordinateResolution)) deallocate(cs%coordinateResolution)
123 if (allocated(cs%target_density)) deallocate(cs%target_density)
124 if (allocated(cs%coordinateResolution_3d)) deallocate(cs%coordinateResolution_3d)
125 if (allocated(cs%target_density_3d)) deallocate(cs%target_density_3d)
126 if (allocated(cs%max_interface_depths)) deallocate(cs%max_interface_depths)
127 if (allocated(cs%max_layer_thickness)) deallocate(cs%max_layer_thickness)
128 deallocate(cs)
129end subroutine end_coord_hycom
130
131!> This subroutine can be used to set the parameters for the coord_hycom module
132subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, only_improves, interp_CS)
133 type(hycom_cs), pointer :: cs !< Coordinate control structure
134 real, dimension(:), optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces [H ~> m or kg m-2]
135 real, dimension(:), optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers [H ~> m or kg m-2]
136 logical, optional, intent(in) :: only_improves !< If true, an interface only moves if it improves the density fit
137 type(interp_cs_type), optional, intent(in) :: interp_cs !< Controls for interpolation
138
139 if (.not. associated(cs)) call mom_error(fatal, "set_hycom_params: CS not associated")
140
141 if (present(max_interface_depths)) then
142 if (size(max_interface_depths) /= cs%nk+1) &
143 call mom_error(fatal, "set_hycom_params: max_interface_depths inconsistent size")
144 allocate(cs%max_interface_depths(cs%nk+1))
145 cs%max_interface_depths(:) = max_interface_depths(:)
146 endif
147
148 if (present(max_layer_thickness)) then
149 if (size(max_layer_thickness) /= cs%nk) &
150 call mom_error(fatal, "set_hycom_params: max_layer_thickness inconsistent size")
151 allocate(cs%max_layer_thickness(cs%nk))
152 cs%max_layer_thickness(:) = max_layer_thickness(:)
153 endif
154
155 if (present(only_improves)) cs%only_improves = only_improves
156
157 if (present(interp_cs)) cs%interp_CS = interp_cs
158end subroutine set_hycom_params
159
160!> Build a HyCOM coordinate column
161subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, ix, jy, depth, h, T, S, p_col, &
162 z_col, z_col_new, zScale, h_neglect, h_neglect_edge)
163 type(hycom_cs), intent(in) :: cs !< Coordinate control structure
164 type(remapping_cs), intent(in) :: remapcs !< Remapping parameters and options
165 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
166 integer, intent(in) :: nz !< Number of levels
167 integer, intent(in) :: ix !< x direction array index
168 integer, intent(in) :: jy !< y direction array index
169 real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2])
170 real, dimension(nz), intent(in) :: t !< Temperature of column [C ~> degC]
171 real, dimension(nz), intent(in) :: s !< Salinity of column [S ~> ppt]
172 real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
173 real, dimension(nz), intent(in) :: p_col !< Layer pressure [R L2 T-2 ~> Pa]
174 real, dimension(nz+1), intent(in) :: z_col !< Interface positions relative to the surface [H ~> m or kg m-2]
175 real, dimension(CS%nk+1), intent(inout) :: z_col_new !< Absolute positions of interfaces [H ~> m or kg m-2]
176 real, optional, intent(in) :: zscale !< Scaling factor from the input coordinate thicknesses in [Z ~> m]
177 !! to desired units for zInterface, perhaps GV%Z_to_H in which
178 !! case this has units of [H Z-1 ~> nondim or kg m-3]
179 real, intent(in) :: h_neglect !< A negligibly small width for the purpose of
180 !! cell reconstruction [H ~> m or kg m-2]
181 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of
182 !! edge value calculation [H ~> m or kg m-2]
183
184 ! Local variables
185 integer :: k
186 real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3]
187 real, dimension(CS%nk) :: h_col_new ! New layer thicknesses [H ~> m or kg m-2]
188 real, dimension(CS%nk) :: r_col_new ! New layer densities [R ~> kg m-3]
189 real, dimension(CS%nk) :: t_col_new ! New layer temperatures [C ~> degC]
190 real, dimension(CS%nk) :: s_col_new ! New layer salinities [S ~> ppt]
191 real, dimension(CS%nk) :: p_col_new ! New layer pressure [R L2 T-2 ~> Pa]
192 real, dimension(CS%nk+1) :: ria_ini ! Initial nk+1 interface density anomaly w.r.t. the
193 ! interface target densities [R ~> kg m-3]
194 real, dimension(CS%nk+1) :: ria_new ! New interface density anomaly w.r.t. the
195 ! interface target densities [R ~> kg m-3]
196 real :: z_1, z_nz ! mid point of 1st and last layers [H ~> m or kg m-2]
197 real :: z_scale ! A scaling factor from the input thicknesses to the target thicknesses,
198 ! perhaps 1 or a factor in [H Z-1 ~> 1 or kg m-3]
199 real :: stretching ! z* stretching, converts z* to z [nondim].
200 real :: nominal_z ! Nominal depth of interface when using z* [H ~> m or kg m-2]
201 logical :: maximum_depths_set ! If true, the maximum depths of interface have been set.
202 logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set.
203
204 maximum_depths_set = allocated(cs%max_interface_depths)
205 maximum_h_set = allocated(cs%max_layer_thickness)
206
207 z_scale = 1.0 ; if (present(zscale)) z_scale = zscale
208
209 if (cs%only_improves .and. nz == cs%nk) then
210 call build_hycom1_target_anomaly(cs, remapcs, eqn_of_state, cs%nk, ix, jy, depth, &
211 h, t, s, p_col, rho_col, ria_ini, h_neglect, h_neglect_edge)
212 else
213 ! Work bottom recording potential density
214 call calculate_density(t, s, p_col, rho_col, eqn_of_state)
215 ! This ensures the potential density profile is monotonic
216 ! although not necessarily single valued.
217 do k = nz-1, 1, -1
218 rho_col(k) = min( rho_col(k), rho_col(k+1) )
219 enddo
220 endif
221
222 ! Interpolates for the target interface position with the rho_col profile
223 ! Based on global density profile, interpolate to generate a new grid
224 if (cs%use_3d) then
225 call build_and_interpolate_grid(cs%interp_CS, rho_col, nz, h(:), z_col, &
226 cs%target_density_3d(:,ix,jy), cs%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge)
227 else
228 call build_and_interpolate_grid(cs%interp_CS, rho_col, nz, h(:), z_col, &
229 cs%target_density, cs%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge)
230 endif
231 if (cs%only_improves .and. nz == cs%nk) then
232 ! Only move an interface if it improves the density fit
233 z_1 = 0.5 * ( z_col(1) + z_col(2) )
234 z_nz = 0.5 * ( z_col(nz) + z_col(nz+1) )
235 do k = 1,cs%nk
236 p_col_new(k) = p_col(1) + ( 0.5 * ( z_col_new(k) + z_col_new(k+1) ) - z_1 ) &
237 / ( z_nz - z_1 ) * ( p_col(nz) - p_col(1) )
238 enddo
239 ! Remap from original h and T,S to get T,S_col_new
240 call remapping_core_h(remapcs, nz, h(:), t, cs%nk, h_col_new, t_col_new)
241 call remapping_core_h(remapcs, nz, h(:), s, cs%nk, h_col_new, s_col_new)
242 call build_hycom1_target_anomaly(cs, remapcs, eqn_of_state, cs%nk, ix, jy, depth, &
243 h_col_new, t_col_new, s_col_new, p_col_new, r_col_new, ria_new, h_neglect, h_neglect_edge)
244 do k= 2,cs%nk
245 if ( abs(ria_ini(k)) <= abs(ria_new(k)) .and. z_col(k) > z_col_new(k-1) .and. &
246 z_col(k) < z_col_new(k+1)) then
247 z_col_new(k) = z_col(k)
248 endif
249 enddo
250 endif !only_improves
251
252 ! Sweep down the interfaces and make sure that the interface is at least
253 ! as deep as a nominal target z* grid
254 nominal_z = 0.
255 stretching = z_col(nz+1) / depth ! Stretches z* to z
256 if (cs%use_3d) then
257 do k = 2, cs%nk+1
258 nominal_z = nominal_z + (z_scale * cs%coordinateResolution_3d(k-1,ix,jy)) * stretching
259 z_col_new(k) = max( z_col_new(k), nominal_z )
260 z_col_new(k) = min( z_col_new(k), z_col(nz+1) )
261 enddo
262 else
263 do k = 2, cs%nk+1
264 nominal_z = nominal_z + (z_scale * cs%coordinateResolution(k-1)) * stretching
265 z_col_new(k) = max( z_col_new(k), nominal_z )
266 z_col_new(k) = min( z_col_new(k), z_col(nz+1) )
267 enddo
268 endif
269
270 if (maximum_depths_set .and. maximum_h_set) then ; do k=2,cs%nk
271 ! The loop bounds are 2 & nz so the top and bottom interfaces do not move.
272 ! Recall that z_col_new is positive downward.
273 z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k), &
274 z_col_new(k-1) + cs%max_layer_thickness(k-1))
275 enddo ; elseif (maximum_depths_set) then ; do k=2,cs%nk
276 z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k))
277 enddo ; elseif (maximum_h_set) then ; do k=2,cs%nk
278 z_col_new(k) = min(z_col_new(k), z_col_new(k-1) + cs%max_layer_thickness(k-1))
279 enddo ; endif
280end subroutine build_hycom1_column
281
282!> Calculate interface density anomaly w.r.t. the target.
283subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, ix, jy, depth, h, T, S, &
284 p_col, R, RiAnom, h_neglect, h_neglect_edge)
285 type(hycom_cs), intent(in) :: CS !< Coordinate control structure
286 type(remapping_cs), intent(in) :: remapCS !< Remapping parameters and options
287 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
288 integer, intent(in) :: nz !< Number of levels
289 integer, intent(in) :: ix !< x direction array index
290 integer, intent(in) :: jy !< y direction array index
291 real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2])
292 real, dimension(nz), intent(in) :: T !< Temperature of column [C ~> degC]
293 real, dimension(nz), intent(in) :: S !< Salinity of column [S ~> ppt]
294 real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
295 real, dimension(nz), intent(in) :: p_col !< Layer pressure [R L2 T-2 ~> Pa]
296 real, dimension(nz), intent(out) :: R !< Layer density [R ~> kg m-3]
297 real, dimension(nz+1), intent(out) :: RiAnom !< The interface density anomaly
298 !! w.r.t. the interface target
299 !! densities [R ~> kg m-3]
300 real, intent(in) :: h_neglect !< A negligibly small width for the purpose of
301 !! cell reconstruction [H ~> m or kg m-2]
302 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of
303 !! edge value calculation [H ~> m or kg m-2]
304 ! Local variables
305 integer :: degree,k
306 real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3]
307 real, dimension(nz,2) :: ppoly_E ! Polynomial edge values [R ~> kg m-3]
308 real, dimension(nz,2) :: ppoly_S ! Polynomial edge slopes [R H-1]
309 real, dimension(nz,DEGREE_MAX+1) :: ppoly_C ! Polynomial interpolant coeficients on the local 0-1 grid [R ~> kg m-3]
310
311 ! Work bottom recording potential density
312 call calculate_density(t, s, p_col, rho_col, eqn_of_state)
313 ! This ensures the potential density profile is monotonic
314 ! although not necessarily single valued.
315 do k = nz-1, 1, -1
316 rho_col(k) = min( rho_col(k), rho_col(k+1) )
317 enddo
318
319 call regridding_set_ppolys(cs%interp_CS, rho_col, nz, h, ppoly_e, ppoly_s, ppoly_c, &
320 degree, h_neglect, h_neglect_edge)
321
322 if (cs%use_3d) then
323 r(1) = rho_col(1)
324 rianom(1) = ppoly_e(1,1) - cs%target_density_3d(1,ix,jy)
325 do k= 2,nz
326 r(k) = rho_col(k)
327 if (ppoly_e(k-1,2) > cs%target_density_3d(k,ix,jy)) then
328 rianom(k) = ppoly_e(k-1,2) - cs%target_density_3d(k,ix,jy) !interface is heavier than target
329 elseif (ppoly_e(k,1) < cs%target_density_3d(k,ix,jy)) then
330 rianom(k) = ppoly_e(k,1) - cs%target_density_3d(k,ix,jy) !interface is lighter than target
331 else
332 rianom(k) = 0.0 !interface spans the target
333 endif
334 enddo
335 rianom(nz+1) = ppoly_e(nz,2) - cs%target_density_3d(nz+1,ix,jy)
336 else
337 r(1) = rho_col(1)
338 rianom(1) = ppoly_e(1,1) - cs%target_density(1)
339 do k= 2,nz
340 r(k) = rho_col(k)
341 if (ppoly_e(k-1,2) > cs%target_density(k)) then
342 rianom(k) = ppoly_e(k-1,2) - cs%target_density(k) !interface is heavier than target
343 elseif (ppoly_e(k,1) < cs%target_density(k)) then
344 rianom(k) = ppoly_e(k,1) - cs%target_density(k) !interface is lighter than target
345 else
346 rianom(k) = 0.0 !interface spans the target
347 endif
348 enddo
349 rianom(nz+1) = ppoly_e(nz,2) - cs%target_density(nz+1)
350 endif !use_3d:else
351
352end subroutine build_hycom1_target_anomaly
353
354end module coord_hycom