MOM6
mom_geothermal Module Reference

Detailed Description

Implemented geothermal heating at the ocean bottom.

Geothermal heating can be added either in a layered isopycnal mode, in which the heating raises the density of the layer to the target density of the layer above, and then moves the water into that layer, or in a simple Eulerian mode, in which the bottommost GEOTHERMAL_THICKNESS are heated. Geothermal heating will also provide a buoyant source of bottom TKE that can be used to further mix the near-bottom water. In cold fresh water lakes where heating increases density, water should be moved into deeper layers, but this is not implemented yet.

Data Types

type  geothermal_cs
 Control structure for geothermal heating. More...
 

Functions/Subroutines

subroutine, public geothermal (h, tv, dt, ea, eb, G, GV, US, CS, halo)
 Applies geothermal heating, including the movement of water between isopycnal layers to match the target densities. The heating is applied to the bottommost layers that occur within ### of the bottom. If the partial derivative of the coordinate density with temperature is positive or very small, the layers are simply heated in place. Any heat that can not be applied to the ocean is returned (WHERE)? More...
 
subroutine, public geothermal_init (Time, G, GV, US, param_file, diag, CS)
 Initialize parameters and allocate memory associated with the geothermal heating module. More...
 
subroutine, public geothermal_end (CS)
 Clean up and deallocate memory associated with the geothermal heating module. More...
 

Function/Subroutine Documentation

◆ geothermal()

subroutine, public mom_geothermal::geothermal ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, intent(in)  dt,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  ea,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  eb,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(geothermal_cs), pointer  CS,
integer, intent(in), optional  halo 
)

Applies geothermal heating, including the movement of water between isopycnal layers to match the target densities. The heating is applied to the bottommost layers that occur within ### of the bottom. If the partial derivative of the coordinate density with temperature is positive or very small, the layers are simply heated in place. Any heat that can not be applied to the ocean is returned (WHERE)?

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thicknesses [H ~> m or kg m-2]
[in,out]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]dtTime increment [T ~> s].
[in,out]eaThe amount of fluid moved downward into a layer; this should be increased due to mixed layer detrainment [H ~> m or kg m-2]
[in,out]ebThe amount of fluid moved upward into a layer; this should be increased due to mixed layer entrainment [H ~> m or kg m-2].
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to geothermal_init.
[in]haloHalo width over which to work

Definition at line 54 of file MOM_geothermal.F90.

