MOM6
MOM_geothermal.F90
Go to the documentation of this file.
1 !> Implemented geothermal heating at the ocean bottom.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
8 use mom_domains, only : pass_var
9 use mom_error_handler, only : mom_error, fatal, warning
11 use mom_io, only : mom_read_data, slasher
12 use mom_grid, only : ocean_grid_type
17 
18 implicit none ; private
19 
20 #include <MOM_memory.h>
21 
23 
24 !> Control structure for geothermal heating
25 type, public :: geothermal_cs ; private
26  real :: drcv_dt_inplace !< The value of dRcv_dT above which (dRcv_dT is
27  !! negative) the water is heated in place instead
28  !! of moving upward between layers [R degC-1 ~> kg m-3 degC-1].
29  real, pointer :: geo_heat(:,:) => null() !< The geothermal heat flux [J m-2 T-1 ~> W m-2].
30  real :: geothermal_thick !< The thickness over which geothermal heating is
31  !! applied [m] (not [H]).
32  logical :: apply_geothermal !< If true, geothermal heating will be applied
33  !! otherwise GEOTHERMAL_SCALE has been set to 0 and
34  !! there is no heat to apply.
35 
36  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
37  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
38  !! regulate the timing of diagnostic output.
39  integer :: id_internal_heat_heat_tendency = -1 !< ID for diagnostic of heat tendency
40  integer :: id_internal_heat_temp_tendency = -1 !< ID for diagnostic of temperature tendency
41  integer :: id_internal_heat_h_tendency = -1 !< ID for diagnostic of thickness tendency
42 
43 end type geothermal_cs
44 
45 contains
46 
47 !> Applies geothermal heating, including the movement of water
48 !! between isopycnal layers to match the target densities. The heating is
49 !! applied to the bottommost layers that occur within ### of the bottom. If
50 !! the partial derivative of the coordinate density with temperature is positive
51 !! or very small, the layers are simply heated in place. Any heat that can not
52 !! be applied to the ocean is returned (WHERE)?
53 subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo)
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 
375 end subroutine geothermal
376 
377 !> Initialize parameters and allocate memory associated with the geothermal heating module.
378 subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS)
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 
478 end subroutine geothermal_init
479 
480 !> Clean up and deallocate memory associated with the geothermal heating module.
481 subroutine geothermal_end(CS)
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)
487 end subroutine geothermal_end
488 
489 !> \namespace mom_geothermal
490 !!
491 !! Geothermal heating can be added either in a layered isopycnal mode, in which the heating raises the density
492 !! of the layer to the target density of the layer above, and then moves the water into that layer, or in a
493 !! simple Eulerian mode, in which the bottommost GEOTHERMAL_THICKNESS are heated. Geothermal heating will also
494 !! provide a buoyant source of bottom TKE that can be used to further mix the near-bottom water. In cold fresh
495 !! water lakes where heating increases density, water should be moved into deeper layers, but this is not
496 !! implemented yet.
497 
498 end module mom_geothermal
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
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_geothermal::geothermal
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 tar...
Definition: MOM_geothermal.F90:54
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_geothermal::geothermal_init
subroutine, public geothermal_init(Time, G, GV, US, param_file, diag, CS)
Initialize parameters and allocate memory associated with the geothermal heating module.
Definition: MOM_geothermal.F90:379
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
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_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_geothermal::geothermal_end
subroutine, public geothermal_end(CS)
Clean up and deallocate memory associated with the geothermal heating module.
Definition: MOM_geothermal.F90:482
mom_geothermal
Implemented geothermal heating at the ocean bottom.
Definition: MOM_geothermal.F90:2
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_diag_mediator::register_static_field
integer function, public register_static_field(module_name, field_name, axes, long_name, units, missing_value, range, mask_variant, standard_name, do_not_log, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, area, x_cell_method, y_cell_method, area_cell_method, conversion)
Registers a static diagnostic, returning an integer handle.
Definition: MOM_diag_mediator.F90:2700
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_geothermal::geothermal_cs
Control structure for geothermal heating.
Definition: MOM_geothermal.F90:25
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.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
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60