coord_zlike.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 a z-like coordinate (z-star, z-level)
6module coord_zlike
7
8use mom_error_handler, only : mom_error, fatal
9
10implicit none ; private
11
12!> Control structure containing required parameters for a z-like coordinate
13type, public :: zlike_cs ; private
14
15 !> Number of levels to be generated
16 integer :: nk
17
18 !> Minimum thickness allowed for layers, in the same thickness units (perhaps [H ~> m or kg m-2])
19 !! that will be used in all subsequent calls to build_zstar_column with this structure.
20 real :: min_thickness
21
22 !> Target coordinate resolution, usually in [Z ~> m]
23 real, allocatable, dimension(:) :: coordinateresolution
24end type zlike_cs
25
27
28contains
29
30!> Initialise a zlike_CS with pointers to parameters
31subroutine init_coord_zlike(CS, nk, coordinateResolution)
32 type(zlike_cs), pointer :: cs !< Unassociated pointer to hold the control structure
33 integer, intent(in) :: nk !< Number of levels in the grid
34 real, dimension(:), intent(in) :: coordinateresolution !< Target coordinate resolution [Z ~> m]
35
36 if (associated(cs)) call mom_error(fatal, "init_coord_zlike: CS already associated!")
37 allocate(cs)
38 allocate(cs%coordinateResolution(nk))
39
40 cs%nk = nk
41 cs%coordinateResolution = coordinateresolution
42end subroutine init_coord_zlike
43
44!> Deallocates the zlike control structure
45subroutine end_coord_zlike(CS)
46 type(zlike_cs), pointer :: cs !< Coordinate control structure
47
48 ! Nothing to do
49 if (.not. associated(cs)) return
50 deallocate(cs%coordinateResolution)
51 deallocate(cs)
52end subroutine end_coord_zlike
53
54!> Set parameters in the zlike structure
55subroutine set_zlike_params(CS, min_thickness)
56 type(zlike_cs), pointer :: cs !< Coordinate control structure
57 real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
58
59 if (.not. associated(cs)) call mom_error(fatal, "set_zlike_params: CS not associated")
60
61 if (present(min_thickness)) cs%min_thickness = min_thickness
62end subroutine set_zlike_params
63
64!> Builds a z* coordinate with a minimum thickness
65subroutine build_zstar_column(CS, depth, total_thickness, zInterface, &
66 z_rigid_top, eta_orig, zScale)
67 type(zlike_cs), intent(in) :: cs !< Coordinate control structure
68 real, intent(in) :: depth !< Depth of ocean bottom (positive downward in the
69 !! output units), units may be [Z ~> m] or [H ~> m or kg m-2]
70 real, intent(in) :: total_thickness !< Column thickness (positive definite in the same
71 !! units as depth) [Z ~> m] or [H ~> m or kg m-2]
72 real, dimension(CS%nk+1), intent(inout) :: zinterface !< Absolute positions of interfaces (in the same
73 !! units as depth) [Z ~> m] or [H ~> m or kg m-2]
74 real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (positive upward in the same
75 !! units as depth) [Z ~> m] or [H ~> m or kg m-2]
76 real, optional, intent(in) :: eta_orig !< The actual original height of the top (in the same
77 !! units as depth) [Z ~> m] or [H ~> m or kg m-2]
78 real, optional, intent(in) :: zscale !< Scaling factor from the target coordinate resolution
79 !! in Z to desired units for zInterface, perhaps Z_to_H,
80 !! often [nondim] or [H Z-1 ~> 1 or kg m-3]
81 ! Local variables
82 real :: eta ! Free surface height [Z ~> m] or [H ~> m or kg m-2]
83 real :: stretching ! A stretching factor for the coordinate [nondim]
84 real :: dh, min_thickness, z0_top, z_star, z_scale ! Thicknesses or heights [Z ~> m] or [H ~> m or kg m-2]
85 integer :: k
86 logical :: new_zstar_def
87
88 z_scale = 1.0 ; if (present(zscale)) z_scale = zscale
89
90 new_zstar_def = .false.
91 min_thickness = min( cs%min_thickness, total_thickness/real(cs%nk) )
92 z0_top = 0.
93 if (present(z_rigid_top)) then
94 z0_top = z_rigid_top
95 new_zstar_def = .true.
96 endif
97
98 ! Position of free-surface (or the rigid top, for which eta ~ z0_top)
99 eta = total_thickness - depth
100 if (present(eta_orig)) eta = eta_orig
101
102 ! Conventional z* coordinate:
103 ! z* = (z-eta) / stretching where stretching = (H+eta)/H
104 ! z = eta + stretching * z*
105 ! The above gives z*(z=eta) = 0, z*(z=-H) = -H.
106 ! With a rigid top boundary at eta = z0_top then
107 ! z* = z0 + (z-eta) / stretching where stretching = (H+eta)/(H+z0)
108 ! z = eta + stretching * (z*-z0) * stretching
109 stretching = total_thickness / ( depth + z0_top )
110
111 if (new_zstar_def) then
112 ! z_star is the notional z* coordinate in absence of upper/lower topography
113 z_star = 0. ! z*=0 at the free-surface
114 zinterface(1) = eta ! The actual position of the top of the column
115 do k = 2,cs%nk
116 z_star = z_star - cs%coordinateResolution(k-1)*z_scale
117 ! This ensures that z is below a rigid upper surface (ice shelf bottom)
118 zinterface(k) = min( eta + stretching * ( z_star - z0_top ), z0_top )
119 ! This ensures that the layer in inflated
120 zinterface(k) = min( zinterface(k), zinterface(k-1) - min_thickness )
121 ! This ensures that z is above or at the topography
122 zinterface(k) = max( zinterface(k), -depth + real(cs%nk+1-k) * min_thickness )
123 enddo
124 zinterface(cs%nk+1) = -depth
125
126 else
127 ! Integrate down from the top for a notional new grid, ignoring topography
128 ! The starting position is offset by z0_top which, if z0_top<0, will place
129 ! interfaces above the rigid boundary.
130 zinterface(1) = eta
131 do k = 1,cs%nk
132 dh = stretching * cs%coordinateResolution(k)*z_scale ! Notional grid spacing
133 zinterface(k+1) = zinterface(k) - dh
134 enddo
135
136 ! Integrating up from the bottom adjusting interface position to accommodate
137 ! inflating layers without disturbing the interface above
138 zinterface(cs%nk+1) = -depth
139 do k = cs%nk,1,-1
140 if ( zinterface(k) < (zinterface(k+1) + min_thickness) ) then
141 zinterface(k) = zinterface(k+1) + min_thickness
142 endif
143 enddo
144 endif
145
146end subroutine build_zstar_column
147
148end module coord_zlike