portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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 | |
| 6 | module MOM_verticalGrid | |
| 7 | ||
| 8 | use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL | |
| 9 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 10 | use MOM_unit_scaling, only : unit_scale_type | |
| 11 | ||
| 12 | implicit none ; private | |
| 13 | ||
| 14 | #include <MOM_memory.h> | |
| 15 | ||
| 16 | public verticalGridInit, verticalGridEnd | |
| 17 | public setVerticalGridAxes | |
| 18 | public get_flux_units, get_thickness_units, get_tr_flux_units | |
| 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 | |
| 26 | type, 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 | ||
| 101 | end type verticalGrid_type | |
| 102 | ||
| 103 | contains | |
| 104 | ||
| 105 | !> Allocates and initializes the ocean model vertical grid structure. | |
| 106 | 1 | subroutine 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 | 1 | if (associated(GV)) call MOM_error(FATAL, & |
| 123 | 0 | 'verticalGridInit: called with an associated GV pointer.') |
| 124 | 1 | 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 | 1 | 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 | 1 | 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 | 1 | units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) |
| 139 | call get_param(param_file, mdl, "BOUSSINESQ", GV%Boussinesq, & | |
| 140 | 1 | "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 | 1 | default=.true., do_not_log=GV%Boussinesq) |
| 146 | 1 | 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 | 1 | 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 | 1 | 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 | 1 | units="nondim", default=0, debuggingParam=.true.) |
| 163 | 1 | if (abs(H_power) > 300) call MOM_error(FATAL, "verticalGridInit: "//& |
| 164 | 0 | "H_RESCALE_POWER is outside of the valid range of -300 to 300.") |
| 165 | 1 | H_rescale_factor = 1.0 |
| 166 | 1 | if (H_power /= 0) H_rescale_factor = 2.0**H_power |
| 167 | 1 | 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 | 0 | default=1.0) |
| 172 | 0 | 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 | 1 | "units of thickness into m.", units="m H-1", default=1.0) |
| 177 | 1 | 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 | 1 | 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 | 1 | "The number of model layers.", units="nondim", fail_if_missing=.true.) |
| 192 | #endif | |
| 193 | 1 | GV%ke = nk |
| 194 | ||
| 195 | 1 | if (GV%Boussinesq) then |
| 196 | 1 | GV%H_to_kg_m2 = US%R_to_kg_m3*GV%Rho0 * GV%H_to_m |
| 197 | 1 | GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2 |
| 198 | 1 | GV%m_to_H = 1.0 / GV%H_to_m |
| 199 | 1 | GV%H_to_MKS = GV%H_to_m |
| 200 | 1 | GV%m2_s_to_HZ_T = GV%m_to_H * US%m_to_Z * US%T_to_s |
| 201 | ||
| 202 | 1 | GV%H_to_Z = GV%H_to_m * US%m_to_Z |
| 203 | 1 | GV%Z_to_H = US%Z_to_m * GV%m_to_H |
| 204 | else | |
| 205 | 0 | 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 | 0 | GV%m_to_H = US%R_to_kg_m3*rho_Kv * GV%kg_m2_to_H |
| 208 | 0 | GV%H_to_MKS = GV%H_to_kg_m2 |
| 209 | 0 | 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 | 0 | GV%H_to_m = 1.0 / GV%m_to_H |
| 211 | ||
| 212 | 0 | GV%H_to_Z = US%m_to_Z * ( GV%H_to_kg_m2 / (US%R_to_kg_m3*GV%Rho0) ) |
| 213 | 0 | GV%Z_to_H = US%Z_to_m * ( US%R_to_kg_m3*GV%Rho0 * GV%kg_m2_to_H ) |
| 214 | endif | |
| 215 | ||
| 216 | 1 | GV%Angstrom_H = (US%Z_to_m * GV%m_to_H) * GV%Angstrom_Z |
| 217 | 1 | GV%Angstrom_m = US%Z_to_m * GV%Angstrom_Z |
| 218 | ||
| 219 | 1 | GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H, GV%m_to_H*1e-17) |
| 220 | 1 | GV%dZ_subroundoff = 1e-20 * max(GV%Angstrom_Z, US%m_to_Z*1e-17) |
| 221 | ||
| 222 | 1 | 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 | 1 | GV%H_to_RZ = GV%H_to_kg_m2 * US%kg_m3_to_R * US%m_to_Z |
| 225 | 1 | GV%RZ_to_H = GV%kg_m2_to_H * US%R_to_kg_m3 * US%Z_to_m |
| 226 | ||
| 227 | 1 | GV%HZ_T_to_m2_s = 1.0 / GV%m2_s_to_HZ_T |
| 228 | 1 | 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 | 1 | call log_param(param_file, mdl, "M to THICKNESS", GV%m_to_H*H_rescale_factor, units="H m-1") |
| 239 | 1 | call log_param(param_file, mdl, "M to THICKNESS rescaled by 2^-n", GV%m_to_H, units="2^n H m-1") |
| 240 | 1 | 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 | 1 | allocate( GV%sInterface(nk+1) ) |
| 243 | 1 | allocate( GV%sLayer(nk) ) |
| 244 | 77 | allocate( GV%g_prime(nk+1), source=0.0 ) |
| 245 | 76 | allocate( GV%Rlay(nk), source=0.0 ) |
| 246 | ||
| 247 | 1 | end subroutine verticalGridInit |
| 248 | ||
| 249 | !> Returns the model's thickness units, usually m or kg/m^2. | |
| 250 | 12 | function 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 | 12 | if (GV%Boussinesq) then |
| 258 | 12 | get_thickness_units = "m" |
| 259 | else | |
| 260 | 0 | get_thickness_units = "kg m-2" |
| 261 | endif | |
| 262 | 12 | end function get_thickness_units |
| 263 | ||
| 264 | !> Returns the model's thickness flux units, usually m^3/s or kg/s. | |
| 265 | 4 | function 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 | 4 | if (GV%Boussinesq) then |
| 273 | 4 | get_flux_units = "m3 s-1" |
| 274 | else | |
| 275 | 0 | get_flux_units = "kg s-1" |
| 276 | endif | |
| 277 | 4 | end function get_flux_units |
| 278 | ||
| 279 | !> Returns the model's tracer flux units. | |
| 280 | 1 | function 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 | 1 | cnt = 0 |
| 302 | 1 | if (present(tr_units)) cnt = cnt+1 |
| 303 | 1 | if (present(tr_vol_conc_units)) cnt = cnt+1 |
| 304 | 1 | if (present(tr_mass_conc_units)) cnt = cnt+1 |
| 305 | ||
| 306 | 1 | 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 | 0 | "must be present.") |
| 309 | 1 | if (cnt > 1) call MOM_error(FATAL, "get_tr_flux_units: Only one of "//& |
| 310 | 0 | "tr_units, tr_vol_conc_units, and tr_mass_conc_units may be present.") |
| 311 | 1 | if (present(tr_units)) then |
| 312 | 1 | if (GV%Boussinesq) then |
| 313 | 1 | get_tr_flux_units = trim(tr_units)//" m3 s-1" |
| 314 | else | |
| 315 | 0 | get_tr_flux_units = trim(tr_units)//" kg s-1" |
| 316 | endif | |
| 317 | endif | |
| 318 | 1 | if (present(tr_vol_conc_units)) then |
| 319 | 0 | if (GV%Boussinesq) then |
| 320 | 0 | get_tr_flux_units = trim(tr_vol_conc_units)//" s-1" |
| 321 | else | |
| 322 | 0 | get_tr_flux_units = trim(tr_vol_conc_units)//" m-3 kg s-1" |
| 323 | endif | |
| 324 | endif | |
| 325 | 1 | if (present(tr_mass_conc_units)) then |
| 326 | 0 | if (GV%Boussinesq) then |
| 327 | 0 | get_tr_flux_units = trim(tr_mass_conc_units)//" kg-1 m3 s-1" |
| 328 | else | |
| 329 | 0 | get_tr_flux_units = trim(tr_mass_conc_units)//" s-1" |
| 330 | endif | |
| 331 | endif | |
| 332 | ||
| 333 | 1 | end function get_tr_flux_units |
| 334 | ||
| 335 | !> This sets the coordinate data for the "layer mode" of the isopycnal model. | |
| 336 | 1 | subroutine 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 | 1 | nk = GV%ke |
| 346 | ||
| 347 | 1 | GV%zAxisLongName = 'Target Potential Density' |
| 348 | 1 | GV%zAxisUnits = 'kg m-3' |
| 349 | 76 | do k=1,nk ; GV%sLayer(k) = scale*Rlay(k) ; enddo |
| 350 | 1 | if (nk > 1) then |
| 351 | 1 | GV%sInterface(1) = scale * (1.5*Rlay(1) - 0.5*Rlay(2)) |
| 352 | 75 | do K=2,nk ; GV%sInterface(K) = scale * 0.5*( Rlay(k-1) + Rlay(k) ) ; enddo |
| 353 | 1 | GV%sInterface(nk+1) = scale * (1.5*Rlay(nk) - 0.5*Rlay(nk-1)) |
| 354 | else | |
| 355 | 0 | GV%sInterface(1) = 0.0 ; GV%sInterface(nk+1) = 2.0*scale*Rlay(nk) |
| 356 | endif | |
| 357 | ||
| 358 | 1 | end subroutine setVerticalGridAxes |
| 359 | ||
| 360 | !> Deallocates the model's vertical grid structure. | |
| 361 | 1 | subroutine verticalGridEnd( GV ) |
| 362 | type(verticalGrid_type), pointer :: GV !< The ocean's vertical grid structure | |
| 363 | ||
| 364 | 1 | deallocate( GV%g_prime, GV%Rlay ) |
| 365 | 1 | deallocate( GV%sInterface , GV%sLayer ) |
| 366 | 1 | deallocate( GV ) |
| 367 | ||
| 368 | 1 | end subroutine verticalGridEnd |
| 369 | ||
| 370 | 0 | end module MOM_verticalGrid |