10 implicit none ;
private
12 #include <MOM_memory.h>
35 character(len=40) :: zaxisunits
36 character(len=40) :: zaxislongname
38 real,
allocatable,
dimension(:) :: slayer
39 real,
allocatable,
dimension(:) :: sinterface
40 integer :: direction = 1
50 real,
allocatable,
dimension(:) :: &
51 g_prime, & !< The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2].
55 integer :: nk_rho_varies = 0
69 real :: m_to_h_restart = 0.0
83 integer :: nk, h_power
84 real :: h_rescale_factor
86 # include "version_variable.h"
87 character(len=16) :: mdl =
'MOM_verticalGrid'
89 if (
associated(gv))
call mom_error(fatal, &
90 'verticalGridInit: called with an associated GV pointer.')
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", &
123 gv%H_to_kg_m2 = gv%H_to_kg_m2 * h_rescale_factor
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
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_
133 call get_param(param_file, mdl,
"NK", nk, &
134 "The number of model layers.", units=
"nondim", &
136 if (nk /= nk_)
call mom_error(fatal,
"verticalGridInit: " // &
137 "Mismatched number of layers NK_ between MOM_memory.h and param_file")
140 call get_param(param_file, mdl,
"NK", nk, &
141 "The number of model layers.", units=
"nondim", fail_if_missing=.true.)
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
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
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
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
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
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)
173 allocate( gv%sInterface(nk+1) )
174 allocate( gv%sLayer(nk) )
175 allocate( gv%g_prime(nk+1) ) ; gv%g_prime(:) = 0.0
177 allocate( gv%Rlay(nk+1) ) ; gv%Rlay(:) = 0.0
185 gv%m_to_H_restart = gv%m_to_H
196 if (gv%Boussinesq)
then
211 if (gv%Boussinesq)
then
224 character(len=*),
optional,
intent(in) :: tr_units
226 character(len=*),
optional,
intent(in) :: tr_vol_conc_units
230 character(len=*),
optional,
intent(in) :: tr_mass_conc_units
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
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 "//&
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
257 if (
present(tr_vol_conc_units))
then
258 if (gv%Boussinesq)
then
264 if (
present(tr_mass_conc_units))
then
265 if (gv%Boussinesq)
then
277 real,
dimension(GV%ke),
intent(in) :: rlay
278 real,
intent(in) :: scale
284 gv%zAxisLongName =
'Target Potential Density'
285 gv%zAxisUnits =
'kg m-3'
286 do k=1,nk ; gv%sLayer(k) = scale*rlay(k) ;
enddo
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))
292 gv%sInterface(1) = 0.0 ; gv%sInterface(nk+1) = 2.0*scale*rlay(nk)
301 deallocate( gv%g_prime, gv%Rlay )
302 deallocate( gv%sInterface , gv%sLayer )