54  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
55  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
56  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
57  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers
58  !! to any available thermodynamic
59  !! fields. Absent fields have NULL
60  !! ptrs.
61  real, intent(in) :: dt !< Time increment [T ~> s].
62  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: ea !< The amount of fluid moved
63  !! downward into a layer; this
64  !! should be increased due to mixed
65  !! layer detrainment [H ~> m or kg m-2]
66  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: eb !< The amount of fluid moved upward
67  !! into a layer; this should be
68  !! increased due to mixed layer
69  !! entrainment [H ~> m or kg m-2].
70  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
71  type(geothermal_CS), pointer :: CS !< The control structure returned by
72  !! a previous call to
73  !! geothermal_init.
74  integer, optional, intent(in) :: halo !< Halo width over which to work
75  ! Local variables
76  real, dimension(SZI_(G)) :: &
77  heat_rem, & ! remaining heat [H degC ~> m degC or kg degC m-2]
78  h_geo_rem, & ! remaining thickness to apply geothermal heating [H ~> m or kg m-2]
79  Rcv_BL, & ! coordinate density in the deepest variable density layer [R ~> kg m-3]
80  p_ref ! coordiante densities reference pressure [Pa]
81 
82  real, dimension(2) :: &
83  T2, S2, & ! temp and saln in the present and target layers [degC] and [ppt]
84  dRcv_dT_, & ! partial derivative of coordinate density wrt temp [R degC-1 ~> kg m-3 degC-1]
85  dRcv_dS_ ! partial derivative of coordinate density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
86 
87  real :: Angstrom, H_neglect ! small thicknesses [H ~> m or kg m-2]
88  real :: Rcv ! coordinate density of present layer [R ~> kg m-3]
89  real :: Rcv_tgt ! coordinate density of target layer [R ~> kg m-3]
90  real :: dRcv ! difference between Rcv and Rcv_tgt [R ~> kg m-3]
91  real :: dRcv_dT ! partial derivative of coordinate density wrt temp
92  ! in the present layer [R degC-1 ~> kg m-3 degC-1]; usually negative
93  real :: h_heated ! thickness that is being heated [H ~> m or kg m-2]
94  real :: heat_avail ! heating available for the present layer [degC H ~> degC m or degC kg m-2]
95  real :: heat_in_place ! heating to warm present layer w/o movement between layers
96  ! [degC H ~> degC m or degC kg m-2]
97  real :: heat_trans ! heating available to move water from present layer to target
98  ! layer [degC H ~> degC m or degC kg m-2]
99  real :: heating ! heating used to move water from present layer to target layer
100  ! [degC H ~> degC m or degC kg m-2]
101  ! 0 <= heating <= heat_trans
102  real :: h_transfer ! thickness moved between layers [H ~> m or kg m-2]
103  real :: wt_in_place ! relative weighting that goes from 0 to 1 [nondim]
104  real :: I_h ! inverse thickness [H-1 ~> m-1 or m2 kg-1]
105  real :: dTemp ! temperature increase in a layer [degC]
106  real :: Irho_cp ! inverse of heat capacity per unit layer volume
107  ! [degC H m2 J-1 ~> degC m3 J-1 or degC kg J-1]
108 
109  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_old ! Temperature of each layer
110  ! before any heat is added,
111  ! for diagnostics [degC]
112  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_old ! Thickness of each layer
113  ! before any heat is added,
114  ! for diagnostics [m or kg m-2]
115  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d ! Scratch variable used to
116  ! calculate change in heat
117  ! due to geothermal
118  real :: Idt ! inverse of the timestep [s-1]
119 
120  logical :: do_i(SZI_(G))
121  logical :: compute_h_old, compute_T_old
122  integer :: i, j, k, is, ie, js, je, nz, k2, i2
123  integer :: isj, iej, num_start, num_left, nkmb, k_tgt
124 
125  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
126  if (present(halo)) then
127  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
128  endif
129 
130  if (.not. associated(cs)) call mom_error(fatal, "MOM_geothermal: "//&
131  "Module must be initialized before it is used.")
132  if (.not.cs%apply_geothermal) return
133 
134  nkmb = gv%nk_rho_varies
135  irho_cp = 1.0 / (gv%H_to_kg_m2 * tv%C_p)
136  angstrom = gv%Angstrom_H
137  h_neglect = gv%H_subroundoff
138  p_ref(:) = tv%P_Ref
139  idt = 1.0 / dt
140 
141  if (.not.associated(tv%T)) call mom_error(fatal, "MOM geothermal: "//&
142  "Geothermal heating can only be applied if T & S are state variables.")
143 
144 ! do i=is,ie ; do j=js,je
145 ! resid(i,j) = tv%internal_heat(i,j)
146 ! enddo ; enddo
147 
148  ! Conditionals for tracking diagnostic depdendencies
149  compute_h_old = cs%id_internal_heat_h_tendency > 0 &
150  .or. cs%id_internal_heat_heat_tendency > 0 &
151  .or. cs%id_internal_heat_temp_tendency > 0
152 
153  compute_t_old = cs%id_internal_heat_heat_tendency > 0 &
154  .or. cs%id_internal_heat_temp_tendency > 0
155 
156  if (cs%id_internal_heat_heat_tendency > 0) work_3d(:,:,:) = 0.0
157  if (compute_h_old) h_old(:,:,:) = 0.0
158  if (compute_t_old) t_old(:,:,:) = 0.0
159 
160 !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,US,CS,dt,Irho_cp,nkmb,tv, &
161 !$OMP p_Ref,h,Angstrom,nz,H_neglect,eb, &
162 !$OMP compute_h_old,compute_T_old,h_old,T_old, &
163 !$OMP work_3d,Idt) &
164 !$OMP private(num_start,heat_rem,do_i,h_geo_rem,num_left,&
165 !$OMP isj,iej,Rcv_BL,h_heated,heat_avail,k_tgt, &
166 !$OMP Rcv_tgt,Rcv,dRcv_dT,T2,S2,dRcv_dT_, &
167 !$OMP dRcv_dS_,heat_in_place,heat_trans, &
168 !$OMP wt_in_place,dTemp,dRcv,h_transfer,heating, &
169 !$OMP I_h)
170 
171  do j=js,je
172  ! 1. Only work on columns that are being heated.
173  ! 2. Find the deepest layer with any mass.
174  ! 3. Find the partial derivative of locally referenced potential density
175  ! and coordinate density with temperature, and the density of the layer
176  ! and the layer above.
177  ! 4. Heat a portion of the bottommost layer until it matches the target
178  ! density of the layer above, and move it.
179  ! 4a. In the case of variable density layers, heat but do not move.
180  ! 5. If there is still heat left over, repeat for the next layer up.
181  ! This subroutine updates thickness, T & S, and increments eb accordingly.
182 
183  ! 6. If there is not enough mass in the ocean, pass some of the heat up
184  ! from the ocean via the frazil field?
185 
186  num_start = 0
187  do i=is,ie
188  heat_rem(i) = g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp))
189  do_i(i) = .true. ; if (heat_rem(i) <= 0.0) do_i(i) = .false.
190  if (do_i(i)) num_start = num_start + 1
191  h_geo_rem(i) = cs%Geothermal_thick * gv%m_to_H
192  enddo
193  if (num_start == 0) cycle
194  num_left = num_start
195 
196  ! Find the first and last columns that need to be worked on.
197  isj = ie+1 ; do i=is,ie ; if (do_i(i)) then ; isj = i ; exit ; endif ; enddo
198  iej = is-1 ; do i=ie,is,-1 ; if (do_i(i)) then ; iej = i ; exit ; endif ; enddo
199 
200  if (nkmb > 0) then
201  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref(:), &
202  rcv_bl(:), isj, iej-isj+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
203  else
204  rcv_bl(:) = -1.0
205  endif
206 
207  do k=nz,1,-1
208  do i=isj,iej ; if (do_i(i)) then
209 
210  ! Save temperature and thickness before any changes are made (for diagnostic)
211  if (compute_h_old) then
212  h_old(i,j,k) = h(i,j,k)
213  endif
214  if (compute_t_old) then
215  t_old(i,j,k) = tv%T(i,j,k)
216  endif
217 
218  if (h(i,j,k) > angstrom) then
219  if ((h(i,j,k)-angstrom) >= h_geo_rem(i)) then
220  h_heated = h_geo_rem(i)
221  heat_avail = heat_rem(i)
222  h_geo_rem(i) = 0.0
223  else
224  h_heated = (h(i,j,k)-angstrom)
225  heat_avail = heat_rem(i) * (h_heated / &
226  (h_geo_rem(i) + h_neglect))
227  h_geo_rem(i) = h_geo_rem(i) - h_heated
228  endif
229 
230  if (k<=nkmb .or. nkmb<=0) then
231  ! Simply heat the layer; convective adjustment occurs later
232  ! if necessary.
233  k_tgt = k
234  elseif ((k==nkmb+1) .or. (gv%Rlay(k-1) < rcv_bl(i))) then
235  ! Add enough heat to match the lowest buffer layer density.
236  k_tgt = nkmb
237  rcv_tgt = rcv_bl(i)
238  else
239  ! Add enough heat to match the target density of layer k-1.
240  k_tgt = k-1
241  rcv_tgt = gv%Rlay(k-1)
242  endif
243 
244  if (k<=nkmb .or. nkmb<=0) then
245  rcv = 0.0 ; drcv_dt = 0.0 ! Is this OK?
246  else
247  call calculate_density(tv%T(i,j,k), tv%S(i,j,k), tv%P_Ref, &
248  rcv, tv%eqn_of_state, scale=us%kg_m3_to_R)
249  t2(1) = tv%T(i,j,k) ; s2(1) = tv%S(i,j,k)
250  t2(2) = tv%T(i,j,k_tgt) ; s2(2) = tv%S(i,j,k_tgt)
251  call calculate_density_derivs(t2(:), s2(:), p_ref(:), &
252  drcv_dt_, drcv_ds_, 1, 2, tv%eqn_of_state, scale=us%kg_m3_to_R)
253  drcv_dt = 0.5*(drcv_dt_(1) + drcv_dt_(2))
254  endif
255 
256  if ((drcv_dt >= 0.0) .or. (k<=nkmb .or. nkmb<=0)) then
257  ! This applies to variable density layers.
258  heat_in_place = heat_avail
259  heat_trans = 0.0
260  elseif (drcv_dt <= cs%dRcv_dT_inplace) then
261  ! This is the option that usually applies in isopycnal coordinates.
262  heat_in_place = min(heat_avail, max(0.0, h(i,j,k) * &
263  ((gv%Rlay(k)-rcv) / drcv_dt)))
264  heat_trans = heat_avail - heat_in_place
265  else
266  ! wt_in_place should go from 0 to 1.
267  wt_in_place = (cs%dRcv_dT_inplace - drcv_dt) / cs%dRcv_dT_inplace
268  heat_in_place = max(wt_in_place*heat_avail, &
269  h(i,j,k) * ((gv%Rlay(k)-rcv) / drcv_dt) )
270  heat_trans = heat_avail - heat_in_place
271  endif
272 
273  if (heat_in_place > 0.0) then
274  ! This applies to variable density layers. In isopycnal coordinates
275  ! this only arises for relatively fresh water near the freezing
276  ! point, in which case heating in place will eventually cause things
277  ! to sort themselves out, if only because the water will warm to
278  ! the temperature of maximum density.
279  dtemp = heat_in_place / (h(i,j,k) + h_neglect)
280  tv%T(i,j,k) = tv%T(i,j,k) + dtemp
281  heat_rem(i) = heat_rem(i) - heat_in_place
282  rcv = rcv + drcv_dt * dtemp
283  endif
284 
285  if (heat_trans > 0.0) then
286  ! The second expression might never be used, but will avoid
287  ! division by 0.
288  drcv = max(rcv - rcv_tgt, 0.0)
289 
290  ! dTemp = -dRcv / dRcv_dT
291  ! h_transfer = min(heat_rem(i) / dTemp, h(i,j,k)-Angstrom)
292  if ((-drcv_dt * heat_trans) >= drcv * (h(i,j,k)-angstrom)) then
293  h_transfer = h(i,j,k) - angstrom
294  heating = (h_transfer * drcv) / (-drcv_dt)
295  ! Since not all the heat has been applied, return the fraction
296  ! of the layer thickness that has not yet been fully heated to
297  ! the h_geo_rem.
298  h_geo_rem(i) = h_geo_rem(i) + h_heated * &
299  ((heat_avail - (heating + heat_in_place)) / heat_avail)
300  else
301  h_transfer = (-drcv_dt * heat_trans) / drcv
302  heating = heat_trans
303  endif
304  heat_rem(i) = heat_rem(i) - heating
305 
306  i_h = 1.0 / (h(i,j,k_tgt) + h_transfer + h_neglect)
307  tv%T(i,j,k_tgt) = (h(i,j,k_tgt) * tv%T(i,j,k_tgt) + &
308  (h_transfer * tv%T(i,j,k) + heating)) * i_h
309  tv%S(i,j,k_tgt) = (h(i,j,k_tgt) * tv%S(i,j,k_tgt) + &
310  h_transfer * tv%S(i,j,k)) * i_h
311 
312  h(i,j,k) = h(i,j,k) - h_transfer
313  h(i,j,k_tgt) = h(i,j,k_tgt) + h_transfer
314  eb(i,j,k_tgt) = eb(i,j,k_tgt) + h_transfer
315  if (k_tgt < k-1) then
316  do k2 = k_tgt+1,k-1
317  eb(i,j,k2) = eb(i,j,k2) + h_transfer
318  enddo
319  endif
320  endif
321 
322  if (heat_rem(i) <= 0.0) then
323  do_i(i) = .false. ; num_left = num_left-1
324  ! For efficiency, uncomment these?
325  ! if ((i==isj) .and. (num_left > 0)) then
326  ! do i2=isj+1,iej ; if (do_i(i2)) then
327  ! isj = i2 ; exit ! Set the new starting value.
328  ! endif ; enddo
329  ! endif
330  ! if ((i==iej) .and. (num_left > 0)) then
331  ! do i2=iej-1,isj,-1 ; if (do_i(i2)) then
332  ! iej = i2 ; exit ! Set the new ending value.
333  ! endif ; enddo
334  ! endif
335  endif
336  endif
337 
338  ! Calculate heat tendency due to addition and transfer of internal heat
339  if (cs%id_internal_heat_heat_tendency > 0) then
340  work_3d(i,j,k) = ((gv%H_to_kg_m2 * tv%C_p) * idt) * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * t_old(i,j,k))
341  endif
342 
343  endif ; enddo
344  if (num_left <= 0) exit
345  enddo ! k-loop
346 
347  if (associated(tv%internal_heat)) then ; do i=is,ie
348  tv%internal_heat(i,j) = tv%internal_heat(i,j) + gv%H_to_kg_m2 * &
349  (g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp)) - heat_rem(i))
350  enddo ; endif
351  enddo ! j-loop
352 
353  ! Post diagnostic of 3D tendencies (heat, temperature, and thickness) due to internal heat
354  if (cs%id_internal_heat_heat_tendency > 0) then
355  call post_data(cs%id_internal_heat_heat_tendency, work_3d, cs%diag, alt_h=h_old)
356  endif
357  if (cs%id_internal_heat_temp_tendency > 0) then
358  do j=js,je; do i=is,ie; do k=nz,1,-1
359  work_3d(i,j,k) = idt * (tv%T(i,j,k) - t_old(i,j,k))
360  enddo; enddo; enddo
361  call post_data(cs%id_internal_heat_temp_tendency, work_3d, cs%diag, alt_h=h_old)
362  endif
363  if (cs%id_internal_heat_h_tendency > 0) then
364  do j=js,je; do i=is,ie; do k=nz,1,-1
365  work_3d(i,j,k) = idt * (h(i,j,k) - h_old(i,j,k))
366  enddo; enddo; enddo
367  call post_data(cs%id_internal_heat_h_tendency, work_3d, cs%diag, alt_h=h_old)
368  endif
369 
370 ! do i=is,ie ; do j=js,je
371 ! resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_kg_m2 * &
372 ! (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
373 ! enddo ; enddo
374 

References mom_error_handler::mom_error().

Referenced by mom_diabatic_driver::diabatic_ale(), mom_diabatic_driver::diabatic_ale_legacy(), and mom_diabatic_driver::layered_diabatic().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ geothermal_end()

subroutine, public mom_geothermal::geothermal_end ( type(geothermal_cs), pointer  CS)

Clean up and deallocate memory associated with the geothermal heating module.

Parameters
csGeothermal heating control structure that will be deallocated in this subroutine.

Definition at line 482 of file MOM_geothermal.F90.

482  type(geothermal_CS), pointer :: CS !< Geothermal heating control structure that
483  !! will be deallocated in this subroutine.
484 
485  if (associated(cs%geo_heat)) deallocate(cs%geo_heat)
486  if (associated(cs)) deallocate(cs)

◆ geothermal_init()

subroutine, public mom_geothermal::geothermal_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(geothermal_cs), pointer  CS 
)

