18 implicit none ;
private
20 #include <MOM_memory.h>
26 real :: drcv_dt_inplace
29 real,
pointer :: geo_heat(:,:) => null()
30 real :: geothermal_thick
32 logical :: apply_geothermal
36 type(time_type),
pointer :: time => null()
39 integer :: id_internal_heat_heat_tendency = -1
40 integer :: id_internal_heat_temp_tendency = -1
41 integer :: id_internal_heat_h_tendency = -1
53 subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo)
56 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
61 real,
intent(in) :: dt
62 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: ea
66 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: eb
74 integer,
optional,
intent(in) :: halo
76 real,
dimension(SZI_(G)) :: &
82 real,
dimension(2) :: &
87 real :: angstrom, h_neglect
109 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_old
112 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_old
115 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
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
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
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
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
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.")
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
153 compute_t_old = cs%id_internal_heat_heat_tendency > 0 &
154 .or. cs%id_internal_heat_temp_tendency > 0
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
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
193 if (num_start == 0) cycle
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
202 rcv_bl(:), isj, iej-isj+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
208 do i=isj,iej ;
if (do_i(i))
then
211 if (compute_h_old)
then
212 h_old(i,j,k) = h(i,j,k)
214 if (compute_t_old)
then
215 t_old(i,j,k) = tv%T(i,j,k)
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)
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
230 if (k<=nkmb .or. nkmb<=0)
then
234 elseif ((k==nkmb+1) .or. (gv%Rlay(k-1) < rcv_bl(i)))
then
241 rcv_tgt = gv%Rlay(k-1)
244 if (k<=nkmb .or. nkmb<=0)
then
245 rcv = 0.0 ; drcv_dt = 0.0
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)
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))
256 if ((drcv_dt >= 0.0) .or. (k<=nkmb .or. nkmb<=0))
then
258 heat_in_place = heat_avail
260 elseif (drcv_dt <= cs%dRcv_dT_inplace)
then
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
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
273 if (heat_in_place > 0.0)
then
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
285 if (heat_trans > 0.0)
then
288 drcv = max(rcv - rcv_tgt, 0.0)
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)
298 h_geo_rem(i) = h_geo_rem(i) + h_heated * &
299 ((heat_avail - (heating + heat_in_place)) / heat_avail)
301 h_transfer = (-drcv_dt * heat_trans) / drcv
304 heat_rem(i) = heat_rem(i) - heating
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
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
317 eb(i,j,k2) = eb(i,j,k2) + h_transfer
322 if (heat_rem(i) <= 0.0)
then
323 do_i(i) = .false. ; num_left = num_left-1
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))
344 if (num_left <= 0)
exit
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))
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)
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))
361 call post_data(cs%id_internal_heat_temp_tendency, work_3d, cs%diag, alt_h=h_old)
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))
367 call post_data(cs%id_internal_heat_h_tendency, work_3d, cs%diag, alt_h=h_old)
379 type(time_type),
target,
intent(in) :: time
385 type(
diag_ctrl),
target,
intent(inout) :: diag
389 #include "version_variable.h"
390 character(len=40) :: mdl =
"MOM_geothermal"
391 character(len=48) :: thickness_units
393 character(len=200) :: inputdir, geo_file, filename, geotherm_var
396 integer :: i, j, isd, ied, jsd, jed, id
397 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
399 if (
associated(cs))
then
400 call mom_error(warning,
"geothermal_init called with an associated"// &
401 "associated control structure.")
403 else ;
allocate(cs) ;
endif
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
418 call safe_alloc_ptr(cs%geo_heat, isd, ied, jsd, jed) ; cs%geo_heat(:,:) = 0.0
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.")
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)
447 do j=jsd,jed ;
do i=isd,ied
448 cs%geo_heat(i,j) = g%mask2dT(i,j) * scale
451 call pass_var(cs%geo_heat, g%domain)
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.)
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.)
485 if (
associated(cs%geo_heat))
deallocate(cs%geo_heat)
486 if (
associated(cs))
deallocate(cs)