More work on iloadtide=3 option (Stepanov & Hughes 2004)

sample_inputs/param.nml

@@ -60,15 +60,18 @@
 
 !-----------------------------------------------------------------------
 ! Option to add self-attrection loading tide (SAL) into tidal potential 
-! (usually for large-scale model). 
+! (usually for basin-scale model). 
+! If iloadtide=0, no SAL.
 ! If iloadtide=1, needs inputs: loadtide_[FREQ].gr3, 
 ! where [FREQ] are freq names (shared with tidal potential, in upper cases) 
 ! and the _two_ 'depths' inside are amplitude (m) and phases (degrees behind GMT), 
 ! interpolated from global tide model (e.g. FES2014). In this option, SAL is 
 ! lumped into tidal potential so it shares some parameters with tidal potential 
 ! in bctides.in (cut-off depth, frequencies).
-! If iloadtide=2, use a simple scaling for gravity approach (in this option,
+! If iloadtide=2 or 3, use a simple scaling for gravity approach (in this option,
 ! SAL is applied everywhere and does not share parameters with tidal potential).
+! For iloadtide=2, a const scaling (0.1) is used; for iloadtide=3, the scaling is 
+! dependent on depth (Stepanov & Hughes 2004).
 !-----------------------------------------------------------------------
   iloadtide = 0
 

src/Hydro/schism_step.F90

@@ -285,16 +285,23 @@ subroutine schism_step(it)
 #endif
 
 !     SAL option: scale gravity
-      if(iloadtide==2) then
+      if(iloadtide==2) then !simple const
+        grav2=grav*0.9d0
+      else if(iloadtide==3) then !Stepanov & Hughes (2004)
         do i=1,npa
           dp1=max(dp(i),0.d0)
           tmp1=sqrt(dp1)
           beta_bar=-9.8169d-3+1.8289d-3*tmp1+4.3787d-4*dp1-2.9042d-5*dp1*tmp1+ &
      &6.6038d-7*dp1*dp1-4.7393d-9*dp1*dp1*tmp1-1.9354d-11*dp1*dp1*dp1+ &
      &2.6969d-13*dp1*dp1*dp1*tmp1
-          grav2(i)=grav*min(0.15d0,beta_bar) !0.9d0
+          beta_bar=max(0.d0,min(0.12d0,beta_bar))
+
+          !Debug
+          !if(it==iths_main+1) write(12,*)'SAL beta=',i,dp(i),beta_bar
+
+          grav2(i)=grav*(1.d0-beta_bar) !0.9d0
         enddo !i
-      else
+      else !iloadtide=0,1
         grav2=grav
       endif