added diagnostic for internal heat in 3D

src/parameterizations/vertical/MOM_geothermal.F90

@@ -35,6 +35,8 @@ module MOM_geothermal
   type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                              !! regulate the timing of diagnostic output.
+  integer :: id_internal_heat_tend_3d = -1   !< ID for 3D diagnostic of internal heat
+
 end type geothermal_CS
 
 contains
@@ -100,6 +102,14 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
   real :: Irho_cp       ! inverse of heat capacity per unit layer volume
                         ! [degC H m2 J-1 ~> degC m3 J-1 or degC kg J-1]
 
+  real :: T_old         ! Temperature of each layer before any heat is added,
+                        ! for diagnostics [degC]
+  real, allocatable, dimension(:,:,:) :: h_old     ! Thickness of each layer before any heat is added,
+                                      ! for diagnostics [m or kg m-2]
+  real, allocatable, dimension(:,:,:) :: work_3d   ! Scratch variable used to calculate change in heat
+                                      ! due to geothermal
+  real :: Idt           ! inverse of the timestep [s-1]
+
   logical :: do_i(SZI_(G))
   integer :: i, j, k, is, ie, js, je, nz, k2, i2
   integer :: isj, iej, num_start, num_left, nkmb, k_tgt
@@ -119,6 +129,7 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
   Angstrom  = GV%Angstrom_H
   H_neglect = GV%H_subroundoff
   p_ref(:)  = tv%P_Ref
+  Idt       = 1.0 / dt
 
   if (.not.associated(tv%T)) call MOM_error(FATAL, "MOM geothermal: "//&
       "Geothermal heating can only be applied if T & S are state variables.")
@@ -136,6 +147,12 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
 !$OMP                                  wt_in_place,dTemp,dRcv,h_transfer,heating, &
 !$OMP                                  I_h)
 
+  ! Allocate diagnostic arrays if required
+  if (CS%id_internal_heat_tend_3d > 0) then
+    allocate(h_old(is:ie,js:je,nz)) ; h_old(:,:,:) = 0.0
+    allocate(work_3d(is:ie,js:je,nz)) ; work_3d(:,:,:) = 0.0
+  endif
+
   do j=js,je
     ! 1. Only work on columns that are being heated.
     ! 2. Find the deepest layer with any mass.
@@ -175,6 +192,12 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
     do k=nz,1,-1
       do i=isj,iej ; if (do_i(i)) then
 
+        ! Save temperature and thickness before any changes are made (for diagnostic)
+        if (CS%id_internal_heat_tend_3d > 0) then
+          T_old = tv%T(i,j,k)
+          h_old(i,j,k) = h(i,j,k)
+        endif
+
         if (h(i,j,k) > Angstrom) then
           if ((h(i,j,k)-Angstrom) >= h_geo_rem(i)) then
             h_heated = h_geo_rem(i)
@@ -294,6 +317,12 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
             ! endif
           endif
         endif
+
+        ! Calculate heat tendency due to addition and transfer of internal heat
+        if (CS%id_internal_heat_tend_3d > 0) then
+          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)
+        endif
+
       endif ; enddo
       if (num_left <= 0) exit
     enddo ! k-loop
@@ -304,6 +333,13 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
     enddo ; endif
   enddo ! j-loop
 
+  ! Post diagnostic of internal heat tendency in 3D
+  if (CS%id_internal_heat_tend_3d > 0) then
+    call post_data(CS%id_internal_heat_tend_3d, work_3d, CS%diag, alt_h = h_old)
+    deallocate(h_old)
+    deallocate(work_3d)
+  endif
+
 !  do i=is,ie ; do j=js,je
 !    resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_kg_m2 * &
 !           (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
@@ -392,6 +428,12 @@ subroutine geothermal_init(Time, G, param_file, diag, CS)
         x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
   if (id > 0) call post_data(id, CS%geo_heat, diag, .true.)
 
+  ! Diagnostic for tendency due to internal heat (in 3d)
+  CS%id_internal_heat_tend_3d=register_diag_field('ocean_model', &
+        'internal_heat_tend_3d', diag%axesTL, Time,              &
+        '3D heat tendency due to internal (geothermal) sources', &
+        'W m-2', v_extensive = .true.)
+
 end subroutine geothermal_init
 
 !> Clean up and deallocate memory associated with the geothermal heating module.