Initialize parameters and allocate memory associated with the geothermal heating module.

Parameters
[in]timeCurrent model time.
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagStructure used to regulate diagnostic output.
csPointer pointing to the module control structure.

Definition at line 379 of file MOM_geothermal.F90.

379  type(time_type), target, intent(in) :: Time !< Current model time.
380  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
381  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
382  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
383  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
384  !! parameters.
385  type(diag_ctrl), target, intent(inout) :: diag !< Structure used to regulate diagnostic output.
386  type(geothermal_CS), pointer :: CS !< Pointer pointing to the module control
387  !! structure.
388 ! This include declares and sets the variable "version".
389 #include "version_variable.h"
390  character(len=40) :: mdl = "MOM_geothermal" ! module name
391  character(len=48) :: thickness_units
392  ! Local variables
393  character(len=200) :: inputdir, geo_file, filename, geotherm_var
394  real :: scale ! A constant heat flux or dimensionally rescaled scaling factor
395  ! [J m-2 T-1 ~> W m-2] or [s T-1 ~> 1]
396  integer :: i, j, isd, ied, jsd, jed, id
397  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
398 
399  if (associated(cs)) then
400  call mom_error(warning, "geothermal_init called with an associated"// &
401  "associated control structure.")
402  return
403  else ; allocate(cs) ; endif
404 
405  cs%diag => diag
406  cs%Time => time
407 
408  ! write parameters to the model log.
409  call log_version(param_file, mdl, version, "")
410  call get_param(param_file, mdl, "GEOTHERMAL_SCALE", scale, &
411  "The constant geothermal heat flux, a rescaling "//&
412  "factor for the heat flux read from GEOTHERMAL_FILE, or "//&
413  "0 to disable the geothermal heating.", &
414  units="W m-2 or various", default=0.0, scale=us%T_to_s)
415  cs%apply_geothermal = .not.(scale == 0.0)
416  if (.not.cs%apply_geothermal) return
417 
418  call safe_alloc_ptr(cs%geo_heat, isd, ied, jsd, jed) ; cs%geo_heat(:,:) = 0.0
419 
420  call get_param(param_file, mdl, "GEOTHERMAL_FILE", geo_file, &
421  "The file from which the geothermal heating is to be "//&
422  "read, or blank to use a constant heating rate.", default=" ")
423  call get_param(param_file, mdl, "GEOTHERMAL_THICKNESS", cs%geothermal_thick, &
424  "The thickness over which to apply geothermal heating.", &
425  units="m", default=0.1)
426  call get_param(param_file, mdl, "GEOTHERMAL_DRHO_DT_INPLACE", cs%dRcv_dT_inplace, &
427  "The value of drho_dT above which geothermal heating "//&
428  "simply heats water in place instead of moving it between "//&
429  "isopycnal layers. This must be negative.", &
430  units="kg m-3 K-1", scale=us%kg_m3_to_R, default=-0.01)
431  if (cs%dRcv_dT_inplace >= 0.0) call mom_error(fatal, "geothermal_init: "//&
432  "GEOTHERMAL_DRHO_DT_INPLACE must be negative.")
433 
434  if (len_trim(geo_file) >= 1) then
435  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
436  inputdir = slasher(inputdir)
437  filename = trim(inputdir)//trim(geo_file)
438  call log_param(param_file, mdl, "INPUTDIR/GEOTHERMAL_FILE", filename)
439  call get_param(param_file, mdl, "GEOTHERMAL_VARNAME", geotherm_var, &
440  "The name of the geothermal heating variable in "//&
441  "GEOTHERMAL_FILE.", default="geo_heat")
442  call mom_read_data(filename, trim(geotherm_var), cs%geo_heat, g%Domain)
443  do j=jsd,jed ; do i=isd,ied
444  cs%geo_heat(i,j) = (g%mask2dT(i,j) * scale) * cs%geo_heat(i,j)
445  enddo ; enddo
446  else
447  do j=jsd,jed ; do i=isd,ied
448  cs%geo_heat(i,j) = g%mask2dT(i,j) * scale
449  enddo ; enddo
450  endif
451  call pass_var(cs%geo_heat, g%domain)
452 
453  thickness_units = get_thickness_units(gv)
454 
455  ! post the static geothermal heating field
456  id = register_static_field('ocean_model', 'geo_heat', diag%axesT1, &
457  'Geothermal heat flux into ocean', 'W m-2', conversion=us%s_to_T, &
458  cmor_field_name='hfgeou', cmor_units='W m-2', &
459  cmor_standard_name='upward_geothermal_heat_flux_at_sea_floor', &
460  cmor_long_name='Upward geothermal heat flux at sea floor', &
461  x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
462  if (id > 0) call post_data(id, cs%geo_heat, diag, .true.)
463 
464  ! Diagnostic for tendencies due to internal heat (in 3d)
465  cs%id_internal_heat_heat_tendency=register_diag_field('ocean_model', &
466  'internal_heat_heat_tendency', diag%axesTL, time, &
467  'Heat tendency (in 3D) due to internal (geothermal) sources', &
468  'W m-2', conversion=us%s_to_T, v_extensive=.true.)
469  cs%id_internal_heat_temp_tendency=register_diag_field('ocean_model', &
470  'internal_heat_temp_tendency', diag%axesTL, time, &
471  'Temperature tendency (in 3D) due to internal (geothermal) sources', &
472  'degC s-1', conversion=us%s_to_T, v_extensive=.true.)
473  cs%id_internal_heat_h_tendency=register_diag_field('ocean_model', &
474  'internal_heat_h_tendency', diag%axesTL, time, &
475  'Thickness tendency (in 3D) due to internal (geothermal) sources', &
476  trim(thickness_units), conversion=gv%H_to_MKS*us%s_to_T, v_extensive=.true.)
477 

References mom_verticalgrid::get_thickness_units(), mom_error_handler::mom_error(), and mom_diag_mediator::register_static_field().

Referenced by mom_diabatic_driver::diabatic_driver_init().

Here is the call graph for this function:
Here is the caller graph for this function: