MOM6
MOM_verticalGrid.F90
Go to the documentation of this file.
1 !> Provides a transparent vertical ocean grid type and supporting routines
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_error_handler, only : mom_error, mom_mesg, fatal
9 
10 implicit none ; private
11 
12 #include <MOM_memory.h>
13 
17 
18 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
19 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
20 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
21 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
22 
23 !> Describes the vertical ocean grid, including unit conversion factors
24 type, public :: verticalgrid_type
25 
26  ! Commonly used parameters
27  integer :: ke !< The number of layers/levels in the vertical
28  real :: max_depth !< The maximum depth of the ocean [Z ~> m].
29  real :: mks_g_earth !< The gravitational acceleration in unscaled MKS units [m s-2].
30  real :: g_earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
31  real :: rho0 !< The density used in the Boussinesq approximation or nominal
32  !! density used to convert depths into mass units [kg m-3].
33 
34  ! Vertical coordinate descriptions for diagnostics and I/O
35  character(len=40) :: zaxisunits !< The units that vertical coordinates are written in
36  character(len=40) :: zaxislongname !< Coordinate name to appear in files,
37  !! e.g. "Target Potential Density" or "Height"
38  real, allocatable, dimension(:) :: slayer !< Coordinate values of layer centers
39  real, allocatable, dimension(:) :: sinterface !< Coordinate values on interfaces
40  integer :: direction = 1 !< Direction defaults to 1, positive up.
41 
42  ! The following variables give information about the vertical grid.
43  logical :: boussinesq !< If true, make the Boussinesq approximation.
44  real :: angstrom_h !< A one-Angstrom thickness in the model thickness units [H ~> m or kg m-2].
45  real :: angstrom_z !< A one-Angstrom thickness in the model depth units [Z ~> m].
46  real :: angstrom_m !< A one-Angstrom thickness [m].
47  real :: h_subroundoff !< A thickness that is so small that it can be added to a thickness of
48  !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2].
49  !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m.
50  real, allocatable, dimension(:) :: &
51  g_prime, & !< The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2].
52  rlay !< The target coordinate value (potential density) in each layer [R ~> kg m-3].
53  integer :: nkml = 0 !< The number of layers at the top that should be treated
54  !! as parts of a homogeneous region.
55  integer :: nk_rho_varies = 0 !< The number of layers at the top where the
56  !! density does not track any target density.
57  real :: h_to_kg_m2 !< A constant that translates thicknesses from the units of thickness to kg m-2.
58  real :: kg_m2_to_h !< A constant that translates thicknesses from kg m-2 to the units of thickness.
59  real :: m_to_h !< A constant that translates distances in m to the units of thickness.
60  real :: h_to_m !< A constant that translates distances in the units of thickness to m.
61  real :: h_to_pa !< A constant that translates the units of thickness to pressure [Pa].
62  real :: h_to_z !< A constant that translates thickness units to the units of depth.
63  real :: z_to_h !< A constant that translates depth units to thickness units.
64  real :: h_to_rz !< A constant that translates thickness units to the units of mass per unit area.
65  real :: rz_to_h !< A constant that translates mass per unit area units to thickness units.
66  real :: h_to_mks !< A constant that translates thickness units to its
67  !! MKS unit (m or kg m-2) based on GV%Boussinesq
68 
69  real :: m_to_h_restart = 0.0 !< A copy of the m_to_H that is used in restart files.
70 end type verticalgrid_type
71 
72 contains
73 
74 !> Allocates and initializes the ocean model vertical grid structure.
75 subroutine verticalgridinit( param_file, GV, US )
76  type(param_file_type), intent(in) :: param_file !< Parameter file handle/type
77  type(verticalgrid_type), pointer :: gv !< The container for vertical grid data
78  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
79  ! This routine initializes the verticalGrid_type structure (GV).
80  ! All memory is allocated but not necessarily set to meaningful values until later.
81 
82  ! Local variables
83  integer :: nk, h_power
84  real :: h_rescale_factor
85  ! This include declares and sets the variable "version".
86 # include "version_variable.h"
87  character(len=16) :: mdl = 'MOM_verticalGrid'
88 
89  if (associated(gv)) call mom_error(fatal, &
90  'verticalGridInit: called with an associated GV pointer.')
91  allocate(gv)
92 
93  ! Read all relevant parameters and write them to the model log.
94  call log_version(param_file, mdl, version, &
95  "Parameters providing information about the vertical grid.")
96  call get_param(param_file, mdl, "G_EARTH", gv%mks_g_Earth, &
97  "The gravitational acceleration of the Earth.", &
98  units="m s-2", default = 9.80)
99  call get_param(param_file, mdl, "RHO_0", gv%Rho0, &
100  "The mean ocean density used with BOUSSINESQ true to "//&
101  "calculate accelerations and the mass for conservation "//&
102  "properties, or with BOUSSINSEQ false to convert some "//&
103  "parameters from vertical units of m to kg m-2.", &
104  units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
105  call get_param(param_file, mdl, "BOUSSINESQ", gv%Boussinesq, &
106  "If true, make the Boussinesq approximation.", default=.true.)
107  call get_param(param_file, mdl, "ANGSTROM", gv%Angstrom_m, &
108  "The minimum layer thickness, usually one-Angstrom.", &
109  units="m", default=1.0e-10)
110  call get_param(param_file, mdl, "H_RESCALE_POWER", h_power, &
111  "An integer power of 2 that is used to rescale the model's "//&
112  "intenal units of thickness. Valid values range from -300 to 300.", &
113  units="nondim", default=0, debuggingparam=.true.)
114  if (abs(h_power) > 300) call mom_error(fatal, "verticalGridInit: "//&
115  "H_RESCALE_POWER is outside of the valid range of -300 to 300.")
116  h_rescale_factor = 1.0
117  if (h_power /= 0) h_rescale_factor = 2.0**h_power
118  if (.not.gv%Boussinesq) then
119  call get_param(param_file, mdl, "H_TO_KG_M2", gv%H_to_kg_m2,&
120  "A constant that translates thicknesses from the model's "//&
121  "internal units of thickness to kg m-2.", units="kg m-2 H-1", &
122  default=1.0)
123  gv%H_to_kg_m2 = gv%H_to_kg_m2 * h_rescale_factor
124  else
125  call get_param(param_file, mdl, "H_TO_M", gv%H_to_m, &
126  "A constant that translates the model's internal "//&
127  "units of thickness into m.", units="m H-1", default=1.0)
128  gv%H_to_m = gv%H_to_m * h_rescale_factor
129  endif
130  gv%g_Earth = us%m_to_L**2*us%Z_to_m*us%T_to_s**2 * gv%mks_g_Earth
131 #ifdef STATIC_MEMORY_
132  ! Here NK_ is a macro, while nk is a variable.
133  call get_param(param_file, mdl, "NK", nk, &
134  "The number of model layers.", units="nondim", &
135  static_value=nk_)
136  if (nk /= nk_) call mom_error(fatal, "verticalGridInit: " // &
137  "Mismatched number of layers NK_ between MOM_memory.h and param_file")
138 
139 #else
140  call get_param(param_file, mdl, "NK", nk, &
141  "The number of model layers.", units="nondim", fail_if_missing=.true.)
142 #endif
143  gv%ke = nk
144 
145  if (gv%Boussinesq) then
146  gv%H_to_kg_m2 = us%R_to_kg_m3*gv%Rho0 * gv%H_to_m
147  gv%kg_m2_to_H = 1.0 / gv%H_to_kg_m2
148  gv%m_to_H = 1.0 / gv%H_to_m
149  gv%Angstrom_H = gv%m_to_H * gv%Angstrom_m
150  gv%H_to_MKS = gv%H_to_m
151  else
152  gv%kg_m2_to_H = 1.0 / gv%H_to_kg_m2
153  gv%m_to_H = us%R_to_kg_m3*gv%Rho0 * gv%kg_m2_to_H
154  gv%H_to_m = gv%H_to_kg_m2 / (us%R_to_kg_m3*gv%Rho0)
155  gv%Angstrom_H = gv%Angstrom_m*1000.0*gv%kg_m2_to_H
156  gv%H_to_MKS = gv%H_to_kg_m2
157  endif
158  gv%H_subroundoff = 1e-20 * max(gv%Angstrom_H,gv%m_to_H*1e-17)
159  gv%H_to_Pa = gv%mks_g_Earth * gv%H_to_kg_m2
160 
161  gv%H_to_Z = gv%H_to_m * us%m_to_Z
162  gv%Z_to_H = us%Z_to_m * gv%m_to_H
163  gv%Angstrom_Z = us%m_to_Z * gv%Angstrom_m
164 
165  gv%H_to_RZ = gv%H_to_kg_m2 * us%kg_m3_to_R * us%m_to_Z
166  gv%RZ_to_H = gv%kg_m2_to_H * us%R_to_kg_m3 * us%Z_to_m
167 
168 ! Log derivative values.
169  call log_param(param_file, mdl, "M to THICKNESS", gv%m_to_H*h_rescale_factor)
170  call log_param(param_file, mdl, "M to THICKNESS rescaled by 2^-n", gv%m_to_H)
171  call log_param(param_file, mdl, "THICKNESS to M rescaled by 2^n", gv%H_to_m)
172 
173  allocate( gv%sInterface(nk+1) )
174  allocate( gv%sLayer(nk) )
175  allocate( gv%g_prime(nk+1) ) ; gv%g_prime(:) = 0.0
176  ! The extent of Rlay should be changed to nk?
177  allocate( gv%Rlay(nk+1) ) ; gv%Rlay(:) = 0.0
178 
179 end subroutine verticalgridinit
180 
181 !> Set the scaling factors for restart files to the scaling factors for this run.
182 subroutine fix_restart_scaling(GV)
183  type(verticalgrid_type), intent(inout) :: gv !< The ocean's vertical grid structure
184 
185  gv%m_to_H_restart = gv%m_to_H
186 end subroutine fix_restart_scaling
187 
188 !> Returns the model's thickness units, usually m or kg/m^2.
189 function get_thickness_units(GV)
190  character(len=48) :: get_thickness_units !< The vertical thickness units
191  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
192  ! This subroutine returns the appropriate units for thicknesses,
193  ! depending on whether the model is Boussinesq or not and the scaling for
194  ! the vertical thickness.
195 
196  if (gv%Boussinesq) then
197  get_thickness_units = "m"
198  else
199  get_thickness_units = "kg m-2"
200  endif
201 end function get_thickness_units
202 
203 !> Returns the model's thickness flux units, usually m^3/s or kg/s.
204 function get_flux_units(GV)
205  character(len=48) :: get_flux_units !< The thickness flux units
206  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
207  ! This subroutine returns the appropriate units for thickness fluxes,
208  ! depending on whether the model is Boussinesq or not and the scaling for
209  ! the vertical thickness.
210 
211  if (gv%Boussinesq) then
212  get_flux_units = "m3 s-1"
213  else
214  get_flux_units = "kg s-1"
215  endif
216 end function get_flux_units
217 
218 !> Returns the model's tracer flux units.
219 function get_tr_flux_units(GV, tr_units, tr_vol_conc_units, tr_mass_conc_units)
220  character(len=48) :: get_tr_flux_units !< The model's flux units
221  !! for a tracer.
222  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical
223  !! grid structure.
224  character(len=*), optional, intent(in) :: tr_units !< Units for a tracer, for example
225  !! Celsius or PSU.
226  character(len=*), optional, intent(in) :: tr_vol_conc_units !< The concentration units per unit
227  !! volume, for example if the units are
228  !! umol m-3, tr_vol_conc_units would
229  !! be umol.
230  character(len=*), optional, intent(in) :: tr_mass_conc_units !< The concentration units per unit
231  !! mass of sea water, for example if
232  !! the units are mol kg-1,
233  !! tr_vol_conc_units would be mol.
234 
235  ! This subroutine returns the appropriate units for thicknesses and fluxes,
236  ! depending on whether the model is Boussinesq or not and the scaling for
237  ! the vertical thickness.
238  integer :: cnt
239 
240  cnt = 0
241  if (present(tr_units)) cnt = cnt+1
242  if (present(tr_vol_conc_units)) cnt = cnt+1
243  if (present(tr_mass_conc_units)) cnt = cnt+1
244 
245  if (cnt == 0) call mom_error(fatal, "get_tr_flux_units: One of the three "//&
246  "arguments tr_units, tr_vol_conc_units, or tr_mass_conc_units "//&
247  "must be present.")
248  if (cnt > 1) call mom_error(fatal, "get_tr_flux_units: Only one of "//&
249  "tr_units, tr_vol_conc_units, and tr_mass_conc_units may be present.")
250  if (present(tr_units)) then
251  if (gv%Boussinesq) then
252  get_tr_flux_units = trim(tr_units)//" m3 s-1"
253  else
254  get_tr_flux_units = trim(tr_units)//" kg s-1"
255  endif
256  endif
257  if (present(tr_vol_conc_units)) then
258  if (gv%Boussinesq) then
259  get_tr_flux_units = trim(tr_vol_conc_units)//" s-1"
260  else
261  get_tr_flux_units = trim(tr_vol_conc_units)//" m-3 kg s-1"
262  endif
263  endif
264  if (present(tr_mass_conc_units)) then
265  if (gv%Boussinesq) then
266  get_tr_flux_units = trim(tr_mass_conc_units)//" kg-1 m3 s-1"
267  else
268  get_tr_flux_units = trim(tr_mass_conc_units)//" s-1"
269  endif
270  endif
271 
272 end function get_tr_flux_units
273 
274 !> This sets the coordinate data for the "layer mode" of the isopycnal model.
275 subroutine setverticalgridaxes( Rlay, GV, scale )
276  type(verticalgrid_type), intent(inout) :: gv !< The container for vertical grid data
277  real, dimension(GV%ke), intent(in) :: rlay !< The layer target density [R ~> kg m-3]
278  real, intent(in) :: scale !< A unit scaling factor for Rlay
279  ! Local variables
280  integer :: k, nk
281 
282  nk = gv%ke
283 
284  gv%zAxisLongName = 'Target Potential Density'
285  gv%zAxisUnits = 'kg m-3'
286  do k=1,nk ; gv%sLayer(k) = scale*rlay(k) ; enddo
287  if (nk > 1) then
288  gv%sInterface(1) = scale * (1.5*rlay(1) - 0.5*rlay(2))
289  do k=2,nk ; gv%sInterface(k) = scale * 0.5*( rlay(k-1) + rlay(k) ) ; enddo
290  gv%sInterface(nk+1) = scale * (1.5*rlay(nk) - 0.5*rlay(nk-1))
291  else
292  gv%sInterface(1) = 0.0 ; gv%sInterface(nk+1) = 2.0*scale*rlay(nk)
293  endif
294 
295 end subroutine setverticalgridaxes
296 
297 !> Deallocates the model's vertical grid structure.
298 subroutine verticalgridend( GV )
299  type(verticalgrid_type), pointer :: gv !< The ocean's vertical grid structure
300 
301  deallocate( gv%g_prime, gv%Rlay )
302  deallocate( gv%sInterface , gv%sLayer )
303  deallocate( gv )
304 
305 end subroutine verticalgridend
306 
307 end module mom_verticalgrid
mom_verticalgrid::setverticalgridaxes
subroutine, public setverticalgridaxes(Rlay, GV, scale)
This sets the coordinate data for the "layer mode" of the isopycnal model.
Definition: MOM_verticalGrid.F90:276
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_verticalgrid::get_tr_flux_units
character(len=48) function, public get_tr_flux_units(GV, tr_units, tr_vol_conc_units, tr_mass_conc_units)
Returns the model's tracer flux units.
Definition: MOM_verticalGrid.F90:220
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_verticalgrid::get_thickness_units
character(len=48) function, public get_thickness_units(GV)
Returns the model's thickness units, usually m or kg/m^2.
Definition: MOM_verticalGrid.F90:190
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_verticalgrid::fix_restart_scaling
subroutine, public fix_restart_scaling(GV)
Set the scaling factors for restart files to the scaling factors for this run.
Definition: MOM_verticalGrid.F90:183
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_verticalgrid::verticalgridinit
subroutine, public verticalgridinit(param_file, GV, US)
Allocates and initializes the ocean model vertical grid structure.
Definition: MOM_verticalGrid.F90:76
mom_verticalgrid::verticalgridend
subroutine, public verticalgridend(GV)
Deallocates the model's vertical grid structure.
Definition: MOM_verticalGrid.F90:299
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_verticalgrid::get_flux_units
character(len=48) function, public get_flux_units(GV)
Returns the model's thickness flux units, usually m^3/s or kg/s.
Definition: MOM_verticalGrid.F90:205
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2