Added an optional to impose min shapiro with ishapiro=2.

src/Core/schism_glbl.F90

@@ -435,6 +435,7 @@ module schism_glbl
   real(rkind),save,allocatable :: rho_mean(:,:)         ! mean density
   real(rkind),save,allocatable :: Cdp(:)         ! drag at node
   real(rkind),save,allocatable :: rmanning(:)         ! Manning's n at node
+  real(rkind),save,allocatable :: shapiro_min(:)      !min of Shapiro filter strength (used with some ishapiro options)
   real(rkind),save,allocatable,target :: windx(:),windy(:) !wind vector
   real(rkind),save,allocatable,target :: sdbt(:,:,:),shapiro(:), &
                                   &windx1(:),windy1(:),windx2(:),windy2(:), &
@@ -451,7 +452,7 @@ module schism_glbl
                                   &tr_nudge(:,:),fun_lat(:,:), &
                                   &elev_nudge(:),uv_nudge(:),fluxprc(:),fluxevp(:), &
                                   &dav(:,:),elevmax(:),dav_max(:,:),dav_maxmag(:), & 
-                                  &etaic(:),diffmax(:),diffmin(:),dfq1(:,:),dfq2(:,:) 
+                                  &etaic(:),diffmax(:),diffmin(:),dfq1(:,:),dfq2(:,:)
 
   !(2,npa). ocean-ice stress (junk if no ice) [m^2/s/s]
   real(rkind),save,allocatable :: tau_oi(:,:)

src/Hydro/schism_init.F90

@@ -1340,7 +1340,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
          &  diffmax(npa),diffmin(npa),dfq1(nvrt,npa),dfq2(nvrt,npa), & 
          &  iwater_type(npa),rho_mean(nvrt,nea),erho(nvrt,nea),& 
          & surf_t1(npa),surf_t2(npa),surf_t(npa),etaic(npa),sav_alpha(npa), &
-         & sav_h(npa),sav_nv(npa),sav_di(npa),sav_cd(npa),stat=istat)
+         & sav_h(npa),sav_nv(npa),sav_di(npa),sav_cd(npa),shapiro_min(npa),stat=istat)
       if(istat/=0) call parallel_abort('INIT: other allocation failure')
 
 !     Tracers
@@ -2721,6 +2721,22 @@ subroutine schism_init(iorder,indir,iths,ntime)
           if(shapiro(i)<0.d0.or.shapiro(i)>0.5d0) call parallel_abort('INIT: check shapiro')
 !'
         enddo !i
+      else if(ishapiro==2) then !read in optional shapiro_min.gr3
+        shapiro_min=0.d0 !init min in case shapiro_min.gr3 does not exist
+        inquire(file=in_dir(1:len_in_dir)//'shapiro_min.gr3', exist=lexist)
+        if(lexist) then
+          open(32,file=in_dir(1:len_in_dir)//'shapiro_min.gr3',status='old')
+          read(32,*)
+          read(32,*) itmp1,itmp2
+          if(itmp1/=ne_global.or.itmp2/=np_global) &
+     &call parallel_abort('Check shapiro_min.gr3')
+          do i=1,np_global
+            read(32,*)j,xtmp,ytmp,tmp
+            if(tmp<0.d0.or.tmp>0.5d0) call parallel_abort('INIT: check shapiro_min')
+            if(ipgl(i)%rank==myrank) shapiro_min(ipgl(i)%id)=tmp
+          enddo !i
+          close(32)
+        endif !lexist
       endif !ishapiro==-1
 
 !...  Horizontal viscosity option

src/Hydro/schism_step.F90

@@ -3421,7 +3421,7 @@ subroutine schism_step(it)
 !...  ishapiro=2: Smag-like filter
       if(ishapiro==2) then
 !$OMP   parallel default(shared) private(j,k,l,ie,i,jsj,swild,ibelow,swild10,ll, &
-!$OMP   in1,in2,in3,swild2,swild4,delta_wc,vmax,dudx,dudy,dvdx,dvdy)
+!$OMP   in1,in2,in3,swild2,swild4,delta_wc,vmax,vmin,dudx,dudy,dvdx,dvdy)
 
 !$OMP   workshare
         shapiro=0.d0
@@ -3502,7 +3502,9 @@ subroutine schism_step(it)
           enddo !k=kbs(j)+1,nvrt 
 
           shapiro(j)=0.5d0*tanh(dt*vmax*shapiro0)
-
+          !min value
+          vmin=0.5d0*(shapiro_min(isidenode(1,j))+shapiro_min(isidenode(2,j)))
+          shapiro(j)=max(shapiro(j),vmin)
         enddo !j=1,ns
 !$OMP   end do
 !$OMP   end parallel
@@ -3519,6 +3521,7 @@ subroutine schism_step(it)
             if(isdel(2,j)==0) then !isidenei2 not defined
               bcc(1,1,j)=shapiro(j)
             else
+              !Weighted average so positivity is guaranteed
               bcc(1,1,j)=shapiro(j)+0.5d0/4.d0*(sum(shapiro(isidenei2(1:4,j)))-4.d0*shapiro(j)) 
             endif
           enddo !j=1,ns