10 implicit none ;
private
12 #include <MOM_memory.h>
14 public verticalgridinit, verticalgridend
15 public setverticalgridaxes, fix_restart_scaling
16 public get_flux_units, get_thickness_units, get_tr_flux_units
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
65 real :: m_to_h_restart = 0.0
71 subroutine verticalgridinit( param_file, GV, US )
79 integer :: nk, h_power
80 real :: h_rescale_factor
82 # include "version_variable.h"
83 character(len=16) :: mdl =
'MOM_verticalGrid'
85 if (
associated(gv))
call mom_error(fatal, &
86 'verticalGridInit: called with an associated GV pointer.')
91 "Parameters providing information about the vertical grid.")
92 call get_param(param_file, mdl,
"G_EARTH", gv%mks_g_Earth, &
93 "The gravitational acceleration of the Earth.", &
94 units=
"m s-2", default = 9.80)
95 call get_param(param_file, mdl,
"RHO_0", gv%Rho0, &
96 "The mean ocean density used with BOUSSINESQ true to "//&
97 "calculate accelerations and the mass for conservation "//&
98 "properties, or with BOUSSINSEQ false to convert some "//&
99 "parameters from vertical units of m to kg m-2.", &
100 units=
"kg m-3", default=1035.0)
101 call get_param(param_file, mdl,
"BOUSSINESQ", gv%Boussinesq, &
102 "If true, make the Boussinesq approximation.", default=.true.)
103 call get_param(param_file, mdl,
"ANGSTROM", gv%Angstrom_m, &
104 "The minimum layer thickness, usually one-Angstrom.", &
105 units=
"m", default=1.0e-10)
106 call get_param(param_file, mdl,
"H_RESCALE_POWER", h_power, &
107 "An integer power of 2 that is used to rescale the model's "//&
108 "intenal units of thickness. Valid values range from -300 to 300.", &
109 units=
"nondim", default=0, debuggingparam=.true.)
110 if (abs(h_power) > 300)
call mom_error(fatal,
"verticalGridInit: "//&
111 "H_RESCALE_POWER is outside of the valid range of -300 to 300.")
112 h_rescale_factor = 1.0
113 if (h_power /= 0) h_rescale_factor = 2.0**h_power
114 if (.not.gv%Boussinesq)
then
115 call get_param(param_file, mdl,
"H_TO_KG_M2", gv%H_to_kg_m2,&
116 "A constant that translates thicknesses from the model's "//&
117 "internal units of thickness to kg m-2.", units=
"kg m-2 H-1", &
119 gv%H_to_kg_m2 = gv%H_to_kg_m2 * h_rescale_factor
121 call get_param(param_file, mdl,
"H_TO_M", gv%H_to_m, &
122 "A constant that translates the model's internal "//&
123 "units of thickness into m.", units=
"m H-1", default=1.0)
124 gv%H_to_m = gv%H_to_m * h_rescale_factor
126 gv%g_Earth = us%m_to_L**2*us%Z_to_m*us%T_to_s**2 * gv%mks_g_Earth
127 #ifdef STATIC_MEMORY_
129 call get_param(param_file, mdl,
"NK", nk, &
130 "The number of model layers.", units=
"nondim", &
132 if (nk /= nk_)
call mom_error(fatal,
"verticalGridInit: " // &
133 "Mismatched number of layers NK_ between MOM_memory.h and param_file")
136 call get_param(param_file, mdl,
"NK", nk, &
137 "The number of model layers.", units=
"nondim", fail_if_missing=.true.)
141 if (gv%Boussinesq)
then
142 gv%H_to_kg_m2 = gv%Rho0 * gv%H_to_m
143 gv%kg_m2_to_H = 1.0 / gv%H_to_kg_m2
144 gv%m_to_H = 1.0 / gv%H_to_m
145 gv%Angstrom_H = gv%m_to_H * gv%Angstrom_m
147 gv%kg_m2_to_H = 1.0 / gv%H_to_kg_m2
148 gv%m_to_H = gv%Rho0 * gv%kg_m2_to_H
149 gv%H_to_m = gv%H_to_kg_m2 / gv%Rho0
150 gv%Angstrom_H = gv%Angstrom_m*1000.0*gv%kg_m2_to_H
152 gv%H_subroundoff = 1e-20 * max(gv%Angstrom_H,gv%m_to_H*1e-17)
153 gv%H_to_Pa = gv%mks_g_Earth * gv%H_to_kg_m2
155 gv%H_to_Z = gv%H_to_m * us%m_to_Z
156 gv%Z_to_H = us%Z_to_m * gv%m_to_H
157 gv%Angstrom_Z = us%m_to_Z * gv%Angstrom_m
160 call log_param(param_file, mdl,
"M to THICKNESS", gv%m_to_H*h_rescale_factor)
161 call log_param(param_file, mdl,
"M to THICKNESS rescaled by 2^-n", gv%m_to_H)
162 call log_param(param_file, mdl,
"THICKNESS to M rescaled by 2^n", gv%H_to_m)
164 allocate( gv%sInterface(nk+1) )
165 allocate( gv%sLayer(nk) )
166 allocate( gv%g_prime(nk+1) ) ; gv%g_prime(:) = 0.0
168 allocate( gv%Rlay(nk+1) ) ; gv%Rlay(:) = 0.0
170 end subroutine verticalgridinit
173 subroutine fix_restart_scaling(GV)
176 gv%m_to_H_restart = gv%m_to_H
177 end subroutine fix_restart_scaling
180 function get_thickness_units(GV)
181 character(len=48) :: get_thickness_units
187 if (gv%Boussinesq)
then
188 get_thickness_units =
"m"
190 get_thickness_units =
"kg m-2"
192 end function get_thickness_units
195 function get_flux_units(GV)
196 character(len=48) :: get_flux_units
202 if (gv%Boussinesq)
then
203 get_flux_units =
"m3 s-1"
205 get_flux_units =
"kg s-1"
207 end function get_flux_units
210 function get_tr_flux_units(GV, tr_units, tr_vol_conc_units, tr_mass_conc_units)
211 character(len=48) :: get_tr_flux_units
215 character(len=*),
optional,
intent(in) :: tr_units
217 character(len=*),
optional,
intent(in) :: tr_vol_conc_units
221 character(len=*),
optional,
intent(in) :: tr_mass_conc_units
232 if (
present(tr_units)) cnt = cnt+1
233 if (
present(tr_vol_conc_units)) cnt = cnt+1
234 if (
present(tr_mass_conc_units)) cnt = cnt+1
236 if (cnt == 0)
call mom_error(fatal,
"get_tr_flux_units: One of the three "//&
237 "arguments tr_units, tr_vol_conc_units, or tr_mass_conc_units "//&
239 if (cnt > 1)
call mom_error(fatal,
"get_tr_flux_units: Only one of "//&
240 "tr_units, tr_vol_conc_units, and tr_mass_conc_units may be present.")
241 if (
present(tr_units))
then
242 if (gv%Boussinesq)
then
243 get_tr_flux_units = trim(tr_units)//
" m3 s-1"
245 get_tr_flux_units = trim(tr_units)//
" kg s-1"
248 if (
present(tr_vol_conc_units))
then
249 if (gv%Boussinesq)
then
250 get_tr_flux_units = trim(tr_vol_conc_units)//
" s-1"
252 get_tr_flux_units = trim(tr_vol_conc_units)//
" m-3 kg s-1"
255 if (
present(tr_mass_conc_units))
then
256 if (gv%Boussinesq)
then
257 get_tr_flux_units = trim(tr_mass_conc_units)//
" kg-1 m3 s-1"
259 get_tr_flux_units = trim(tr_mass_conc_units)//
" s-1"
263 end function get_tr_flux_units
266 subroutine setverticalgridaxes( Rlay, GV )
268 real,
dimension(GV%ke),
intent(in) :: rlay
274 gv%zAxisLongName =
'Target Potential Density'
275 gv%zAxisUnits =
'kg m-3'
276 do k=1,nk ; gv%sLayer(k) = rlay(k) ;
enddo
278 gv%sInterface(1) = 1.5*rlay(1) - 0.5*rlay(2)
279 do k=2,nk ; gv%sInterface(k) = 0.5*( rlay(k-1) + rlay(k) ) ;
enddo
280 gv%sInterface(nk+1) = 1.5*rlay(nk) - 0.5*rlay(nk-1)
282 gv%sInterface(1) = 0.0 ; gv%sInterface(nk+1) = 2.0*rlay(nk)
285 end subroutine setverticalgridaxes
288 subroutine verticalgridend( GV )
291 deallocate( gv%g_prime, gv%Rlay )
292 deallocate( gv%sInterface , gv%sLayer )
295 end subroutine verticalgridend