MOM_verticalGrid.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!> Provides a transparent vertical ocean grid type and supporting routines
7
8use mom_error_handler, only : mom_error, mom_mesg, fatal
9use mom_file_parser, only : get_param, log_param, log_version, param_file_type
11
12implicit none ; private
13
14#include <MOM_memory.h>
15
19
20! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
21! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
22! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
23! vary with the Boussinesq approximation, the Boussinesq variant is given first.
24
25!> Describes the vertical ocean grid, including unit conversion factors
26type, public :: verticalgrid_type
27
28 ! Commonly used parameters
29 integer :: ke !< The number of layers/levels in the vertical
30 real :: max_depth !< The maximum depth of the ocean [Z ~> m].
31! real :: mks_g_Earth !< The gravitational acceleration in unscaled MKS units [m s-2]. This might not be used.
32 real :: g_earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
33 real :: g_earth_z_t2 !< The gravitational acceleration in alternatively rescaled units [Z T-2 ~> m s-2]
34 real :: rho0 !< The density used in the Boussinesq approximation or nominal
35 !! density used to convert depths into mass units [R ~> kg m-3].
36
37 ! Vertical coordinate descriptions for diagnostics and I/O
38 character(len=40) :: zaxisunits !< The units that vertical coordinates are written in
39 character(len=40) :: zaxislongname !< Coordinate name to appear in files,
40 !! e.g. "Target Potential Density" or "Height"
41 real, allocatable, dimension(:) :: slayer !< Coordinate values of layer centers, in unscaled
42 !! units that depend on the vertical coordinate, such as [kg m-3] for an
43 !! isopycnal or some hybrid coordinates, [m] for a Z* coordinate,
44 !! or [nondim] for a sigma coordinate.
45 real, allocatable, dimension(:) :: sinterface !< Coordinate values on interfaces, in the same
46 !! unscale units as sLayer [various].
47 integer :: direction = 1 !< Direction defaults to 1, positive up.
48
49 ! The following variables give information about the vertical grid.
50 logical :: boussinesq !< If true, make the Boussinesq approximation.
51 logical :: semi_boussinesq !< If true, do non-Boussinesq pressure force calculations and
52 !! use mass-based "thicknesses, but use Rho0 to convert layer thicknesses
53 !! into certain height changes. This only applies if BOUSSINESQ is false.
54 real :: angstrom_h !< A one-Angstrom thickness in the model thickness units [H ~> m or kg m-2].
55 real :: angstrom_z !< A one-Angstrom thickness in the model depth units [Z ~> m].
56 real :: angstrom_m !< A one-Angstrom thickness [m].
57 real :: h_subroundoff !< A thickness that is so small that it can be added to a thickness of
58 !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2].
59 !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m.
60 real :: dz_subroundoff !< A thickness in height units that is so small that it can be added to a
61 !! vertical distance of Angstrom_Z or 1e-17 m without changing it at the bit
62 !! level [Z ~> m]. This is the height equivalent of H_subroundoff.
63 real, allocatable, dimension(:) :: &
64 g_prime, & !< The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2].
65 rlay !< The target coordinate value (potential density) in each layer [R ~> kg m-3].
66 integer :: nkml = 0 !< The number of layers at the top that should be treated
67 !! as parts of a homogeneous region.
68 integer :: nk_rho_varies = 0 !< The number of layers at the top where the
69 !! density does not track any target density.
70 real :: h_to_kg_m2 !< A constant that translates thicknesses from the units of thickness
71 !! to kg m-2 [kg m-2 H-1 ~> kg m-3 or 1].
72 real :: kg_m2_to_h !< A constant that translates thicknesses from kg m-2 to the units
73 !! of thickness [H m2 kg-1 ~> m3 kg-1 or 1].
74 real :: m_to_h !< A constant that translates distances in m to the units of
75 !! thickness [H m-1 ~> 1 or kg m-3].
76 real :: h_to_m !< A constant that translates distances in the units of thickness
77 !! to m [m H-1 ~> 1 or m3 kg-1].
78 real :: h_to_pa !< A constant that translates the units of thickness to pressure
79 !! [Pa H-1 ~> kg m-2 s-2 or m s-2].
80 real :: h_to_z !< A constant that translates thickness units to the units of
81 !! depth [Z H-1 ~> 1 or m3 kg-1].
82 real :: z_to_h !< A constant that translates depth units to thickness units
83 !! depth [H Z-1 ~> 1 or kg m-3].
84 real :: h_to_rz !< A constant that translates thickness units to the units of
85 !! mass per unit area [R Z H-1 ~> kg m-3 or 1].
86 real :: rz_to_h !< A constant that translates mass per unit area units to
87 !! thickness units [H R-1 Z-1 ~> m3 kg-2 or 1].
88 real :: h_to_mks !< A constant that translates thickness units to its MKS unit
89 !! (m or kg m-2) based on GV%Boussinesq [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
90 real :: m2_s_to_hz_t !< The combination of conversion factors that converts kinematic viscosities
91 !! in m2 s-1 to the internal units of the kinematic (in Boussinesq mode)
92 !! or dynamic viscosity [H Z s T-1 m-2 ~> 1 or kg m-3]
93 real :: hz_t_to_m2_s !< The combination of conversion factors that converts the viscosities from
94 !! their internal representation into a kinematic viscosity in m2 s-1
95 !! [T m2 H-1 Z-1 s-1 ~> 1 or m3 kg-1]
96 real :: hz_t_to_mks !< The combination of conversion factors that converts the viscosities from
97 !! their internal representation into their unnscaled MKS units
98 !! (m2 s-1 or Pa s), depending on whether the model is Boussinesq
99 !! [T m2 H-1 Z-1 s-1 ~> 1] or [T Pa s H-1 Z-1 ~> 1]
100
101end type verticalgrid_type
102
103contains
104
105!> Allocates and initializes the ocean model vertical grid structure.
106subroutine verticalgridinit( param_file, GV, US )
107 type(param_file_type), intent(in) :: param_file !< Parameter file handle/type
108 type(verticalgrid_type), pointer :: gv !< The container for vertical grid data
109 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
110 ! This routine initializes the verticalGrid_type structure (GV).
111 ! All memory is allocated but not necessarily set to meaningful values until later.
112
113 ! Local variables
114 integer :: nk, h_power
115 real :: h_rescale_factor ! The integer power of 2 by which thicknesses are rescaled [nondim]
116 real :: rho_kv ! The density used convert input kinematic viscosities into dynamic viscosities
117 ! when in non-Boussinesq mode [R ~> kg m-3]
118 ! This include declares and sets the variable "version".
119# include "version_variable.h"
120 character(len=16) :: mdl = 'MOM_verticalGrid'
121
122 if (associated(gv)) call mom_error(fatal, &
123 'verticalGridInit: called with an associated GV pointer.')
124 allocate(gv)
125
126 ! Read all relevant parameters and write them to the model log.
127 call log_version(param_file, mdl, version, &
128 "Parameters providing information about the vertical grid.", &
129 log_to_all=.true., debugging=.true.)
130 call get_param(param_file, mdl, "G_EARTH", gv%g_Earth, &
131 "The gravitational acceleration of the Earth.", &
132 units="m s-2", default=9.80, scale=us%Z_to_m*us%m_s_to_L_T**2)
133 call get_param(param_file, mdl, "RHO_0", gv%Rho0, &
134 "The mean ocean density used with BOUSSINESQ true to "//&
135 "calculate accelerations and the mass for conservation "//&
136 "properties, or with BOUSSINESQ false to convert some "//&
137 "parameters from vertical units of m to kg m-2.", &
138 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
139 call get_param(param_file, mdl, "BOUSSINESQ", gv%Boussinesq, &
140 "If true, make the Boussinesq approximation.", default=.true.)
141 call get_param(param_file, mdl, "SEMI_BOUSSINESQ", gv%semi_Boussinesq, &
142 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
143 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
144 "height changes. This only applies if BOUSSINESQ is false.", &
145 default=.true., do_not_log=gv%Boussinesq)
146 if (gv%Boussinesq) gv%semi_Boussinesq = .true.
147 call get_param(param_file, mdl, "RHO_KV_CONVERT", rho_kv, &
148 "The density used to convert input vertical distances into thickesses in "//&
149 "non-BOUSSINESQ mode, and to convert kinematic viscosities into dynamic "//&
150 "viscosities and similarly for vertical diffusivities. GV%m_to_H is set "//&
151 "using this value, whereas GV%Z_to_H is set using RHO_0. The default is "//&
152 "RHO_0, but this can be set separately to demonstrate the independence of the "//&
153 "non-Boussinesq solutions of the value of RHO_0.", &
154 units="kg m-3", default=gv%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &
155 do_not_log=gv%Boussinesq)
156 call get_param(param_file, mdl, "ANGSTROM", gv%Angstrom_Z, &
157 "The minimum layer thickness, usually one-Angstrom.", &
158 units="m", default=1.0e-10, scale=us%m_to_Z)
159 call get_param(param_file, mdl, "H_RESCALE_POWER", h_power, &
160 "An integer power of 2 that is used to rescale the model's "//&
161 "intenal units of thickness. Valid values range from -300 to 300.", &
162 units="nondim", default=0, debuggingparam=.true.)
163 if (abs(h_power) > 300) call mom_error(fatal, "verticalGridInit: "//&
164 "H_RESCALE_POWER is outside of the valid range of -300 to 300.")
165 h_rescale_factor = 1.0
166 if (h_power /= 0) h_rescale_factor = 2.0**h_power
167 if (.not.gv%Boussinesq) then
168 call get_param(param_file, mdl, "H_TO_KG_M2", gv%H_to_kg_m2,&
169 "A constant that translates thicknesses from the model's "//&
170 "internal units of thickness to kg m-2.", units="kg m-2 H-1", &
171 default=1.0)
172 gv%H_to_kg_m2 = gv%H_to_kg_m2 * h_rescale_factor
173 else
174 call get_param(param_file, mdl, "H_TO_M", gv%H_to_m, &
175 "A constant that translates the model's internal "//&
176 "units of thickness into m.", units="m H-1", default=1.0)
177 gv%H_to_m = gv%H_to_m * h_rescale_factor
178 endif
179 ! This is not used: GV%mks_g_Earth = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth
180 gv%g_Earth_Z_T2 = us%L_to_Z**2 * gv%g_Earth ! This would result from scale=US%m_to_Z*US%T_to_s**2.
181#ifdef STATIC_MEMORY_
182 ! Here NK_ is a macro, while nk is a variable.
183 call get_param(param_file, mdl, "NK", nk, &
184 "The number of model layers.", units="nondim", &
185 default=nk_)
186 if (nk /= nk_) call mom_error(fatal, "verticalGridInit: " // &
187 "Mismatched number of layers NK_ between MOM_memory.h and param_file")
188
189#else
190 call get_param(param_file, mdl, "NK", nk, &
191 "The number of model layers.", units="nondim", fail_if_missing=.true.)
192#endif
193 gv%ke = nk
194
195 if (gv%Boussinesq) then
196 gv%H_to_kg_m2 = us%R_to_kg_m3*gv%Rho0 * gv%H_to_m
197 gv%kg_m2_to_H = 1.0 / gv%H_to_kg_m2
198 gv%m_to_H = 1.0 / gv%H_to_m
199 gv%H_to_MKS = gv%H_to_m
200 gv%m2_s_to_HZ_T = gv%m_to_H * us%m_to_Z * us%T_to_s
201
202 gv%H_to_Z = gv%H_to_m * us%m_to_Z
203 gv%Z_to_H = us%Z_to_m * gv%m_to_H
204 else
205 gv%kg_m2_to_H = 1.0 / gv%H_to_kg_m2
206 ! GV%m_to_H = US%R_to_kg_m3*GV%Rho0 * GV%kg_m2_to_H
207 gv%m_to_H = us%R_to_kg_m3*rho_kv * gv%kg_m2_to_H
208 gv%H_to_MKS = gv%H_to_kg_m2
209 gv%m2_s_to_HZ_T = us%R_to_kg_m3*rho_kv * gv%kg_m2_to_H * us%m_to_Z * us%T_to_s
210 gv%H_to_m = 1.0 / gv%m_to_H
211
212 gv%H_to_Z = us%m_to_Z * ( gv%H_to_kg_m2 / (us%R_to_kg_m3*gv%Rho0) )
213 gv%Z_to_H = us%Z_to_m * ( us%R_to_kg_m3*gv%Rho0 * gv%kg_m2_to_H )
214 endif
215
216 gv%Angstrom_H = (us%Z_to_m * gv%m_to_H) * gv%Angstrom_Z
217 gv%Angstrom_m = us%Z_to_m * gv%Angstrom_Z
218
219 gv%H_subroundoff = 1e-20 * max(gv%Angstrom_H, gv%m_to_H*1e-17)
220 gv%dZ_subroundoff = 1e-20 * max(gv%Angstrom_Z, us%m_to_Z*1e-17)
221
222 gv%H_to_Pa = us%L_T_to_m_s**2*us%m_to_Z * gv%g_Earth * gv%H_to_kg_m2
223
224 gv%H_to_RZ = gv%H_to_kg_m2 * us%kg_m3_to_R * us%m_to_Z
225 gv%RZ_to_H = gv%kg_m2_to_H * us%R_to_kg_m3 * us%Z_to_m
226
227 gv%HZ_T_to_m2_s = 1.0 / gv%m2_s_to_HZ_T
228 gv%HZ_T_to_MKS = gv%H_to_MKS * us%Z_to_m * us%s_to_T
229
230 ! Note based on the above that for both Boussinsq and non-Boussinesq cases that:
231 ! GV%Rho0 = GV%Z_to_H * GV%H_to_RZ
232 ! 1.0/GV%Rho0 = GV%H_to_Z * GV%RZ_to_H
233 ! This is exact for power-of-2 scaling of the units, regardless of the value of Rho0, but
234 ! the first term on the right hand side is invertable in Boussinesq mode, but the second
235 ! is invertable when non-Boussinesq.
236
237 ! Log derivative values.
238 call log_param(param_file, mdl, "M to THICKNESS", gv%m_to_H*h_rescale_factor, units="H m-1")
239 call log_param(param_file, mdl, "M to THICKNESS rescaled by 2^-n", gv%m_to_H, units="2^n H m-1")
240 call log_param(param_file, mdl, "THICKNESS to M rescaled by 2^n", gv%H_to_m, units="2^-n m H-1")
241
242 allocate( gv%sInterface(nk+1) )
243 allocate( gv%sLayer(nk) )
244 allocate( gv%g_prime(nk+1), source=0.0 )
245 allocate( gv%Rlay(nk), source=0.0 )
246
247end subroutine verticalgridinit
248
249!> Returns the model's thickness units, usually m or kg/m^2.
250function get_thickness_units(GV)
251 character(len=48) :: get_thickness_units !< The vertical thickness units
252 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
253 ! This subroutine returns the appropriate units for thicknesses,
254 ! depending on whether the model is Boussinesq or not and the scaling for
255 ! the vertical thickness.
256
257 if (gv%Boussinesq) then
259 else
260 get_thickness_units = "kg m-2"
261 endif
262end function get_thickness_units
263
264!> Returns the model's thickness flux units, usually m^3/s or kg/s.
265function get_flux_units(GV)
266 character(len=48) :: get_flux_units !< The thickness flux units
267 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
268 ! This subroutine returns the appropriate units for thickness fluxes,
269 ! depending on whether the model is Boussinesq or not and the scaling for
270 ! the vertical thickness.
271
272 if (gv%Boussinesq) then
273 get_flux_units = "m3 s-1"
274 else
275 get_flux_units = "kg s-1"
276 endif
277end function get_flux_units
278
279!> Returns the model's tracer flux units.
280function get_tr_flux_units(GV, tr_units, tr_vol_conc_units, tr_mass_conc_units)
281 character(len=48) :: get_tr_flux_units !< The model's flux units
282 !! for a tracer.
283 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical
284 !! grid structure.
285 character(len=*), optional, intent(in) :: tr_units !< Units for a tracer, for example
286 !! Celsius or PSU.
287 character(len=*), optional, intent(in) :: tr_vol_conc_units !< The concentration units per unit
288 !! volume, for example if the units are
289 !! umol m-3, tr_vol_conc_units would
290 !! be umol.
291 character(len=*), optional, intent(in) :: tr_mass_conc_units !< The concentration units per unit
292 !! mass of sea water, for example if
293 !! the units are mol kg-1,
294 !! tr_vol_conc_units would be mol.
295
296 ! This subroutine returns the appropriate units for thicknesses and fluxes,
297 ! depending on whether the model is Boussinesq or not and the scaling for
298 ! the vertical thickness.
299 integer :: cnt
300
301 cnt = 0
302 if (present(tr_units)) cnt = cnt+1
303 if (present(tr_vol_conc_units)) cnt = cnt+1
304 if (present(tr_mass_conc_units)) cnt = cnt+1
305
306 if (cnt == 0) call mom_error(fatal, "get_tr_flux_units: One of the three "//&
307 "arguments tr_units, tr_vol_conc_units, or tr_mass_conc_units "//&
308 "must be present.")
309 if (cnt > 1) call mom_error(fatal, "get_tr_flux_units: Only one of "//&
310 "tr_units, tr_vol_conc_units, and tr_mass_conc_units may be present.")
311 if (present(tr_units)) then
312 if (gv%Boussinesq) then
313 get_tr_flux_units = trim(tr_units)//" m3 s-1"
314 else
315 get_tr_flux_units = trim(tr_units)//" kg s-1"
316 endif
317 endif
318 if (present(tr_vol_conc_units)) then
319 if (gv%Boussinesq) then
320 get_tr_flux_units = trim(tr_vol_conc_units)//" s-1"
321 else
322 get_tr_flux_units = trim(tr_vol_conc_units)//" m-3 kg s-1"
323 endif
324 endif
325 if (present(tr_mass_conc_units)) then
326 if (gv%Boussinesq) then
327 get_tr_flux_units = trim(tr_mass_conc_units)//" kg-1 m3 s-1"
328 else
329 get_tr_flux_units = trim(tr_mass_conc_units)//" s-1"
330 endif
331 endif
332
333end function get_tr_flux_units
334
335!> This sets the coordinate data for the "layer mode" of the isopycnal model.
336subroutine setverticalgridaxes( Rlay, GV, scale )
337 type(verticalgrid_type), intent(inout) :: gv !< The container for vertical grid data
338 real, dimension(GV%ke), intent(in) :: rlay !< The layer target density [R ~> kg m-3]
339 real, intent(in) :: scale !< A unit scaling factor for Rlay to convert
340 !! it into the units of sInterface, usually
341 !! [kg m-3 R-1 ~> 1] when used in layer mode.
342 ! Local variables
343 integer :: k, nk
344
345 nk = gv%ke
346
347 gv%zAxisLongName = 'Target Potential Density'
348 gv%zAxisUnits = 'kg m-3'
349 do k=1,nk ; gv%sLayer(k) = scale*rlay(k) ; enddo
350 if (nk > 1) then
351 gv%sInterface(1) = scale * (1.5*rlay(1) - 0.5*rlay(2))
352 do k=2,nk ; gv%sInterface(k) = scale * 0.5*( rlay(k-1) + rlay(k) ) ; enddo
353 gv%sInterface(nk+1) = scale * (1.5*rlay(nk) - 0.5*rlay(nk-1))
354 else
355 gv%sInterface(1) = 0.0 ; gv%sInterface(nk+1) = 2.0*scale*rlay(nk)
356 endif
357
358end subroutine setverticalgridaxes
359
360!> Deallocates the model's vertical grid structure.
361subroutine verticalgridend( GV )
362 type(verticalgrid_type), pointer :: gv !< The ocean's vertical grid structure
363
364 deallocate( gv%g_prime, gv%Rlay )
365 deallocate( gv%sInterface , gv%sLayer )
366 deallocate( gv )
367
368end subroutine verticalgridend
369
370end module mom_verticalgrid