MOM_variables.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 transparent structures with groups of MOM6 variables and supporting routines
6module mom_variables
7
8use mom_array_transform, only : rotate_array, rotate_vector
9use mom_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type
10use mom_coupler_types, only : coupler_type_spawn, coupler_type_destructor, coupler_type_initialized
11use mom_coupler_types, only : coupler_type_copy_data
12use mom_debugging, only : hchksum
13use mom_domains, only : mom_domain_type, get_domain_extent, group_pass_type
14use mom_eos, only : eos_type
15use mom_error_handler, only : mom_error, fatal
16use mom_grid, only : ocean_grid_type
20
21implicit none ; private
22
23#include <MOM_memory.h>
24
28
29! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
30! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
31! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
32! vary with the Boussinesq approximation, the Boussinesq variant is given first.
33
34!> A structure for creating arrays of pointers to 3D arrays
35type, public :: p3d
36 real, dimension(:,:,:), pointer :: p => null() !< A pointer to a 3D array [various]
37end type p3d
38!> A structure for creating arrays of pointers to 2D arrays
39type, public :: p2d
40 real, dimension(:,:), pointer :: p => null() !< A pointer to a 2D array [various]
41end type p2d
42
43!> Pointers to various fields which may be used describe the surface state of MOM, and which
44!! will be returned to the calling program
45type, public :: surface
46 real, allocatable, dimension(:,:) :: &
47 sst, & !< The sea surface temperature [C ~> degC].
48 sss, & !< The sea surface salinity [S ~> psu or gSalt/kg].
49 sfc_density, & !< The mixed layer density [R ~> kg m-3].
50 hml, & !< The mixed layer depth [Z ~> m].
51 u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1].
52 v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1].
53 sea_lev, & !< The sea level [Z ~> m]. If a reduced surface gravity is
54 !! used, that is compensated for in sea_lev.
55 frazil, & !< The energy needed to heat the ocean column to the freezing point during
56 !! the call to step_MOM [Q R Z ~> J m-2].
57 melt_potential, & !< Instantaneous amount of heat that can be used to melt sea ice [Q R Z ~> J m-2].
58 !! This is computed w.r.t. surface freezing temperature.
59 ocean_mass, & !< The total mass of the ocean [R Z ~> kg m-2].
60 ocean_heat, & !< The total heat content of the ocean in [C R Z ~> degC kg m-2].
61 ocean_salt, & !< The total salt content of the ocean in [1e-3 S R Z ~> kgSalt m-2].
62 taux_shelf, & !< The zonal stresses on the ocean under shelves [R L Z T-2 ~> Pa].
63 tauy_shelf, & !< The meridional stresses on the ocean under shelves [R L Z T-2 ~> Pa].
64 fco2 !< CO2 flux from the ocean to the atmosphere [R Z T-1 ~> kgCO2 m-2 s-1]
65 logical :: t_is_cont = .false. !< If true, the temperature variable SST is actually the
66 !! conservative temperature in [C ~> degC].
67 logical :: s_is_abss = .false. !< If true, the salinity variable SSS is actually the
68 !! absolute salinity in [S ~> gSalt kg-1].
69 type(coupler_2d_bc_type) :: tr_fields !< A structure that may contain an
70 !! array of named fields describing tracer-related quantities.
71 !### NOTE: ALL OF THE ARRAYS IN TR_FIELDS USE THE COUPLER'S INDEXING CONVENTION AND HAVE NO
72 !### HALOS! THIS IS DONE TO CONFORM TO THE TREATMENT IN MOM4, BUT I DON'T LIKE IT! -RWH
73 logical :: arrays_allocated = .false. !< A flag that indicates whether the surface type
74 !! has had its memory allocated.
75end type surface
76
77!> Pointers to an assortment of thermodynamic fields that may be available, including
78!! potential temperature, salinity, heat capacity, and the equation of state control structure.
79type, public :: thermo_var_ptrs
80 ! If allocated, the following variables have nz layers.
81 real, pointer :: t(:,:,:) => null() !< Potential temperature [C ~> degC].
82 real, pointer :: s(:,:,:) => null() !< Salinity [PSU] or [gSalt/kg], generically [S ~> ppt].
83 real, pointer :: p_surf(:,:) => null() !< Ocean surface pressure used in equation of state
84 !! calculations [R L2 T-2 ~> Pa]
85 type(eos_type), pointer :: eqn_of_state => null() !< Type that indicates the
86 !! equation of state to use.
87 real :: p_ref !< The coordinate-density reference pressure [R L2 T-2 ~> Pa].
88 !! This is the pressure used to calculate Rml from
89 !! T and S when eqn_of_state is associated.
90 real :: c_p !< The heat capacity of seawater [Q C-1 ~> J degC-1 kg-1].
91 !! When conservative temperature is used, this is
92 !! constant and exactly 3991.86795711963 J degC-1 kg-1.
93 logical :: t_is_cont = .false. !< If true, the temperature variable tv%T is
94 !! actually the conservative temperature [C ~> degC].
95 logical :: s_is_abss = .false. !< If true, the salinity variable tv%S is
96 !! actually the absolute salinity in units of [S ~> gSalt kg-1].
97 real :: min_salinity !< The minimum value of salinity when BOUND_SALINITY=True [S ~> ppt].
98 real, allocatable, dimension(:,:,:) :: spv_avg
99 !< The layer averaged in situ specific volume [R-1 ~> m3 kg-1].
100 integer :: valid_spv_halo = -1 !< If positive, the valid halo size for SpV_avg, or if negative
101 !! SpV_avg is not currently set.
102
103 ! These arrays are accumulated fluxes for communication with other components.
104 real, dimension(:,:), pointer :: frazil => null()
105 !< The energy needed to heat the ocean column to the
106 !! freezing point since calculate_surface_state was
107 !! last called [Q Z R ~> J m-2].
108 logical :: frazil_was_reset !< If true, frazil has not accumulated since it was last reset.
109 real, dimension(:,:), pointer :: salt_deficit => null()
110 !< The salt needed to maintain the ocean column
111 !! at a minimum salinity of MIN_SALINITY since the last time
112 !! that calculate_surface_state was called, [S R Z ~> gSalt m-2].
113 real, dimension(:,:), pointer :: tempxpme => null()
114 !< The net inflow of water into the ocean times the
115 !! temperature at which this inflow occurs since the
116 !! last call to calculate_surface_state [C R Z ~> degC kg m-2].
117 !! This should be prescribed in the forcing fields, but
118 !! as it often is not, this is a useful heat budget diagnostic.
119 real, dimension(:,:), pointer :: internal_heat => null()
120 !< Any internal or geothermal heat sources that
121 !! have been applied to the ocean since the last call to
122 !! calculate_surface_state [C R Z ~> degC kg m-2].
123 ! The following variables are most normally not used but when they are they
124 ! will be either set by parameterizations or prognostic.
125 real, pointer :: vart(:,:,:) => null() !< SGS variance of potential temperature [C2 ~> degC2].
126 real, pointer :: vars(:,:,:) => null() !< SGS variance of salinity [S2 ~> ppt2].
127 real, pointer :: covarts(:,:,:) => null() !< SGS covariance of salinity and potential
128 !! temperature [C S ~> degC ppt].
129 type(tracer_type), pointer :: tr_t => null() !< pointer to temp in tracer registry
130 type(tracer_type), pointer :: tr_s => null() !< pointer to salinty in tracer registry
131end type thermo_var_ptrs
132
133!> Pointers to all of the prognostic variables allocated in MOM_variables.F90 and MOM.F90.
134!!
135!! It is useful for sending these variables for diagnostics, and in preparation for ensembles
136!! later on. All variables have the same names as the local (public) variables
137!! they refer to in MOM.F90.
138type, public :: ocean_internal_state
139 real, pointer, dimension(:,:,:) :: &
140 t => null(), & !< Pointer to the temperature state variable [C ~> degC]
141 s => null(), & !< Pointer to the salinity state variable [S ~> ppt] (i.e., PSU or g/kg)
142 u => null(), & !< Pointer to the zonal velocity [L T-1 ~> m s-1]
143 v => null(), & !< Pointer to the meridional velocity [L T-1 ~> m s-1]
144 h => null() !< Pointer to the layer thicknesses [H ~> m or kg m-2]
145 real, pointer, dimension(:,:,:) :: &
146 uh => null(), & !< Pointer to zonal transports [H L2 T-1 ~> m3 s-1 or kg s-1]
147 vh => null() !< Pointer to meridional transports [H L2 T-1 ~> m3 s-1 or kg s-1]
148 real, pointer, dimension(:,:,:) :: &
149 cau => null(), & !< Pointer to the zonal Coriolis and Advective acceleration [L T-2 ~> m s-2]
150 cav => null(), & !< Pointer to the meridional Coriolis and Advective acceleration [L T-2 ~> m s-2]
151 pfu => null(), & !< Pointer to the zonal Pressure force acceleration [L T-2 ~> m s-2]
152 pfv => null(), & !< Pointer to the meridional Pressure force acceleration [L T-2 ~> m s-2]
153 diffu => null(), & !< Pointer to the zonal acceleration due to lateral viscosity [L T-2 ~> m s-2]
154 diffv => null(), & !< Pointer to the meridional acceleration due to lateral viscosity [L T-2 ~> m s-2]
155 pbce => null(), & !< Pointer to the baroclinic pressure force dependency on free surface movement
156 !! [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2]
157 u_accel_bt => null(), & !< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
158 v_accel_bt => null() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
159 real, pointer, dimension(:,:,:) :: &
160 u_av => null(), & !< Pointer to zonal velocity averaged over the timestep [L T-1 ~> m s-1]
161 v_av => null(), & !< Pointer to meridional velocity averaged over the timestep [L T-1 ~> m s-1]
162 u_prev => null(), & !< Pointer to zonal velocity at the end of the last timestep [L T-1 ~> m s-1]
163 v_prev => null() !< Pointer to meridional velocity at the end of the last timestep [L T-1 ~> m s-1]
165
166!> Pointers to arrays with accelerations, which can later be used for derived diagnostics, like energy balances.
167type, public :: accel_diag_ptrs
168
169 ! Each of the following fields has nz layers.
170 real, pointer, dimension(:,:,:) :: &
171 diffu => null(), & !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
172 diffv => null(), & !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
173 cau => null(), & !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
174 cav => null(), & !< Meridional Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
175 pfu => null(), & !< Zonal acceleration due to pressure forces [L T-2 ~> m s-2]
176 pfv => null(), & !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2]
177 du_dt_visc => null(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2]
178 dv_dt_visc => null(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2]
179 du_dt_visc_gl90 => null(), &!< Zonal acceleration due to GL90 vertical viscosity
180 ! (is included in du_dt_visc) [L T-2 ~> m s-2]
181 dv_dt_visc_gl90 => null(), &!< Meridional acceleration due to GL90 vertical viscosity
182 ! (is included in dv_dt_visc) [L T-2 ~> m s-2]
183 du_dt_str => null(), & !< Zonal acceleration due to the surface stress (included
184 !! in du_dt_visc) [L T-2 ~> m s-2]
185 dv_dt_str => null(), & !< Meridional acceleration due to the surface stress (included
186 !! in dv_dt_visc) [L T-2 ~> m s-2]
187 du_dt_dia => null(), & !< Zonal acceleration due to diapycnal mixing [L T-2 ~> m s-2]
188 dv_dt_dia => null(), & !< Meridional acceleration due to diapycnal mixing [L T-2 ~> m s-2]
189 u_accel_bt => null(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
190 v_accel_bt => null(), &!< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
191
192 ! sal_[uv] and tide_[uv] are 3D fields because of their baroclinic component in Boussinesq mode.
193 sal_u => null(), & !< Zonal acceleration due to self-attraction and loading [L T-2 ~> m s-2]
194 sal_v => null(), & !< Meridional acceleration due to self-attraction and loading [L T-2 ~> m s-2]
195 tides_u => null(), & !< Zonal acceleration due to astronomical tidal forcing [L T-2 ~> m s-2]
196 tides_v => null() !< Meridional acceleration due to astronomical tidal forcing [L T-2 ~> m s-2]
197 real, pointer, dimension(:,:,:) :: du_other => null()
198 !< Zonal velocity changes due to any other processes that are
199 !! not due to any explicit accelerations [L T-1 ~> m s-1].
200 real, pointer, dimension(:,:,:) :: dv_other => null()
201 !< Meridional velocity changes due to any other processes that are
202 !! not due to any explicit accelerations [L T-1 ~> m s-1].
203
204 ! Sub-terms of [uv]_accel_bt
205 real, pointer :: bt_pgf_u(:,:,:) => null() !< Zonal acceleration due to anomalous pressure gradient from
206 !! barotropic solver, a 3D component of u_accel_bt that includes both
207 !! PFuBT and the offset term for central differencing timestepping
208 !! [L T-2 ~> m s-2]
209 real, pointer :: bt_pgf_v(:,:,:) => null() !< Meridional acceleration due to anomalous pressure gradient from
210 !! barotropic solver, a 3D component of v_accel_bt that includes both
211 !! PFvBT and the offset term for central differencing timestepping
212 !! [L T-2 ~> m s-2]
213 real, pointer :: bt_cor_u(:,:) => null() !< Zonal acceleration due to anomalous Coriolis force from barotropic
214 !! solver, a 2D component of u_accel_bt [L T-2 ~> m s-2]
215 real, pointer :: bt_cor_v(:,:) => null() !< Meridional acceleration due to anomalous Coriolis force from barotropic
216 !! solver, a 2D component of v_accel_bt [L T-2 ~> m s-2]
217 real, pointer :: bt_lwd_u(:,:) => null() !< Zonal acceleration due to linear wave drag from barotropic solver,
218 !! a 2D component of u_accel_bt [L T-2 ~> m s-2]
219 real, pointer :: bt_lwd_v(:,:) => null() !< Meridional acceleration due to linear wave drag from barotropic solver,
220 !! a 2D component of v_accel_bt [L T-2 ~> m s-2]
221
222 ! These accelerations are sub-terms included in the accelerations above.
223 real, pointer :: gradkeu(:,:,:) => null() !< gradKEu = - d/dx(u2) [L T-2 ~> m s-2]
224 real, pointer :: gradkev(:,:,:) => null() !< gradKEv = - d/dy(u2) [L T-2 ~> m s-2]
225 real, pointer :: rv_x_v(:,:,:) => null() !< rv_x_v = rv * v at u [L T-2 ~> m s-2]
226 real, pointer :: rv_x_u(:,:,:) => null() !< rv_x_u = rv * u at v [L T-2 ~> m s-2]
227
228 real, pointer :: diag_hfrac_u(:,:,:) => null() !< Fractional layer thickness at u points [nondim]
229 real, pointer :: diag_hfrac_v(:,:,:) => null() !< Fractional layer thickness at v points [nondim]
230 real, pointer :: diag_hu(:,:,:) => null() !< layer thickness at u points, modulated by the viscous
231 !! remnant and fractional open areas [H ~> m or kg m-2]
232 real, pointer :: diag_hv(:,:,:) => null() !< layer thickness at v points, modulated by the viscous
233 !! remnant and fractional open areas [H ~> m or kg m-2]
234
235 real, pointer :: visc_rem_u(:,:,:) => null() !< viscous remnant at u points [nondim]
236 real, pointer :: visc_rem_v(:,:,:) => null() !< viscous remnant at v points [nondim]
237
238end type accel_diag_ptrs
239
240!> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.
241type, public :: cont_diag_ptrs
242
243! Each of the following fields has nz layers.
244 real, pointer, dimension(:,:,:) :: &
245 uh => null(), & !< Resolved zonal layer thickness fluxes, [H L2 T-1 ~> m3 s-1 or kg s-1]
246 vh => null(), & !< Resolved meridional layer thickness fluxes, [H L2 T-1 ~> m3 s-1 or kg s-1]
247 uh_smooth => null(), & !< Interface height smoothing induced zonal volume fluxes [H L2 T-1 ~> m3 s-1 or kg s-1]
248 vh_smooth => null(), & !< Interface height smoothing induced meridional volume fluxes [H L2 T-1 ~> m3 s-1 or kg s-1]
249 uhgm => null(), & !< Isopycnal height diffusion induced zonal volume fluxes [H L2 T-1 ~> m3 s-1 or kg s-1]
250 vhgm => null() !< Isopycnal height diffusion induced meridional volume fluxes [H L2 T-1 ~> m3 s-1 or kg s-1]
251
252! Each of the following fields is found at nz+1 interfaces.
253 real, pointer :: diapyc_vel(:,:,:) => null() !< The net diapycnal velocity [H T-1 ~> m s-1 or kg m-2 s-1]
254
255end type cont_diag_ptrs
256
257!> Vertical viscosities, drag coefficients, and related fields.
258type, public :: vertvisc_type
259 real, allocatable, dimension(:,:) :: &
260 bbl_thick_u, & !< The bottom boundary layer thickness at the u-points [Z ~> m].
261 bbl_thick_v, & !< The bottom boundary layer thickness at the v-points [Z ~> m].
262 kv_bbl_u, & !< The bottom boundary layer viscosity at the u-points [H Z T-1 ~> m2 s-1 or Pa s]
263 kv_bbl_v, & !< The bottom boundary layer viscosity at the v-points [H Z T-1 ~> m2 s-1 or Pa s]
264 ustar_bbl, & !< The turbulence velocity in the bottom boundary layer at
265 !! h points [H T-1 ~> m s-1 or kg m-2 s-1].
266 bbl_meanke_loss, & !< The viscous loss of mean kinetic energy in the bottom boundary layer
267 !! [H L2 T-3 ~> m3 s-3 or W m-2].
268 bbl_meanke_loss_sqrtcd, & !< The viscous loss of mean kinetic energy in the bottom boundary layer
269 !! divided by the square root of the drag coefficient [H L2 T-3 ~> m3 s-3 or W m-2].
270 !! This is being set only to retain old answers, and should be phased out.
271 taux_shelf, & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa].
272 tauy_shelf !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa].
273 real, allocatable, dimension(:,:) :: tbl_thick_shelf_u
274 !< Thickness of the viscous top boundary layer under ice shelves at u-points [Z ~> m].
275 real, allocatable, dimension(:,:) :: tbl_thick_shelf_v
276 !< Thickness of the viscous top boundary layer under ice shelves at v-points [Z ~> m].
277 real, allocatable, dimension(:,:) :: kv_tbl_shelf_u
278 !< Viscosity in the viscous top boundary layer under ice shelves at
279 !! u-points [H Z T-1 ~> m2 s-1 or Pa s]
280 real, allocatable, dimension(:,:) :: kv_tbl_shelf_v
281 !< Viscosity in the viscous top boundary layer under ice shelves at
282 !! v-points [H Z T-1 ~> m2 s-1 or Pa s]
283 real, allocatable, dimension(:,:) :: nkml_visc_u
284 !< The number of layers in the viscous surface mixed layer at u-points [nondim].
285 !! This is not an integer because there may be fractional layers, and it is stored in
286 !! terms of layers, not depth, to facilitate the movement of the viscous boundary layer
287 !! with the flow.
288 real, allocatable, dimension(:,:) :: nkml_visc_v
289 !< The number of layers in the viscous surface mixed layer at v-points [nondim].
290 real, allocatable, dimension(:,:,:) :: &
291 ray_u, & !< The Rayleigh drag velocity to be applied to each layer at u-points [H T-1 ~> m s-1 or Pa s m-1].
292 ray_v !< The Rayleigh drag velocity to be applied to each layer at v-points [H T-1 ~> m s-1 or Pa s m-1].
293
294 ! The following elements are pointers so they can be used as targets for pointers in the restart registry.
295 real, pointer, dimension(:,:) :: mld => null()
296 !< Instantaneous active mixing layer depth as used by the mixed layer restratification
297 !! parameterization [Z ~> m].
298 real, pointer, dimension(:,:) :: lam2 => null()
299 !< (Langmuir Number)^-2 [nondim].
300 real, pointer, dimension(:,:) :: h_ml => null()
301 !< Instantaneous active mixing layer thickness as used by the mixed layer restratification
302 !! parameterization [H ~> m or kg m-2]
303 real, pointer, dimension(:,:) :: mld_param => null()
304 !< Instantaneous active mixed or mixing layer depth as used by the brine plume parameterization.
305 !! It could be coordinated with MLD above but we may want the ability to use different scales
306 !! in different parameterizations [Z ~> m].
307 real, pointer, dimension(:,:) :: h_ml_param => null()
308 !< Instantaneous active mixed or mixing layer thickness as used by the brine plume
309 !! parameterization [H ~> m or kg m-2].
310 real, pointer, dimension(:,:) :: sfc_buoy_flx => null()
311 !< Surface buoyancy flux (derived) [Z2 T-3 ~> m2 s-3].
312 real, pointer, dimension(:,:,:) :: kd_shear => null()
313 !< The shear-driven turbulent diapycnal diffusivity at the interfaces between layers
314 !! in tracer columns [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
315 real, pointer, dimension(:,:,:) :: kv_shear => null()
316 !< The shear-driven turbulent vertical viscosity at the interfaces between layers
317 !! in tracer columns [H Z T-1 ~> m2 s-1 or Pa s]
318 real, pointer, dimension(:,:,:) :: kv_shear_bu => null()
319 !< The shear-driven turbulent vertical viscosity at the interfaces between layers in
320 !! corner columns [H Z T-1 ~> m2 s-1 or Pa s]
321 real, pointer, dimension(:,:,:) :: kv_slow => null()
322 !< The turbulent vertical viscosity component due to "slow" processes (e.g., tidal,
323 !! background, convection etc) [H Z T-1 ~> m2 s-1 or Pa s]
324 real, pointer, dimension(:,:,:) :: tke_turb => null()
325 !< The turbulent kinetic energy per unit mass at the interfaces [Z2 T-2 ~> m2 s-2].
326 !! This may be at the tracer or corner points
327end type vertvisc_type
328
329!> Container for information about the summed layer transports
330!! and how they will vary as the barotropic velocity is changed.
331type, public :: bt_cont_type
332 real, allocatable :: fa_u_ee(:,:) !< The effective open face area for zonal barotropic transport
333 !! drawing from locations far to the east [H L ~> m2 or kg m-1].
334 real, allocatable :: fa_u_e0(:,:) !< The effective open face area for zonal barotropic transport
335 !! drawing from nearby to the east [H L ~> m2 or kg m-1].
336 real, allocatable :: fa_u_w0(:,:) !< The effective open face area for zonal barotropic transport
337 !! drawing from nearby to the west [H L ~> m2 or kg m-1].
338 real, allocatable :: fa_u_ww(:,:) !< The effective open face area for zonal barotropic transport
339 !! drawing from locations far to the west [H L ~> m2 or kg m-1].
340 real, allocatable :: ubt_ww(:,:) !< uBT_WW is the barotropic velocity [L T-1 ~> m s-1], beyond which the
341 !! marginal open face area is FA_u_WW. uBT_WW must be non-negative.
342 real, allocatable :: ubt_ee(:,:) !< uBT_EE is a barotropic velocity [L T-1 ~> m s-1], beyond which the
343 !! marginal open face area is FA_u_EE. uBT_EE must be non-positive.
344 real, allocatable :: fa_v_nn(:,:) !< The effective open face area for meridional barotropic transport
345 !! drawing from locations far to the north [H L ~> m2 or kg m-1].
346 real, allocatable :: fa_v_n0(:,:) !< The effective open face area for meridional barotropic transport
347 !! drawing from nearby to the north [H L ~> m2 or kg m-1].
348 real, allocatable :: fa_v_s0(:,:) !< The effective open face area for meridional barotropic transport
349 !! drawing from nearby to the south [H L ~> m2 or kg m-1].
350 real, allocatable :: fa_v_ss(:,:) !< The effective open face area for meridional barotropic transport
351 !! drawing from locations far to the south [H L ~> m2 or kg m-1].
352 real, allocatable :: vbt_ss(:,:) !< vBT_SS is the barotropic velocity, [L T-1 ~> m s-1], beyond which the
353 !! marginal open face area is FA_v_SS. vBT_SS must be non-negative.
354 real, allocatable :: vbt_nn(:,:) !< vBT_NN is the barotropic velocity, [L T-1 ~> m s-1], beyond which the
355 !! marginal open face area is FA_v_NN. vBT_NN must be non-positive.
356 real, allocatable :: h_u(:,:,:) !< An effective thickness at zonal faces, taking into account the effects
357 !! of vertical viscosity and fractional open areas [H ~> m or kg m-2].
358 !! This is primarily used as a non-normalized weight in determining
359 !! the depth averaged accelerations for the barotropic solver.
360 real, allocatable :: h_v(:,:,:) !< An effective thickness at meridional faces, taking into account the effects
361 !! of vertical viscosity and fractional open areas [H ~> m or kg m-2].
362 !! This is primarily used as a non-normalized weight in determining
363 !! the depth averaged accelerations for the barotropic solver.
364 type(group_pass_type) :: pass_polarity_bt !< Structure for polarity group halo updates
365 type(group_pass_type) :: pass_fa_uv !< Structure for face area group halo updates
366end type bt_cont_type
367
368!> Container for grids modifying cell metric at porous barriers
369type, public :: porous_barrier_type
370 ! Each of the following fields has nz layers.
371 real, allocatable :: por_face_areau(:,:,:) !< fractional open area of U-faces [nondim]
372 real, allocatable :: por_face_areav(:,:,:) !< fractional open area of V-faces [nondim]
373 ! Each of the following fields is found at nz+1 interfaces.
374 real, allocatable :: por_layer_widthu(:,:,:) !< fractional open width of U-faces [nondim]
375 real, allocatable :: por_layer_widthv(:,:,:) !< fractional open width of V-faces [nondim]
376end type porous_barrier_type
377
378contains
379
380!> Allocates the fields for the surface (return) properties of
381!! the ocean model. Unused fields are unallocated.
382subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
383 gas_fields_ocn, use_meltpot, use_iceshelves, &
384 omit_frazil, sfc_state_in, turns, use_MARBL_tracers)
385 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
386 type(surface), intent(inout) :: sfc_state !< ocean surface state type to be allocated.
387 logical, optional, intent(in) :: use_temperature !< If true, allocate the space for thermodynamic variables.
388 logical, optional, intent(in) :: do_integrals !< If true, allocate the space for vertically
389 !! integrated fields.
390 type(coupler_1d_bc_type), &
391 optional, intent(in) :: gas_fields_ocn !< If present, this type describes the
392 !! ocean and surface-ice fields that will participate
393 !! in the calculation of additional gas or other
394 !! tracer fluxes, and can be used to spawn related
395 !! internal variables in the ice model.
396 logical, optional, intent(in) :: use_meltpot !< If true, allocate the space for melt potential
397 logical, optional, intent(in) :: use_iceshelves !< If true, allocate the space for the stresses
398 !! under ice shelves.
399 logical, optional, intent(in) :: omit_frazil !< If present and false, do not allocate the space to
400 !! pass frazil fluxes to the coupler
401 type(surface), &
402 optional, intent(in) :: sfc_state_in !< If present and its tr_fields are initialized,
403 !! this type describes the ocean and surface-ice fields that
404 !! will participate in the calculation of additional gas or
405 !! other tracer fluxes, and can be used to spawn related
406 !! internal variables in the ice model. If gas_fields_ocn
407 !! is present, it is used and tr_fields_in is ignored.
408 integer, optional, intent(in) :: turns !< If present, the number of counterclockwise quarter
409 !! turns to use on the new grid.
410 logical, optional, intent(in) :: use_marbl_tracers !< If true, allocate the space for CO2 flux from MARBL
411
412 ! local variables
413 logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil, alloc_fco2
414 logical :: even_turns ! True if turns is absent or even
415 integer :: tr_field_i_mem(4), tr_field_j_mem(4)
416 integer :: is, ie, js, je, isd, ied, jsd, jed
417 integer :: isdb, iedb, jsdb, jedb
418
419 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
420 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
421 isdb = g%isdB ; iedb = g%iedB ; jsdb = g%jsdB ; jedb = g%jedB
422
423 use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature
424 alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals
425 use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot
426 alloc_iceshelves = .false. ; if (present(use_iceshelves)) alloc_iceshelves = use_iceshelves
427 alloc_frazil = .true. ; if (present(omit_frazil)) alloc_frazil = .not.omit_frazil
428 alloc_fco2 = .false. ; if (present(use_marbl_tracers)) alloc_fco2 = use_marbl_tracers
429
430 if (sfc_state%arrays_allocated) return
431
432 if (use_temp) then
433 allocate(sfc_state%SST(isd:ied,jsd:jed), source=0.0)
434 allocate(sfc_state%SSS(isd:ied,jsd:jed), source=0.0)
435 else
436 allocate(sfc_state%sfc_density(isd:ied,jsd:jed), source=0.0)
437 endif
438 if (use_temp .and. alloc_frazil) then
439 allocate(sfc_state%frazil(isd:ied,jsd:jed), source=0.0)
440 endif
441 allocate(sfc_state%sea_lev(isd:ied,jsd:jed), source=0.0)
442 allocate(sfc_state%Hml(isd:ied,jsd:jed), source=0.0)
443 allocate(sfc_state%u(isdb:iedb,jsd:jed), source=0.0)
444 allocate(sfc_state%v(isd:ied,jsdb:jedb), source=0.0)
445
446 if (use_melt_potential) then
447 allocate(sfc_state%melt_potential(isd:ied,jsd:jed), source=0.0)
448 endif
449
450 if (alloc_integ) then
451 ! Allocate structures for the vertically integrated ocean_mass, ocean_heat, and ocean_salt.
452 allocate(sfc_state%ocean_mass(isd:ied,jsd:jed), source=0.0)
453 if (use_temp) then
454 allocate(sfc_state%ocean_heat(isd:ied,jsd:jed), source=0.0)
455 allocate(sfc_state%ocean_salt(isd:ied,jsd:jed), source=0.0)
456 endif
457 endif
458
459 if (alloc_iceshelves) then
460 allocate(sfc_state%taux_shelf(isdb:iedb,jsd:jed), source=0.0)
461 allocate(sfc_state%tauy_shelf(isd:ied,jsdb:jedb), source=0.0)
462 endif
463
464 ! The data fields in the coupler_2d_bc_type are never rotated.
465 even_turns = .true. ; if (present(turns)) even_turns = (modulo(turns, 2) == 0)
466 if (even_turns) then
467 tr_field_i_mem(1:4) = (/is,is,ie,ie/) ; tr_field_j_mem(1:4) = (/js,js,je,je/)
468 else
469 tr_field_i_mem(1:4) = (/js,js,je,je/) ; tr_field_j_mem(1:4) = (/is,is,ie,ie/)
470 endif
471 if (present(gas_fields_ocn)) then
472 call coupler_type_spawn(gas_fields_ocn, sfc_state%tr_fields, &
473 tr_field_i_mem, tr_field_j_mem, as_needed=.true.)
474 elseif (present(sfc_state_in)) then
475 if (coupler_type_initialized(sfc_state_in%tr_fields)) then
476 call coupler_type_spawn(sfc_state_in%tr_fields, sfc_state%tr_fields, &
477 tr_field_i_mem, tr_field_j_mem, as_needed=.true.)
478 endif
479 endif
480
481 if (alloc_fco2) then
482 allocate(sfc_state%fco2(isd:ied,jsd:jed), source=0.0)
483 endif
484
485 sfc_state%arrays_allocated = .true.
486
487end subroutine allocate_surface_state
488
489!> Deallocates the elements of a surface state type.
490subroutine deallocate_surface_state(sfc_state)
491 type(surface), intent(inout) :: sfc_state !< ocean surface state type to be deallocated here.
492
493 if (.not.sfc_state%arrays_allocated) return
494
495 if (allocated(sfc_state%melt_potential)) deallocate(sfc_state%melt_potential)
496 if (allocated(sfc_state%SST)) deallocate(sfc_state%SST)
497 if (allocated(sfc_state%SSS)) deallocate(sfc_state%SSS)
498 if (allocated(sfc_state%sfc_density)) deallocate(sfc_state%sfc_density)
499 if (allocated(sfc_state%sea_lev)) deallocate(sfc_state%sea_lev)
500 if (allocated(sfc_state%Hml)) deallocate(sfc_state%Hml)
501 if (allocated(sfc_state%u)) deallocate(sfc_state%u)
502 if (allocated(sfc_state%v)) deallocate(sfc_state%v)
503 if (allocated(sfc_state%ocean_mass)) deallocate(sfc_state%ocean_mass)
504 if (allocated(sfc_state%ocean_heat)) deallocate(sfc_state%ocean_heat)
505 if (allocated(sfc_state%ocean_salt)) deallocate(sfc_state%ocean_salt)
506 if (allocated(sfc_state%fco2)) deallocate(sfc_state%fco2)
507 call coupler_type_destructor(sfc_state%tr_fields)
508
509 sfc_state%arrays_allocated = .false.
510
511end subroutine deallocate_surface_state
512
513!> Rotate the surface state fields from the input to the model indices.
514subroutine rotate_surface_state(sfc_state_in, sfc_state, G, turns)
515 type(surface), intent(in) :: sfc_state_in !< The input unrotated surface state type that is the data source.
516 type(surface), intent(inout) :: sfc_state !< The rotated surface state type whose arrays will be filled in
517 type(ocean_grid_type), intent(in) :: g !< The ocean grid structure
518 integer, intent(in) :: turns !< The number of counterclockwise quarter turns to use on the rotated grid.
519
520 logical :: use_temperature, do_integrals, use_melt_potential, use_iceshelves
521
522 ! NOTE: Many of these are weak tests, since only one is checked
523 use_temperature = allocated(sfc_state_in%SST) &
524 .and. allocated(sfc_state_in%SSS)
525 use_melt_potential = allocated(sfc_state_in%melt_potential)
526 do_integrals = allocated(sfc_state_in%ocean_mass)
527 use_iceshelves = allocated(sfc_state_in%taux_shelf) &
528 .and. allocated(sfc_state_in%tauy_shelf)
529
530 if (.not. sfc_state%arrays_allocated) then
531 call allocate_surface_state(sfc_state, g, use_temperature=use_temperature, &
532 do_integrals=do_integrals, use_meltpot=use_melt_potential, &
533 use_iceshelves=use_iceshelves, sfc_state_in=sfc_state_in, turns=turns)
534 endif
535
536 if (use_temperature) then
537 call rotate_array(sfc_state_in%SST, turns, sfc_state%SST)
538 call rotate_array(sfc_state_in%SSS, turns, sfc_state%SSS)
539 else
540 call rotate_array(sfc_state_in%sfc_density, turns, sfc_state%sfc_density)
541 endif
542
543 call rotate_array(sfc_state_in%Hml, turns, sfc_state%Hml)
544 call rotate_vector(sfc_state_in%u, sfc_state_in%v, turns, &
545 sfc_state%u, sfc_state%v)
546 call rotate_array(sfc_state_in%sea_lev, turns, sfc_state%sea_lev)
547
548 if (use_melt_potential) then
549 call rotate_array(sfc_state_in%melt_potential, turns, sfc_state%melt_potential)
550 endif
551
552 if (do_integrals) then
553 call rotate_array(sfc_state_in%ocean_mass, turns, sfc_state%ocean_mass)
554 if (use_temperature) then
555 call rotate_array(sfc_state_in%ocean_heat, turns, sfc_state%ocean_heat)
556 call rotate_array(sfc_state_in%ocean_salt, turns, sfc_state%ocean_salt)
557 call rotate_array(sfc_state_in%SSS, turns, sfc_state%SSS)
558 endif
559 endif
560
561 if (use_iceshelves) then
562 call rotate_vector(sfc_state_in%taux_shelf, sfc_state_in%tauy_shelf, turns, &
563 sfc_state%taux_shelf, sfc_state%tauy_shelf)
564 endif
565
566 if (use_temperature .and. allocated(sfc_state_in%frazil)) &
567 call rotate_array(sfc_state_in%frazil, turns, sfc_state%frazil)
568
569 ! Scalar transfers
570 sfc_state%T_is_conT = sfc_state_in%T_is_conT
571 sfc_state%S_is_absS = sfc_state_in%S_is_absS
572
573 ! NOTE: Tracer fields are handled by FMS, so are left unrotated. Any
574 ! reads/writes to tr_fields must be appropriately rotated.
575 if (coupler_type_initialized(sfc_state_in%tr_fields)) then
576 call coupler_type_copy_data(sfc_state_in%tr_fields, sfc_state%tr_fields)
577 endif
578end subroutine rotate_surface_state
579
580!> Allocates the arrays contained within a BT_cont_type and initializes them to 0.
581subroutine alloc_bt_cont_type(BT_cont, G, GV, alloc_faces)
582 type(bt_cont_type), pointer :: bt_cont !< The BT_cont_type whose elements will be allocated
583 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
584 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
585 logical, optional, intent(in) :: alloc_faces !< If present and true, allocate
586 !! memory for effective face thicknesses.
587
588 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
589 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
590 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
591
592 if (associated(bt_cont)) call mom_error(fatal, &
593 "alloc_BT_cont_type called with an associated BT_cont_type pointer.")
594
595 allocate(bt_cont)
596 allocate(bt_cont%FA_u_WW(isdb:iedb,jsd:jed), source=0.0)
597 allocate(bt_cont%FA_u_W0(isdb:iedb,jsd:jed), source=0.0)
598 allocate(bt_cont%FA_u_E0(isdb:iedb,jsd:jed), source=0.0)
599 allocate(bt_cont%FA_u_EE(isdb:iedb,jsd:jed), source=0.0)
600 allocate(bt_cont%uBT_WW(isdb:iedb,jsd:jed), source=0.0)
601 allocate(bt_cont%uBT_EE(isdb:iedb,jsd:jed), source=0.0)
602
603 allocate(bt_cont%FA_v_SS(isd:ied,jsdb:jedb), source=0.0)
604 allocate(bt_cont%FA_v_S0(isd:ied,jsdb:jedb), source=0.0)
605 allocate(bt_cont%FA_v_N0(isd:ied,jsdb:jedb), source=0.0)
606 allocate(bt_cont%FA_v_NN(isd:ied,jsdb:jedb), source=0.0)
607 allocate(bt_cont%vBT_SS(isd:ied,jsdb:jedb), source=0.0)
608 allocate(bt_cont%vBT_NN(isd:ied,jsdb:jedb), source=0.0)
609
610 if (present(alloc_faces)) then ; if (alloc_faces) then
611 allocate(bt_cont%h_u(isdb:iedb,jsd:jed,1:nz), source=0.0)
612 allocate(bt_cont%h_v(isd:ied,jsdb:jedb,1:nz), source=0.0)
613 endif ; endif
614
615end subroutine alloc_bt_cont_type
616
617!> Deallocates the arrays contained within a BT_cont_type.
618subroutine dealloc_bt_cont_type(BT_cont)
619 type(bt_cont_type), pointer :: bt_cont !< The BT_cont_type whose elements will be deallocated.
620
621 if (.not.associated(bt_cont)) return
622
623 deallocate(bt_cont%FA_u_WW) ; deallocate(bt_cont%FA_u_W0)
624 deallocate(bt_cont%FA_u_E0) ; deallocate(bt_cont%FA_u_EE)
625 deallocate(bt_cont%uBT_WW) ; deallocate(bt_cont%uBT_EE)
626
627 deallocate(bt_cont%FA_v_SS) ; deallocate(bt_cont%FA_v_S0)
628 deallocate(bt_cont%FA_v_N0) ; deallocate(bt_cont%FA_v_NN)
629 deallocate(bt_cont%vBT_SS) ; deallocate(bt_cont%vBT_NN)
630
631 if (allocated(bt_cont%h_u)) deallocate(bt_cont%h_u)
632 if (allocated(bt_cont%h_v)) deallocate(bt_cont%h_v)
633
634 deallocate(bt_cont)
635
636end subroutine dealloc_bt_cont_type
637
638!> Diagnostic checksums on various elements of a thermo_var_ptrs type for debugging.
639subroutine mom_thermovar_chksum(mesg, tv, G, US)
640 character(len=*), intent(in) :: mesg !< A message that appears in the checksum lines
641 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
642 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
643 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
644
645 ! Note that for the chksum calls to be useful for reproducing across PE
646 ! counts, there must be no redundant points, so all variables use is..ie
647 ! and js...je as their extent.
648 if (associated(tv%T)) &
649 call hchksum(tv%T, mesg//" tv%T", g%HI, unscale=us%C_to_degC)
650 if (associated(tv%S)) &
651 call hchksum(tv%S, mesg//" tv%S", g%HI, unscale=us%S_to_ppt)
652 if (associated(tv%frazil)) &
653 call hchksum(tv%frazil, mesg//" tv%frazil", g%HI, unscale=us%Q_to_J_kg*us%RZ_to_kg_m2)
654 if (associated(tv%salt_deficit)) &
655 call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", g%HI, unscale=us%RZ_to_kg_m2*us%S_to_ppt)
656 if (associated(tv%TempxPmE)) &
657 call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", g%HI, unscale=us%RZ_to_kg_m2*us%C_to_degC)
658end subroutine mom_thermovar_chksum
659
660end module mom_variables