Merge stemp branch

Squashed commit of the following:

commit 82ff748
Author: Joseph Zhang <[email protected]>
Date:   Wed Jan 3 15:11:59 2024 -0500

    More edits for OMP

commit 045d68f
Author: Joseph Zhang <[email protected]>
Date:   Wed Jan 3 15:10:34 2024 -0500

    Edit ZG's comments on sediment-water heat exchange

commit 5ed3103
Author: Zhengui Wang <[email protected]>
Date:   Tue Jan 2 12:41:02 2024 -0500

    change default stemp values

commit ad12e02
Author: Zhengui Wang <[email protected]>
Date:   Tue Jan 2 12:31:07 2024 -0500

    add scheme for the heat exchange between sediment layer and bottom water

sample_inputs/param.nml

@@ -592,6 +592,14 @@
   hmin_salt_ex = 0.2 ![m], no effect if isconsv=0 
   iprecip_off_bnd = 0 !if /=0, precip will be turned off near land bnd
 
+!-----------------------------------------------------------------------
+! Option to add a sediment layer for the buffer effect on temperature
+! stemp_stc=0 recovers previous results
+!-----------------------------------------------------------------------
+  stemp_stc  = 0.0   !heat transfer coefficient W.m-2.K-1
+  stemp_dz(1)= 1.0   !equivalent sediment buffer depth (m) for heat into sediment
+  stemp_dz(2)= 1.0   !equivalent sediment buffer depth (m) for heat out of sediment
+
 !-----------------------------------------------------------------------
 ! Turbulence closure.
 !-----------------------------------------------------------------------

src/Core/schism_glbl.F90

@@ -108,7 +108,8 @@ module schism_glbl
                       &prmsl_ref,hmin_radstress,eos_a,eos_b,eps1_tvd_imp,eps2_tvd_imp, &
                       &xlsc0,rearth_pole,rearth_eq,hvis_coef0,disch_coef(10),hw_depth,hw_ratio, &
                       &slr_rate,rho0,shw,gen_wsett,turbinj,turbinjds,alphaw,h1_bcc,h2_bcc,vclose_surf_frac, &
-                      &hmin_airsea_ex,hmin_salt_ex,shapiro0,loadtide_coef,h_massconsv,rinflation_icm
+                      &hmin_airsea_ex,hmin_salt_ex,shapiro0,loadtide_coef,h_massconsv,rinflation_icm, &
+                      &stemp_stc,stemp_dz(2)
 
   ! Misc. variables shared between routines
   integer,save :: nz_r,ieqstate,kr_co, &
@@ -267,6 +268,7 @@ module schism_glbl
   integer,save,allocatable :: iflux_e(:) !for computing fluxes
   integer,save,allocatable :: ielg2(:)      ! Local-to-global element index table (2-tier augmented)
   integer,save,allocatable :: iegl2(:,:)      ! Global-to-local element index table (2-tier augmented)
+  real(rkind),save,allocatable :: stemp(:)
 
   ! Node geometry data
   integer,save :: mnei  ! Max number of neighboring elements surrounding a node

src/Hydro/schism_init.F90

@@ -209,7 +209,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2, &
      &ielm_transport,max_subcyc,i_hmin_airsea_ex,hmin_airsea_ex,itransport_only,meth_sink, &
      &iloadtide,loadtide_coef,nu_sum_mult,i_hmin_salt_ex,hmin_salt_ex,h_massconsv,lev_tr_source, &
-     &rinflation_icm,iprecip_off_bnd,model_type_pahm
+     &rinflation_icm,iprecip_off_bnd,model_type_pahm,stemp_stc,stemp_dz
 
      namelist /SCHOUT/nc_out,iof_hydro,iof_wwm,iof_gen,iof_age,iof_sed,iof_eco,iof_icm_core, &
      &iof_icm_silica,iof_icm_zb,iof_icm_ph,iof_icm_cbp,iof_icm_sav,iof_icm_veg,iof_icm_sed, &
@@ -496,6 +496,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
       lev_tr_source=-9 !bottom
       iprecip_off_bnd=0
       model_type_pahm=10
+      stemp_stc=0; stemp_dz=1.0 !heat exchange between sediment and bottom water
 
       !Output elev, hvel by detault
       nc_out=1
@@ -1387,7 +1388,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
           &bdy_frc(ntracers,nvrt,nea),flx_sf(ntracers,nea),flx_bt(ntracers,nea), &
           &xlon_el(nea),ylat_el(nea),albedo(npa),flux_adv_vface(nvrt,ntracers,nea), &
           &wsett(ntracers,nvrt,nea),iwsett(ntracers),total_mass_error(ntracers), &
-          &iadjust_mass_consv(ntracers),wind_rotate_angle(npa),lev_tr_source2(ntracers),stat=istat)
+          &iadjust_mass_consv(ntracers),wind_rotate_angle(npa),lev_tr_source2(ntracers),stemp(nea),stat=istat)
       if(istat/=0) call parallel_abort('INIT: dynamical arrays allocation failure')
 !'
 
@@ -6132,6 +6133,11 @@ subroutine schism_init(iorder,indir,iths,ntime)
 !     end hot start section
       endif !ihot/=0
 
+      !Init sediment temp.
+      do i=1,nea
+         stemp(i)=tr_el(1,kbe(i)+1,i)
+      enddo
+
 ! MP from KM
 #ifdef USE_WWM
       ! Computation of the bed slope at nodes

src/Hydro/schism_step.F90

@@ -7677,7 +7677,7 @@ subroutine schism_step(it)
         if(istat/=0) call parallel_abort('STEP: fail to alloc (1.1)')
 
 !$OMP parallel default(shared) private(i,bigv,rat,j,jj,itmp1,itmp2,k,trnu,mm,swild,tmp,zrat, &
-!$OMP ta,ie,kin,swild_m,swild_w,tmp0,vnf,htot,top,dzz1)
+!$OMP ta,ie,kin,swild_m,swild_w,tmp0,vnf,htot,top,dzz1,tmp1)
 
 !       Point sources/sinks using operator splitting (that guarentees max.
 !       principle). Do nothing for net sinks
@@ -7708,6 +7708,25 @@ subroutine schism_step(it)
 
 !       Filter style horizontal diffusion scheme
 
+!       Heat exchange between sediment and bottom water
+        if(abs(stemp_stc)>1.d-16) then
+!$OMP     do
+          do i=1,nea
+             tmp1=(tr_el(1,kbe(i)+1,i)-stemp(i))*dt*stemp_stc !heat transfer budget (J.m-2)
+             if(tmp1>0) then !sediment temp. update; 4.184e6 is the heat capacity of water (J.m-3)
+               stemp(i)=stemp(i)+tmp1/(max(stemp_dz(1),1.d-2)*4.184d6)
+             else
+               stemp(i)=stemp(i)+tmp1/(max(stemp_dz(2),1.d-2)*4.184d6)  ! sediment temp. update
+             endif
+             !Bottom T update
+             tr_el(1,kbe(i)+1,i)=tr_el(1,kbe(i)+1,i)-tmp1/(max((ze(kbe(i)+1,i)-ze(kbe(i),i)),1.d-2)*4.184d6)
+             do k=1,kbe(i) 
+               tr_el(1,k,i)=tr_el(1,kbe(i)+1,i) 
+             enddo !k
+          enddo !i
+!$OMP     enddo
+        endif !abs(stemp_stc)
+
 !       Nudging: sum or product of horizontal & vertical relaxations 
 !$OMP   do 
         do i=1,nea