MOM_geothermal.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Implemented geothermal heating at the ocean bottom.
7
8use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc
10use mom_domains, only : pass_var
11use mom_error_handler, only : mom_error, fatal, warning
12use mom_file_parser, only : get_param, log_param, log_version, param_file_type
13use mom_io, only : mom_read_data, slasher
14use mom_grid, only : ocean_grid_type
18use mom_eos, only : calculate_density, calculate_density_derivs, eos_domain
19use mom_eos, only : calculate_specific_vol_derivs
20
21implicit none ; private
22
23#include <MOM_memory.h>
24
26
27!> Control structure for geothermal heating
28type, public :: geothermal_cs ; private
29 logical :: initialized = .false. !< True if this control structure has been initialized.
30 real :: drcv_dt_inplace !< The value of dRcv_dT above which (dRcv_dT is negative) the
31 !! water is heated in place instead of moving upward between
32 !! layers in non-ALE layered mode [R C-1 ~> kg m-3 degC-1]
33 real, allocatable, dimension(:,:) :: geo_heat !< The geothermal heat flux [Q R Z T-1 ~> W m-2]
34 real :: geothermal_thick !< The thickness over which geothermal heating is
35 !! applied [H ~> m or kg m-2]
36 logical :: apply_geothermal !< If true, geothermal heating will be applied. This is false if
37 !! GEOTHERMAL_SCALE is 0 and there is no heat to apply.
38
39 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock
40 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the timing
41 !! timing of diagnostic output
42 integer :: id_internal_heat_heat_tendency = -1 !< ID for diagnostic of heat tendency
43 integer :: id_internal_heat_temp_tendency = -1 !< ID for diagnostic of temperature tendency
44 integer :: id_internal_heat_h_tendency = -1 !< ID for diagnostic of thickness tendency
45 integer :: id_geothermal_buoyancy_flux = -1 !< ID for diagnostic of bottom buoyancy flux
46end type geothermal_cs
47
48contains
49
50!> Applies geothermal heating, including the movement of water
51!! between isopycnal layers to match the target densities. The heating is
52!! applied to the bottommost layers that occur within GEOTHERMAL_THICKNESS of the bottom. If
53!! the partial derivative of the coordinate density with temperature is positive
54!! or very small, the layers are simply heated in place. Any heat that can not
55!! be applied to the ocean is returned (WHERE)?
56subroutine geothermal_entraining(h, tv, dt, ea, eb, G, GV, US, CS, halo)
57 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
58 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
59 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
60 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers
61 !! to any available thermodynamic fields.
62 real, intent(in) :: dt !< Time increment [T ~> s].
63 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: ea !< The amount of fluid moved
64 !! downward into a layer; this
65 !! should be increased due to mixed
66 !! layer detrainment [H ~> m or kg m-2]
67 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: eb !< The amount of fluid moved upward
68 !! into a layer; this should be
69 !! increased due to mixed layer
70 !! entrainment [H ~> m or kg m-2].
71 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
72 type(geothermal_cs), intent(in) :: cs !< The control structure returned by
73 !! a previous call to
74 !! geothermal_init.
75 integer, optional, intent(in) :: halo !< Halo width over which to work
76 ! Local variables
77 real, dimension(SZI_(G)) :: &
78 heat_rem, & ! remaining heat [H C ~> m degC or kg degC m-2]
79 h_geo_rem, & ! remaining thickness to apply geothermal heating [H ~> m or kg m-2]
80 rcv_bl, & ! coordinate density in the deepest variable density layer [R ~> kg m-3]
81 p_ref ! coordinate densities reference pressure [R L2 T-2 ~> Pa]
82
83 real, dimension(2) :: &
84 t2, s2, & ! temp and saln in the present and target layers [C ~> degC] and [S ~> ppt]
85 drcv_dt_, & ! partial derivative of coordinate density wrt temp [R C-1 ~> kg m-3 degC-1]
86 drcv_ds_ ! partial derivative of coordinate density wrt saln [R S-1 ~> kg m-3 ppt-1]
87
88 real :: angstrom, h_neglect ! small thicknesses [H ~> m or kg m-2]
89 real :: rcv ! coordinate density of present layer [R ~> kg m-3]
90 real :: rcv_tgt ! coordinate density of target layer [R ~> kg m-3]
91 real :: drcv ! difference between Rcv and Rcv_tgt [R ~> kg m-3]
92 real :: drcv_dt ! partial derivative of coordinate density wrt temp
93 ! in the present layer [R C-1 ~> kg m-3 degC-1]; usually negative
94 real :: h_heated ! thickness that is being heated [H ~> m or kg m-2]
95 real :: heat_avail ! heating available for the present layer [C H ~> degC m or degC kg m-2]
96 real :: heat_in_place ! heating to warm present layer w/o movement between layers
97 ! [C H ~> degC m or degC kg m-2]
98 real :: heat_trans ! heating available to move water from present layer to target
99 ! layer [C H ~> degC m or degC kg m-2]
100 real :: heating ! heating used to move water from present layer to target layer
101 ! [C H ~> degC m or degC kg m-2]
102 ! 0 <= heating <= heat_trans
103 real :: h_transfer ! thickness moved between layers [H ~> m or kg m-2]
104 real :: wt_in_place ! relative weighting that goes from 0 to 1 [nondim]
105 real :: i_h ! inverse thickness [H-1 ~> m-1 or m2 kg-1]
106 real :: dtemp ! temperature increase in a layer [C ~> degC]
107 real :: irho_cp ! inverse of heat capacity per unit layer volume
108 ! [C H Q-1 R-1 Z-1 ~> degC m3 J-1 or degC kg J-1]
109
110 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
111 t_old, & ! Temperature of each layer before any heat is added, for diagnostics [C ~> degC]
112 h_old, & ! Thickness of each layer before any heat is added, for diagnostics [H ~> m or kg m-2]
113 work_3d ! Scratch variable used to calculate changes due to geothermal [various]
114 real :: idt ! inverse of the timestep [T-1 ~> s-1]
115
116 logical :: do_i(szi_(g))
117 logical :: compute_h_old, compute_t_old
118 integer :: i, j, k, is, ie, js, je, nz, k2
119 integer :: isj, iej, num_left, nkmb, k_tgt
120
121 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
122 if (present(halo)) then
123 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
124 endif
125
126 if (.not. cs%initialized) call mom_error(fatal, "MOM_geothermal: "//&
127 "Module must be initialized before it is used.")
128 if (.not.cs%apply_geothermal) return
129
130 nkmb = gv%nk_rho_varies
131 irho_cp = 1.0 / (gv%H_to_RZ * tv%C_p)
132 angstrom = gv%Angstrom_H
133 h_neglect = gv%H_subroundoff
134 p_ref(:) = tv%P_Ref
135 idt = 1.0 / dt
136
137 if (.not.associated(tv%T)) call mom_error(fatal, "MOM geothermal_entraining: "//&
138 "Geothermal heating can only be applied if T & S are state variables.")
139
140! do j=js,je ; do i=is,ie
141! resid(i,j) = tv%internal_heat(i,j)
142! enddo ; enddo
143
144 ! Conditionals for tracking diagnostic depdendencies
145 compute_h_old = cs%id_internal_heat_h_tendency > 0 &
146 .or. cs%id_internal_heat_heat_tendency > 0 &
147 .or. cs%id_internal_heat_temp_tendency > 0
148
149 compute_t_old = cs%id_internal_heat_heat_tendency > 0 &
150 .or. cs%id_internal_heat_temp_tendency > 0
151
152 if (cs%id_internal_heat_heat_tendency > 0) work_3d(:,:,:) = 0.0
153
154 if (compute_h_old .or. compute_t_old) then ; do k=1,nz ; do j=js,je ; do i=is,ie
155 ! Save temperature and thickness before any changes are made (for diagnostics)
156 h_old(i,j,k) = h(i,j,k)
157 t_old(i,j,k) = tv%T(i,j,k)
158 enddo ; enddo ; enddo ; endif
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 h_old,T_old,work_3d,Idt) &
163!$OMP private(heat_rem,do_i,h_geo_rem,num_left, &
164!$OMP isj,iej,Rcv_BL,h_heated,heat_avail,k_tgt, &
165!$OMP Rcv_tgt,Rcv,dRcv_dT,T2,S2,dRcv_dT_, &
166!$OMP dRcv_dS_,heat_in_place,heat_trans, &
167!$OMP wt_in_place,dTemp,dRcv,h_transfer,heating, &
168!$OMP I_h)
169
170 do j=js,je
171 ! 1. Only work on columns that are being heated.
172 ! 2. Find the deepest layer with any mass.
173 ! 3. Find the partial derivative of locally referenced potential density
174 ! and coordinate density with temperature, and the density of the layer
175 ! and the layer above.
176 ! 4. Heat a portion of the bottommost layer until it matches the target
177 ! density of the layer above, and move it.
178 ! 4a. In the case of variable density layers, heat but do not move.
179 ! 5. If there is still heat left over, repeat for the next layer up.
180 ! This subroutine updates thickness, T & S, and increments eb accordingly.
181
182 ! 6. If there is not enough mass in the ocean, pass some of the heat up
183 ! from the ocean via the frazil field?
184
185 num_left = 0
186 do i=is,ie
187 heat_rem(i) = g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp))
188 do_i(i) = .true. ; if (heat_rem(i) <= 0.0) do_i(i) = .false.
189 if (do_i(i)) num_left = num_left + 1
190 h_geo_rem(i) = cs%Geothermal_thick
191 enddo
192 if (num_left == 0) cycle
193
194 ! Find the first and last columns that need to be worked on.
195 isj = ie+1 ; do i=is,ie ; if (do_i(i)) then ; isj = i ; exit ; endif ; enddo
196 iej = is-1 ; do i=ie,is,-1 ; if (do_i(i)) then ; iej = i ; exit ; endif ; enddo
197
198 if (nkmb > 0) then
199 call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref(:), rcv_bl(:), &
200 tv%eqn_of_state, (/isj-(g%isd-1),iej-(g%isd-1)/) )
201 else
202 rcv_bl(:) = -1.0
203 endif
204
205 do k=nz,1,-1
206 do i=isj,iej ; if (do_i(i)) then
207
208 if (h(i,j,k) > angstrom) then
209 if ((h(i,j,k)-angstrom) >= h_geo_rem(i)) then
210 h_heated = h_geo_rem(i)
211 heat_avail = heat_rem(i)
212 h_geo_rem(i) = 0.0
213 else
214 h_heated = (h(i,j,k)-angstrom)
215 heat_avail = heat_rem(i) * (h_heated / &
216 (h_geo_rem(i) + h_neglect))
217 h_geo_rem(i) = h_geo_rem(i) - h_heated
218 endif
219
220 if (k<=nkmb .or. nkmb<=0) then
221 ! Simply heat the layer; convective adjustment occurs later
222 ! if necessary.
223 k_tgt = k
224 elseif ((k==nkmb+1) .or. (gv%Rlay(k-1) < rcv_bl(i))) then
225 ! Add enough heat to match the lowest buffer layer density.
226 k_tgt = nkmb
227 rcv_tgt = rcv_bl(i)
228 else
229 ! Add enough heat to match the target density of layer k-1.
230 k_tgt = k-1
231 rcv_tgt = gv%Rlay(k-1)
232 endif
233
234 if (k<=nkmb .or. nkmb<=0) then
235 rcv = 0.0 ; drcv_dt = 0.0 ! Is this OK?
236 else
237 call calculate_density(tv%T(i,j,k), tv%S(i,j,k), tv%P_Ref, &
238 rcv, tv%eqn_of_state)
239 t2(1) = tv%T(i,j,k) ; s2(1) = tv%S(i,j,k)
240 t2(2) = tv%T(i,j,k_tgt) ; s2(2) = tv%S(i,j,k_tgt)
241 call calculate_density_derivs(t2(:), s2(:), p_ref(:), drcv_dt_, drcv_ds_, &
242 tv%eqn_of_state, (/1,2/) )
243 drcv_dt = 0.5*(drcv_dt_(1) + drcv_dt_(2))
244 endif
245
246 if ((drcv_dt >= 0.0) .or. (k<=nkmb .or. nkmb<=0)) then
247 ! This applies to variable density layers.
248 heat_in_place = heat_avail
249 heat_trans = 0.0
250 elseif (drcv_dt <= cs%dRcv_dT_inplace) then
251 ! This is the option that usually applies in isopycnal coordinates.
252 heat_in_place = min(heat_avail, max(0.0, h(i,j,k) * &
253 ((gv%Rlay(k)-rcv) / drcv_dt)))
254 heat_trans = heat_avail - heat_in_place
255 else
256 ! wt_in_place should go from 0 to 1.
257 wt_in_place = (cs%dRcv_dT_inplace - drcv_dt) / cs%dRcv_dT_inplace
258 heat_in_place = max(wt_in_place*heat_avail, &
259 h(i,j,k) * ((gv%Rlay(k)-rcv) / drcv_dt) )
260 heat_trans = heat_avail - heat_in_place
261 endif
262
263 if (heat_in_place > 0.0) then
264 ! This applies to variable density layers. In isopycnal coordinates
265 ! this only arises for relatively fresh water near the freezing
266 ! point, in which case heating in place will eventually cause things
267 ! to sort themselves out, if only because the water will warm to
268 ! the temperature of maximum density.
269 dtemp = heat_in_place / (h(i,j,k) + h_neglect)
270 tv%T(i,j,k) = tv%T(i,j,k) + dtemp
271 heat_rem(i) = heat_rem(i) - heat_in_place
272 rcv = rcv + drcv_dt * dtemp
273 endif
274
275 if (heat_trans > 0.0) then
276 ! The second expression might never be used, but will avoid
277 ! division by 0.
278 drcv = max(rcv - rcv_tgt, 0.0)
279
280 ! dTemp = -dRcv / dRcv_dT
281 ! h_transfer = min(heat_rem(i) / dTemp, h(i,j,k)-Angstrom)
282 if ((-drcv_dt * heat_trans) >= drcv * (h(i,j,k)-angstrom)) then
283 h_transfer = h(i,j,k) - angstrom
284 heating = (h_transfer * drcv) / (-drcv_dt)
285 ! Since not all the heat has been applied, return the fraction
286 ! of the layer thickness that has not yet been fully heated to
287 ! the h_geo_rem.
288 h_geo_rem(i) = h_geo_rem(i) + h_heated * &
289 ((heat_avail - (heating + heat_in_place)) / heat_avail)
290 else
291 h_transfer = (-drcv_dt * heat_trans) / drcv
292 heating = heat_trans
293 endif
294 heat_rem(i) = heat_rem(i) - heating
295
296 i_h = 1.0 / ((h(i,j,k_tgt) + h_neglect) + h_transfer)
297 tv%T(i,j,k_tgt) = ((h(i,j,k_tgt) + h_neglect) * tv%T(i,j,k_tgt) + &
298 (h_transfer * tv%T(i,j,k) + heating)) * i_h
299 tv%S(i,j,k_tgt) = ((h(i,j,k_tgt) + h_neglect) * tv%S(i,j,k_tgt) + &
300 h_transfer * tv%S(i,j,k)) * i_h
301
302 h(i,j,k) = h(i,j,k) - h_transfer
303 h(i,j,k_tgt) = h(i,j,k_tgt) + h_transfer
304 eb(i,j,k_tgt) = eb(i,j,k_tgt) + h_transfer
305 if (k_tgt < k-1) then
306 do k2 = k_tgt+1,k-1
307 eb(i,j,k2) = eb(i,j,k2) + h_transfer
308 enddo
309 endif
310 endif
311
312 if (heat_rem(i) <= 0.0) then
313 do_i(i) = .false. ; num_left = num_left-1
314 ! For efficiency, uncomment these?
315 ! if ((i==isj) .and. (num_left > 0)) then ; do i2=isj+1,iej ; if (do_i(i2)) then
316 ! isj = i2 ; exit ! Set the new starting value.
317 ! endif ; enddo ; endif
318 ! if ((i==iej) .and. (num_left > 0)) then ; do i2=iej-1,isj,-1 ; if (do_i(i2)) then
319 ! iej = i2 ; exit ! Set the new ending value.
320 ! endif ; enddo ; endif
321 endif
322 endif
323
324 ! Calculate heat tendency due to addition and transfer of internal heat
325 if (cs%id_internal_heat_heat_tendency > 0) then
326 work_3d(i,j,k) = ((gv%H_to_RZ*tv%C_p) * idt) * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * t_old(i,j,k))
327 endif
328
329 endif ; enddo
330 if (num_left <= 0) exit
331 enddo ! k-loop
332
333 if (associated(tv%internal_heat)) then ; do i=is,ie
334 tv%internal_heat(i,j) = tv%internal_heat(i,j) + gv%H_to_RZ * &
335 (g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp)) - heat_rem(i))
336 enddo ; endif
337 enddo ! j-loop
338
339 ! Post diagnostic of 3D tendencies (heat, temperature, and thickness) due to internal heat
340 if (cs%id_internal_heat_heat_tendency > 0) then
341 call post_data(cs%id_internal_heat_heat_tendency, work_3d, cs%diag, alt_h=h_old)
342 endif
343 if (cs%id_internal_heat_temp_tendency > 0) then
344 do k=1,nz ; do j=js,je ; do i=is,ie
345 work_3d(i,j,k) = idt * (tv%T(i,j,k) - t_old(i,j,k))
346 enddo ; enddo ; enddo
347 call post_data(cs%id_internal_heat_temp_tendency, work_3d, cs%diag, alt_h=h_old)
348 endif
349 if (cs%id_internal_heat_h_tendency > 0) then
350 do k=1,nz ; do j=js,je ; do i=is,ie
351 work_3d(i,j,k) = idt * (h(i,j,k) - h_old(i,j,k))
352 enddo ; enddo ; enddo
353 call post_data(cs%id_internal_heat_h_tendency, work_3d, cs%diag, alt_h=h_old)
354 endif
355
356! do j=js,je ; do i=is,ie
357! resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_RZ * &
358! (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
359! enddo ; enddo
360
361end subroutine geothermal_entraining
362
363!> Applies geothermal heating to the bottommost layers that occur within GEOTHERMAL_THICKNESS of
364!! the bottom, by simply heating the water in place. Any heat that can not be applied to the ocean
365!! is returned (WHERE)?
366subroutine geothermal_in_place(h, tv, dt, G, GV, US, CS, BFlx_geothermal, halo)
367 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
368 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
369 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
370 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers
371 !! to any available thermodynamic fields.
372 real, intent(in) :: dt !< Time increment [T ~> s].
373 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
374 type(geothermal_cs), intent(in) :: cs !< Geothermal heating control struct
375 real, dimension(SZI_(G), SZJ_(G)), intent(out) :: bflx_geothermal !< Geothermal buoyancy flux
376 !! in [Z2 T-3 ~> m2 s-3]
377 integer, optional, intent(in) :: halo !< Halo width over which to work
378
379
380 ! Local variables
381 real, dimension(SZI_(G)) :: &
382 heat_rem, & ! remaining heat [H C ~> m degC or kg degC m-2]
383 h_geo_rem, & ! remaining thickness to apply geothermal heating [H ~> m or kg m-2]
384 bottom_pressure, & ! Hydrostatic pressure in bottom layer [R L2 T-2 ~> Pa]
385 drhodt, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
386 drhods, & ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
387 dspvdt, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
388 dspvds ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
389
390 real :: angstrom, h_neglect ! small thicknesses [H ~> m or kg m-2]
391 real :: heat_here ! heating applied to the present layer [C H ~> degC m or degC kg m-2]
392 real :: dtemp ! temperature increase in a layer [C ~> degC]
393 real :: irho_cp ! inverse of heat capacity per unit layer volume
394 ! [C H Q-1 R-1 Z-1 ~> degC m3 J-1 or degC kg J-1]
395
396 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
397 dtdt_diag ! Diagnostic of temperature tendency [C T-1 ~> degC s-1] which might be
398 ! converted into a layer-integrated heat tendency [Q R Z T-1 ~> W m-2]
399 real :: idt ! inverse of the timestep [T-1 ~> s-1]
400 real :: h_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
401 real :: i_cp ! 1.0 / C_p [C Q-1 ~> kg degC J-1]
402 real :: i_rho0squared ! 1.0 / rho_0^2 (Boussinesq only) [R-2 ~> m6 kg-2]
403 logical :: do_any ! True if there is more to be done on the current j-row.
404 logical :: calc_diags ! True if diagnostic tendencies are needed.
405 logical :: nonbous ! If true, do not make the Boussinesq approximation.
406 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
407 integer :: i, j, k, is, ie, js, je, nz, isj, iej
408
409 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
410 if (present(halo)) then
411 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
412 endif
413
414 if (.not. cs%initialized) call mom_error(fatal, "MOM_geothermal: "//&
415 "Module must be initialized before it is used.")
416 if (.not.cs%apply_geothermal) return
417
418 nonbous = .not.(gv%Boussinesq .or. gv%semi_Boussinesq)
419 irho_cp = 1.0 / (gv%H_to_RZ * tv%C_p)
420 angstrom = gv%Angstrom_H
421 h_neglect = gv%H_subroundoff
422 idt = 1.0 / dt
423 h_to_pres = gv%H_to_RZ * gv%g_Earth
424 i_cp = 1. /tv%C_p
425 if (.not.nonbous) i_rho0squared = 1. / (gv%Rho0**2)
426 eosdom(:) = eos_domain(g%HI)
427
428 if (.not.associated(tv%T)) call mom_error(fatal, "MOM geothermal_in_place: "//&
429 "Geothermal heating can only be applied if T & S are state variables.")
430
431! do j=js,je ; do i=is,ie
432! resid(i,j) = tv%internal_heat(i,j)
433! enddo ; enddo
434
435 ! Conditionals for tracking diagnostic depdendencies
436 calc_diags = (cs%id_internal_heat_heat_tendency > 0) .or. (cs%id_internal_heat_temp_tendency > 0)
437 bflx_geothermal(:,:) = 0.0
438
439 if (calc_diags) dtdt_diag(:,:,:) = 0.0
440
441 !$OMP parallel do default(shared) private(heat_rem,do_any,h_geo_rem,isj,iej,heat_here,dTemp)
442 do j=js,je
443 bottom_pressure(:) = 0.0
444 do k=1,nz ; do i=is,ie
445 bottom_pressure(i) = bottom_pressure(i) + h_to_pres * h(i,j,k)
446 enddo ; enddo
447 if (nonbous) then
448 dspvdt(:) = 0.0
449 dspvds(:) = 0.0
450 call calculate_specific_vol_derivs(tv%T(:,j,nz), tv%S(:,j,nz), bottom_pressure, dspvdt, dspvds, &
451 tv%eqn_of_state, eosdom)
452 do i=is,ie
453 bflx_geothermal(i,j) = ( (gv%g_Earth_Z_T2 * dspvdt(i)) * (cs%geo_heat(i,j)*i_cp) ) * g%mask2dT(i,j)
454 enddo
455 else
456 drhodt(:) = 0.0
457 drhods(:) = 0.0
458 call calculate_density_derivs(tv%T(:,j,nz), tv%S(:,j,nz), bottom_pressure, drhodt, drhods, &
459 tv%eqn_of_state, eosdom)
460 do i=is,ie
461 bflx_geothermal(i,j) = - ( (gv%g_Earth_Z_T2*i_rho0squared) * ((i_cp*drhodt(i))*cs%geo_heat(i,j)) ) &
462 * g%mask2dT(i,j)
463 enddo
464 endif
465
466
467
468 ! Only work on columns that are being heated, and heat the near-bottom water.
469
470 ! If there is not enough mass in the ocean, pass some of the heat up
471 ! from the ocean via the frazil field?
472
473 do_any = .false.
474 do i=is,ie
475 heat_rem(i) = g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp))
476 if (heat_rem(i) > 0.0) do_any = .true.
477 h_geo_rem(i) = cs%Geothermal_thick
478 enddo
479 if (.not.do_any) cycle
480
481 ! Find the first and last columns that need to be worked on.
482 isj = ie+1 ; do i=is,ie ; if (heat_rem(i) > 0.0) then ; isj = i ; exit ; endif ; enddo
483 iej = is-1 ; do i=ie,is,-1 ; if (heat_rem(i) > 0.0) then ; iej = i ; exit ; endif ; enddo
484
485 do k=nz,1,-1
486 do_any = .false.
487 do i=isj,iej
488 if ((heat_rem(i) > 0.0) .and. (h(i,j,k) > angstrom)) then
489 ! Apply some or all of the remaining heat to this layer.
490 ! Convective adjustment occurs outside of this module if necessary.
491 if ((h(i,j,k)-angstrom) >= h_geo_rem(i)) then
492 heat_here = heat_rem(i)
493 h_geo_rem(i) = 0.0
494 heat_rem(i) = 0.0
495 else
496 heat_here = heat_rem(i) * ((h(i,j,k)-angstrom) / (h_geo_rem(i) + h_neglect))
497 h_geo_rem(i) = h_geo_rem(i) - (h(i,j,k)-angstrom)
498 heat_rem(i) = heat_rem(i) - heat_here
499 endif
500
501 dtemp = heat_here / (h(i,j,k) + h_neglect)
502 tv%T(i,j,k) = tv%T(i,j,k) + dtemp
503 if (calc_diags) dtdt_diag(i,j,k) = dtemp * idt
504 endif
505
506 if (heat_rem(i) > 0.0) do_any= .true.
507 enddo
508
509 if (.not.do_any) exit
510 enddo ! k-loop
511
512 if (associated(tv%internal_heat)) then ; do i=is,ie
513 tv%internal_heat(i,j) = tv%internal_heat(i,j) + gv%H_to_RZ * &
514 (g%mask2dT(i,j) * (cs%geo_heat(i,j) * (dt*irho_cp)) - heat_rem(i))
515 enddo ; endif
516 enddo ! j-loop
517
518 ! Post diagnostics of 3D tendencies of heat and temperature due to geothermal heat
519 if (cs%id_internal_heat_temp_tendency > 0) then
520 call post_data(cs%id_internal_heat_temp_tendency, dtdt_diag, cs%diag, alt_h=h)
521 endif
522 if (cs%id_internal_heat_heat_tendency > 0) then
523 do k=1,nz ; do j=js,je ; do i=is,ie
524 ! Dangerously reuse dTdt_diag for a related variable with different units, going from
525 ! units of [C T-1 ~> degC s-1] to units of [Q R Z T-1 ~> W m-2]
526 dtdt_diag(i,j,k) = (gv%H_to_RZ*tv%C_p) * (h(i,j,k) * dtdt_diag(i,j,k))
527 enddo ; enddo ; enddo
528 call post_data(cs%id_internal_heat_heat_tendency, dtdt_diag, cs%diag, alt_h=h)
529 endif
530 if (cs%id_geothermal_buoyancy_flux > 0) then
531 call post_data(cs%id_geothermal_buoyancy_flux, bflx_geothermal, cs%diag)
532 endif
533! do j=js,je ; do i=is,ie
534! resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_RZ * &
535! (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
536! enddo ; enddo
537
538end subroutine geothermal_in_place
539
540!> Initialize parameters and allocate memory associated with the geothermal heating module.
541subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm)
542 type(time_type), target, intent(in) :: time !< Current model time.
543 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
544 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
545 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
546 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
547 !! parameters.
548 type(diag_ctrl), target, intent(inout) :: diag !< Structure used to regulate diagnostic output.
549 type(geothermal_cs), intent(inout) :: cs !< Geothermal heating control struct
550 logical, optional, intent(in) :: usealealgorithm !< logical for whether to use ALE remapping
551
552! This include declares and sets the variable "version".
553#include "version_variable.h"
554 character(len=40) :: mdl = "MOM_geothermal" ! module name
555 character(len=48) :: thickness_units
556 ! Local variables
557 character(len=200) :: inputdir, geo_file, filename, geotherm_var
558 real :: geo_scale ! A constant heat flux or dimensionally rescaled geothermal flux scaling factor
559 ! [Q R Z T-1 ~> W m-2] or [Q R Z m2 s J-1 T-1 ~> nondim]
560 integer :: i, j, isd, ied, jsd, jed, id
561 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
562
563 cs%initialized = .true.
564 cs%diag => diag
565 cs%Time => time
566
567 ! write parameters to the model log.
568 call log_version(param_file, mdl, version, "")
569 call get_param(param_file, mdl, "GEOTHERMAL_SCALE", geo_scale, &
570 "The constant geothermal heat flux, a rescaling "//&
571 "factor for the heat flux read from GEOTHERMAL_FILE, or "//&
572 "0 to disable the geothermal heating.", &
573 units="W m-2 or various", default=0.0, scale=us%W_m2_to_QRZ_T)
574 cs%apply_geothermal = .not.(geo_scale == 0.0)
575 if (.not.cs%apply_geothermal) return
576
577 call safe_alloc_alloc(cs%geo_heat, isd, ied, jsd, jed) ; cs%geo_heat(:,:) = 0.0
578
579 call get_param(param_file, mdl, "GEOTHERMAL_FILE", geo_file, &
580 "The file from which the geothermal heating is to be "//&
581 "read, or blank to use a constant heating rate.", default=" ")
582 call get_param(param_file, mdl, "GEOTHERMAL_THICKNESS", cs%geothermal_thick, &
583 "The thickness over which to apply geothermal heating.", &
584 units="m", default=0.1, scale=gv%m_to_H)
585 call get_param(param_file, mdl, "GEOTHERMAL_DRHO_DT_INPLACE", cs%dRcv_dT_inplace, &
586 "The value of drho_dT above which geothermal heating "//&
587 "simply heats water in place instead of moving it between "//&
588 "isopycnal layers. This must be negative.", &
589 units="kg m-3 K-1", scale=us%kg_m3_to_R*us%C_to_degC, default=-0.01, &
590 do_not_log=((gv%nk_rho_varies<=0).or.(gv%nk_rho_varies>=gv%ke)) )
591 if (cs%dRcv_dT_inplace >= 0.0) call mom_error(fatal, "geothermal_init: "//&
592 "GEOTHERMAL_DRHO_DT_INPLACE must be negative.")
593
594 if (len_trim(geo_file) >= 1) then
595 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
596 inputdir = slasher(inputdir)
597 filename = trim(inputdir)//trim(geo_file)
598 call log_param(param_file, mdl, "INPUTDIR/GEOTHERMAL_FILE", filename)
599 call get_param(param_file, mdl, "GEOTHERMAL_VARNAME", geotherm_var, &
600 "The name of the geothermal heating variable in GEOTHERMAL_FILE.", &
601 default="geo_heat")
602 call mom_read_data(filename, trim(geotherm_var), cs%geo_heat, g%Domain)
603 do j=jsd,jed ; do i=isd,ied
604 cs%geo_heat(i,j) = (g%mask2dT(i,j) * geo_scale) * cs%geo_heat(i,j)
605 enddo ; enddo
606 else
607 do j=jsd,jed ; do i=isd,ied
608 cs%geo_heat(i,j) = g%mask2dT(i,j) * geo_scale
609 enddo ; enddo
610 endif
611 call pass_var(cs%geo_heat, g%domain)
612
613 thickness_units = get_thickness_units(gv)
614
615 ! post the static geothermal heating field
616 id = register_static_field('ocean_model', 'geo_heat', diag%axesT1, &
617 'Geothermal heat flux into ocean', 'W m-2', conversion=us%QRZ_T_to_W_m2, &
618 cmor_field_name='hfgeou', &
619 cmor_standard_name='upward_geothermal_heat_flux_at_sea_floor', &
620 cmor_long_name='Upward geothermal heat flux at sea floor', &
621 x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
622 if (id > 0) call post_data(id, cs%geo_heat, diag, .true.)
623
624 cs%id_geothermal_buoyancy_flux = register_diag_field('ocean_model', &
625 'geo_bflx', diag%axesT1, time, 'Geothermal buoyancy flux into ocean', &
626 'm2 s-3', conversion=us%Z_to_m**2*us%s_to_T**3)
627
628 ! Diagnostic for tendencies due to internal heat (in 3d)
629 cs%id_internal_heat_heat_tendency = register_diag_field('ocean_model', &
630 'internal_heat_heat_tendency', diag%axesTL, time, &
631 'Heat tendency (in 3D) due to internal (geothermal) sources', &
632 'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)
633 cs%id_internal_heat_temp_tendency = register_diag_field('ocean_model', &
634 'internal_heat_temp_tendency', diag%axesTL, time, &
635 'Temperature tendency (in 3D) due to internal (geothermal) sources', &
636 'degC s-1', conversion=us%C_to_degC*us%s_to_T, v_extensive=.true.)
637 if (.not.usealealgorithm) then
638 ! Do not offer this diagnostic if heating will be in place.
639 cs%id_internal_heat_h_tendency = register_diag_field('ocean_model', &
640 'internal_heat_h_tendency', diag%axesTL, time, &
641 'Thickness tendency (in 3D) due to internal (geothermal) sources', &
642 trim(thickness_units)//' s-1', conversion=gv%H_to_MKS*us%s_to_T, v_extensive=.true.)
643 endif
644
645end subroutine geothermal_init
646
647!> Clean up and deallocate memory associated with the geothermal heating module.
648subroutine geothermal_end(CS)
649 type(geothermal_cs), intent(inout) :: cs !< Geothermal heating control struct
650 if (allocated(cs%geo_heat)) deallocate(cs%geo_heat)
651end subroutine geothermal_end
652
653!> \namespace mom_geothermal
654!!
655!! Geothermal heating can be added either in a layered isopycnal mode, in which the heating raises the density
656!! of the layer to the target density of the layer above, and then moves the water into that layer, or in a
657!! simple Eulerian mode, in which the bottommost GEOTHERMAL_THICKNESS are heated. Geothermal heating will also
658!! provide a buoyant source of bottom TKE that can be used to further mix the near-bottom water. In cold fresh
659!! water lakes where heating increases density, water should be moved into deeper layers, but this is not
660!! implemented yet.
661
662end module mom_geothermal