portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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. | |
| 6 | module MOM_geothermal | |
| 7 | ||
| 8 | use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc | |
| 9 | use MOM_diag_mediator, only : register_static_field, time_type, diag_ctrl | |
| 10 | use MOM_domains, only : pass_var | |
| 11 | use MOM_error_handler, only : MOM_error, FATAL, WARNING | |
| 12 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 13 | use MOM_io, only : MOM_read_data, slasher | |
| 14 | use MOM_grid, only : ocean_grid_type | |
| 15 | use MOM_unit_scaling, only : unit_scale_type | |
| 16 | use MOM_variables, only : thermo_var_ptrs | |
| 17 | use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units | |
| 18 | use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain | |
| 19 | use MOM_EOS, only : calculate_specific_vol_derivs | |
| 20 | ||
| 21 | implicit none ; private | |
| 22 | ||
| 23 | #include <MOM_memory.h> | |
| 24 | ||
| 25 | public geothermal_entraining, geothermal_in_place, geothermal_init, geothermal_end | |
| 26 | ||
| 27 | !> Control structure for geothermal heating | |
| 28 | type, 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 | |
| 46 | end type geothermal_CS | |
| 47 | ||
| 48 | contains | |
| 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)? | |
| 56 | 0 | subroutine 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 | 0 | heat_rem, & ! remaining heat [H C ~> m degC or kg degC m-2] |
| 79 | 0 | h_geo_rem, & ! remaining thickness to apply geothermal heating [H ~> m or kg m-2] |
| 80 | 0 | Rcv_BL, & ! coordinate density in the deepest variable density layer [R ~> kg m-3] |
| 81 | 0 | 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 | 0 | T_old, & ! Temperature of each layer before any heat is added, for diagnostics [C ~> degC] |
| 112 | 0 | h_old, & ! Thickness of each layer before any heat is added, for diagnostics [H ~> m or kg m-2] |
| 113 | 0 | 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 122 | 0 | if (present(halo)) then |
| 123 | 0 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo |
| 124 | endif | |
| 125 | ||
| 126 | 0 | if (.not. CS%initialized) call MOM_error(FATAL, "MOM_geothermal: "//& |
| 127 | 0 | "Module must be initialized before it is used.") |
| 128 | 0 | if (.not.CS%apply_geothermal) return |
| 129 | ||
| 130 | 0 | nkmb = GV%nk_rho_varies |
| 131 | 0 | Irho_cp = 1.0 / (GV%H_to_RZ * tv%C_p) |
| 132 | 0 | Angstrom = GV%Angstrom_H |
| 133 | 0 | H_neglect = GV%H_subroundoff |
| 134 | 0 | p_ref(:) = tv%P_Ref |
| 135 | 0 | Idt = 1.0 / dt |
| 136 | ||
| 137 | 0 | if (.not.associated(tv%T)) call MOM_error(FATAL, "MOM geothermal_entraining: "//& |
| 138 | 0 | "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 | 0 | .or. CS%id_internal_heat_temp_tendency > 0 |
| 148 | ||
| 149 | compute_T_old = CS%id_internal_heat_heat_tendency > 0 & | |
| 150 | 0 | .or. CS%id_internal_heat_temp_tendency > 0 |
| 151 | ||
| 152 | 0 | if (CS%id_internal_heat_heat_tendency > 0) work_3d(:,:,:) = 0.0 |
| 153 | ||
| 154 | 0 | 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 | 0 | h_old(i,j,k) = h(i,j,k) |
| 157 | 0 | 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 | 0 | 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 | 0 | num_left = 0 |
| 186 | 0 | do i=is,ie |
| 187 | 0 | heat_rem(i) = G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)) |
| 188 | 0 | do_i(i) = .true. ; if (heat_rem(i) <= 0.0) do_i(i) = .false. |
| 189 | 0 | if (do_i(i)) num_left = num_left + 1 |
| 190 | 0 | h_geo_rem(i) = CS%Geothermal_thick |
| 191 | enddo | |
| 192 | 0 | if (num_left == 0) cycle |
| 193 | ||
| 194 | ! Find the first and last columns that need to be worked on. | |
| 195 | 0 | isj = ie+1 ; do i=is,ie ; if (do_i(i)) then ; isj = i ; exit ; endif ; enddo |
| 196 | 0 | iej = is-1 ; do i=ie,is,-1 ; if (do_i(i)) then ; iej = i ; exit ; endif ; enddo |
| 197 | ||
| 198 | 0 | if (nkmb > 0) then |
| 199 | call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_Ref(:), Rcv_BL(:), & | |
| 200 | 0 | tv%eqn_of_state, (/isj-(G%isd-1),iej-(G%isd-1)/) ) |
| 201 | else | |
| 202 | 0 | Rcv_BL(:) = -1.0 |
| 203 | endif | |
| 204 | ||
| 205 | 0 | do k=nz,1,-1 |
| 206 | 0 | do i=isj,iej ; if (do_i(i)) then |
| 207 | ||
| 208 | 0 | if (h(i,j,k) > Angstrom) then |
| 209 | 0 | if ((h(i,j,k)-Angstrom) >= h_geo_rem(i)) then |
| 210 | 0 | h_heated = h_geo_rem(i) |
| 211 | 0 | heat_avail = heat_rem(i) |
| 212 | 0 | h_geo_rem(i) = 0.0 |
| 213 | else | |
| 214 | 0 | h_heated = (h(i,j,k)-Angstrom) |
| 215 | heat_avail = heat_rem(i) * (h_heated / & | |
| 216 | 0 | (h_geo_rem(i) + H_neglect)) |
| 217 | 0 | h_geo_rem(i) = h_geo_rem(i) - h_heated |
| 218 | endif | |
| 219 | ||
| 220 | 0 | if (k<=nkmb .or. nkmb<=0) then |
| 221 | ! Simply heat the layer; convective adjustment occurs later | |
| 222 | ! if necessary. | |
| 223 | 0 | k_tgt = k |
| 224 | 0 | 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 | 0 | k_tgt = nkmb |
| 227 | 0 | Rcv_tgt = Rcv_BL(i) |
| 228 | else | |
| 229 | ! Add enough heat to match the target density of layer k-1. | |
| 230 | 0 | k_tgt = k-1 |
| 231 | 0 | Rcv_tgt = GV%Rlay(k-1) |
| 232 | endif | |
| 233 | ||
| 234 | 0 | if (k<=nkmb .or. nkmb<=0) then |
| 235 | 0 | 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 | 0 | Rcv, tv%eqn_of_state) |
| 239 | 0 | T2(1) = tv%T(i,j,k) ; S2(1) = tv%S(i,j,k) |
| 240 | 0 | 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 | 0 | tv%eqn_of_state, (/1,2/) ) |
| 243 | 0 | dRcv_dT = 0.5*(dRcv_dT_(1) + dRcv_dT_(2)) |
| 244 | endif | |
| 245 | ||
| 246 | 0 | if ((dRcv_dT >= 0.0) .or. (k<=nkmb .or. nkmb<=0)) then |
| 247 | ! This applies to variable density layers. | |
| 248 | 0 | heat_in_place = heat_avail |
| 249 | 0 | heat_trans = 0.0 |
| 250 | 0 | 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 | 0 | ((GV%Rlay(k)-Rcv) / dRcv_dT))) |
| 254 | 0 | heat_trans = heat_avail - heat_in_place |
| 255 | else | |
| 256 | ! wt_in_place should go from 0 to 1. | |
| 257 | 0 | wt_in_place = (CS%dRcv_dT_inplace - dRcv_dT) / CS%dRcv_dT_inplace |
| 258 | heat_in_place = max(wt_in_place*heat_avail, & | |
| 259 | 0 | h(i,j,k) * ((GV%Rlay(k)-Rcv) / dRcv_dT) ) |
| 260 | 0 | heat_trans = heat_avail - heat_in_place |
| 261 | endif | |
| 262 | ||
| 263 | 0 | 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 | 0 | dTemp = heat_in_place / (h(i,j,k) + H_neglect) |
| 270 | 0 | tv%T(i,j,k) = tv%T(i,j,k) + dTemp |
| 271 | 0 | heat_rem(i) = heat_rem(i) - heat_in_place |
| 272 | 0 | Rcv = Rcv + dRcv_dT * dTemp |
| 273 | endif | |
| 274 | ||
| 275 | 0 | if (heat_trans > 0.0) then |
| 276 | ! The second expression might never be used, but will avoid | |
| 277 | ! division by 0. | |
| 278 | 0 | 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 | 0 | if ((-dRcv_dT * heat_trans) >= dRcv * (h(i,j,k)-Angstrom)) then |
| 283 | 0 | h_transfer = h(i,j,k) - Angstrom |
| 284 | 0 | 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 | 0 | ((heat_avail - (heating + heat_in_place)) / heat_avail) |
| 290 | else | |
| 291 | 0 | h_transfer = (-dRcv_dT * heat_trans) / dRcv |
| 292 | 0 | heating = heat_trans |
| 293 | endif | |
| 294 | 0 | heat_rem(i) = heat_rem(i) - heating |
| 295 | ||
| 296 | 0 | 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 | 0 | (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 | 0 | h_transfer * tv%S(i,j,k)) * I_h |
| 301 | ||
| 302 | 0 | h(i,j,k) = h(i,j,k) - h_transfer |
| 303 | 0 | h(i,j,k_tgt) = h(i,j,k_tgt) + h_transfer |
| 304 | 0 | eb(i,j,k_tgt) = eb(i,j,k_tgt) + h_transfer |
| 305 | 0 | if (k_tgt < k-1) then |
| 306 | 0 | do k2 = k_tgt+1,k-1 |
| 307 | 0 | eb(i,j,k2) = eb(i,j,k2) + h_transfer |
| 308 | enddo | |
| 309 | endif | |
| 310 | endif | |
| 311 | ||
| 312 | 0 | if (heat_rem(i) <= 0.0) then |
| 313 | 0 | 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 | 0 | if (CS%id_internal_heat_heat_tendency > 0) then |
| 326 | 0 | 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 | 0 | if (num_left <= 0) exit |
| 331 | enddo ! k-loop | |
| 332 | ||
| 333 | 0 | 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 | 0 | (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 | 0 | if (CS%id_internal_heat_heat_tendency > 0) then |
| 341 | 0 | call post_data(CS%id_internal_heat_heat_tendency, work_3d, CS%diag, alt_h=h_old) |
| 342 | endif | |
| 343 | 0 | if (CS%id_internal_heat_temp_tendency > 0) then |
| 344 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 345 | 0 | work_3d(i,j,k) = Idt * (tv%T(i,j,k) - T_old(i,j,k)) |
| 346 | enddo ; enddo ; enddo | |
| 347 | 0 | call post_data(CS%id_internal_heat_temp_tendency, work_3d, CS%diag, alt_h=h_old) |
| 348 | endif | |
| 349 | 0 | if (CS%id_internal_heat_h_tendency > 0) then |
| 350 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 351 | 0 | work_3d(i,j,k) = Idt * (h(i,j,k) - h_old(i,j,k)) |
| 352 | enddo ; enddo ; enddo | |
| 353 | 0 | 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 | ||
| 361 | end 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)? | |
| 366 | 12 | subroutine 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 | 24 | heat_rem, & ! remaining heat [H C ~> m degC or kg degC m-2] |
| 383 | 24 | h_geo_rem, & ! remaining thickness to apply geothermal heating [H ~> m or kg m-2] |
| 384 | 12 | bottom_pressure, & ! Hydrostatic pressure in bottom layer [R L2 T-2 ~> Pa] |
| 385 | 24 | dRhodT, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] |
| 386 | 24 | dRhodS, & ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1] |
| 387 | 24 | dSpVdT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] |
| 388 | 24 | 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 | 12 | 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 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 410 | 12 | if (present(halo)) then |
| 411 | 12 | is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo |
| 412 | endif | |
| 413 | ||
| 414 | 12 | if (.not. CS%initialized) call MOM_error(FATAL, "MOM_geothermal: "//& |
| 415 | 0 | "Module must be initialized before it is used.") |
| 416 | 12 | if (.not.CS%apply_geothermal) return |
| 417 | ||
| 418 | 12 | nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq) |
| 419 | 12 | Irho_cp = 1.0 / (GV%H_to_RZ * tv%C_p) |
| 420 | 12 | Angstrom = GV%Angstrom_H |
| 421 | 12 | H_neglect = GV%H_subroundoff |
| 422 | 12 | Idt = 1.0 / dt |
| 423 | 12 | H_to_pres = GV%H_to_RZ * GV%g_Earth |
| 424 | 12 | I_Cp = 1. /tv%C_p |
| 425 | 12 | if (.not.nonBous) I_Rho0squared = 1. / (GV%Rho0**2) |
| 426 | 12 | EOSdom(:) = EOS_domain(G%HI) |
| 427 | ||
| 428 | 12 | if (.not.associated(tv%T)) call MOM_error(FATAL, "MOM geothermal_in_place: "//& |
| 429 | 0 | "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 | 12 | calc_diags = (CS%id_internal_heat_heat_tendency > 0) .or. (CS%id_internal_heat_temp_tendency > 0) |
| 437 | 105276 | BFlx_geothermal(:,:) = 0.0 |
| 438 | ||
| 439 | 12 | 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 | 756 | do j=js,je |
| 443 | 95976 | bottom_pressure(:) = 0.0 |
| 444 | 6864144 | do k=1,nz ; do i=is,ie |
| 445 | 6863400 | bottom_pressure(i) = bottom_pressure(i) + H_to_pres * h(i,j,k) |
| 446 | enddo ; enddo | |
| 447 | 744 | if (nonBous) then |
| 448 | 0 | dSpVdT(:) = 0.0 |
| 449 | 0 | dSpVdS(:) = 0.0 |
| 450 | call calculate_specific_vol_derivs(tv%T(:,j,nz), tv%S(:,j,nz), bottom_pressure, dSpVdT, dSpVdS, & | |
| 451 | 0 | tv%eqn_of_state, EOSdom) |
| 452 | 0 | do i=is,ie |
| 453 | 0 | 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 | 95976 | dRhodT(:) = 0.0 |
| 457 | 95976 | dRhodS(:) = 0.0 |
| 458 | call calculate_density_derivs(tv%T(:,j,nz), tv%S(:,j,nz), bottom_pressure, dRhodT, dRhodS, & | |
| 459 | 744 | tv%eqn_of_state, EOSdom) |
| 460 | 91512 | 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 | 91512 | * 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 | 744 | do_any = .false. |
| 474 | 91512 | do i=is,ie |
| 475 | 90768 | heat_rem(i) = G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)) |
| 476 | 90768 | if (heat_rem(i) > 0.0) do_any = .true. |
| 477 | 91512 | h_geo_rem(i) = CS%Geothermal_thick |
| 478 | enddo | |
| 479 | 744 | if (.not.do_any) cycle |
| 480 | ||
| 481 | ! Find the first and last columns that need to be worked on. | |
| 482 | 7476 | isj = ie+1 ; do i=is,ie ; if (heat_rem(i) > 0.0) then ; isj = i ; exit ; endif ; enddo |
| 483 | 7476 | iej = is-1 ; do i=ie,is,-1 ; if (heat_rem(i) > 0.0) then ; iej = i ; exit ; endif ; enddo |
| 484 | ||
| 485 | 43164 | do k=nz,1,-1 |
| 486 | 43164 | do_any = .false. |
| 487 | 4428588 | do i=isj,iej |
| 488 | 4385424 | 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 | 1476973 | if ((h(i,j,k)-Angstrom) >= h_geo_rem(i)) then |
| 492 | 60600 | heat_here = heat_rem(i) |
| 493 | 60600 | h_geo_rem(i) = 0.0 |
| 494 | 60600 | heat_rem(i) = 0.0 |
| 495 | else | |
| 496 | 1416373 | heat_here = heat_rem(i) * ((h(i,j,k)-Angstrom) / (h_geo_rem(i) + H_neglect)) |
| 497 | 1416373 | h_geo_rem(i) = h_geo_rem(i) - (h(i,j,k)-Angstrom) |
| 498 | 1416373 | heat_rem(i) = heat_rem(i) - heat_here |
| 499 | endif | |
| 500 | ||
| 501 | 1476973 | dTemp = heat_here / (h(i,j,k) + H_neglect) |
| 502 | 1476973 | tv%T(i,j,k) = tv%T(i,j,k) + dTemp |
| 503 | 1476973 | if (calc_diags) dTdt_diag(i,j,k) = dTemp * Idt |
| 504 | endif | |
| 505 | ||
| 506 | 4428588 | if (heat_rem(i) > 0.0) do_any= .true. |
| 507 | enddo | |
| 508 | ||
| 509 | 43164 | if (.not.do_any) exit |
| 510 | enddo ! k-loop | |
| 511 | ||
| 512 | 88572 | 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 | 88560 | (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 | 12 | if (CS%id_internal_heat_temp_tendency > 0) then |
| 520 | 0 | call post_data(CS%id_internal_heat_temp_tendency, dTdt_diag, CS%diag, alt_h=h) |
| 521 | endif | |
| 522 | 12 | if (CS%id_internal_heat_heat_tendency > 0) then |
| 523 | 0 | 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 | 0 | 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 | 0 | call post_data(CS%id_internal_heat_heat_tendency, dTdt_diag, CS%diag, alt_h=h) |
| 529 | endif | |
| 530 | 12 | if (CS%id_geothermal_buoyancy_flux > 0) then |
| 531 | 0 | 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 | ||
| 538 | end subroutine geothermal_in_place | |
| 539 | ||
| 540 | !> Initialize parameters and allocate memory associated with the geothermal heating module. | |
| 541 | 1 | subroutine 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 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 562 | ||
| 563 | 1 | CS%initialized = .true. |
| 564 | 1 | CS%diag => diag |
| 565 | 1 | CS%Time => Time |
| 566 | ||
| 567 | ! write parameters to the model log. | |
| 568 | 1 | 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 | 1 | units="W m-2 or various", default=0.0, scale=US%W_m2_to_QRZ_T) |
| 574 | 1 | CS%apply_geothermal = .not.(geo_scale == 0.0) |
| 575 | 1 | if (.not.CS%apply_geothermal) return |
| 576 | ||
| 577 | 8773 | 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 | 1 | "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 | 1 | 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 | 1 | do_not_log=((GV%nk_rho_varies<=0).or.(GV%nk_rho_varies>=GV%ke)) ) |
| 591 | 1 | if (CS%dRcv_dT_inplace >= 0.0) call MOM_error(FATAL, "geothermal_init: "//& |
| 592 | 0 | "GEOTHERMAL_DRHO_DT_INPLACE must be negative.") |
| 593 | ||
| 594 | 1 | if (len_trim(geo_file) >= 1) then |
| 595 | 0 | call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") |
| 596 | 0 | inputdir = slasher(inputdir) |
| 597 | 0 | filename = trim(inputdir)//trim(geo_file) |
| 598 | 0 | 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 | 0 | default="geo_heat") |
| 602 | 0 | call MOM_read_data(filename, trim(geotherm_var), CS%geo_heat, G%Domain) |
| 603 | 0 | do j=jsd,jed ; do i=isd,ied |
| 604 | 0 | CS%geo_heat(i,j) = (G%mask2dT(i,j) * geo_scale) * CS%geo_heat(i,j) |
| 605 | enddo ; enddo | |
| 606 | else | |
| 607 | 8773 | do j=jsd,jed ; do i=isd,ied |
| 608 | 8772 | CS%geo_heat(i,j) = G%mask2dT(i,j) * geo_scale |
| 609 | enddo ; enddo | |
| 610 | endif | |
| 611 | 1 | call pass_var(CS%geo_heat, G%domain) |
| 612 | ||
| 613 | 1 | 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 | 1 | x_cell_method='mean', y_cell_method='mean', area_cell_method='mean') |
| 622 | 1 | 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 | 1 | '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 | 1 | '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 | 1 | 'degC s-1', conversion=US%C_to_degC*US%s_to_T, v_extensive=.true.) |
| 637 | 1 | 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 | 0 | trim(thickness_units)//' s-1', conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.) |
| 643 | endif | |
| 644 | ||
| 645 | end subroutine geothermal_init | |
| 646 | ||
| 647 | !> Clean up and deallocate memory associated with the geothermal heating module. | |
| 648 | 1 | subroutine geothermal_end(CS) |
| 649 | type(geothermal_CS), intent(inout) :: CS !< Geothermal heating control struct | |
| 650 | 1 | if (allocated(CS%geo_heat)) deallocate(CS%geo_heat) |
| 651 | 1 | end 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 | ||
| 662 | 0 | end module MOM_geothermal |