start working on Smagorinsky-like filter (ishapiro=2). Need test.

sample_inputs/param.nml

@@ -145,10 +145,10 @@
 !-----------------------------------------------------------------------
 ! 2nd stabilization method via Shapiro filter. This should normally be used 
 ! if indvel=ihorcon=0. To transition between eddying/non-eddying regimes, use
-! indvel=0, ihorcon/=0, and ishapiro=-1 (shapiro.gr3).
+! indvel=0, ihorcon/=0, and ishapiro=-1 (shapiro.gr3) or 2
 !-----------------------------------------------------------------------
-  ishapiro = 1 !on/off flag
-  shapiro0 = 0.5 !Shapiro filter strength, needed only if ishapiro=1; max is 0.5 
+  ishapiro = 1 !options
+  shapiro0 = 0.5 !Shapiro filter strength, needed only if ishapiro=1 (max is 0.5 ) or 2 (coefficient in tanh())
   niter_shap = 1 !needed if ishapiro/=0 - # of iterations with Shapiro filter. Suggested: 1
 
 !-----------------------------------------------------------------------
@@ -172,7 +172,7 @@
 ! force expression.
 !-----------------------------------------------------------------------
   icou_elfe_wwm = 0 
-  nstep_wwm = 1  !call WWM every this many time steps. If /=1, consider using quasi-steady mode in WWM
+  nstep_wwm = 1  !call WWM every this many time steps 
   iwbl = 0 !wave boundary layer formulation (used only if USE_WMM and 
            !icou_elfe_wwm/=0 and nchi=1. If icou_elfe_wwm=0, set iwbl=0): 
            !1-modified Grant-Madsen formulation; 2-Soulsby (1997)

sample_inputs/wwminput.nml.WW3

@@ -54,7 +54,7 @@
 /
 
 &INIT
- LHOTR          = F                  ! Use hotstart file (see &HOTFILE section)
+ LHOTR          = F                  ! Use hotstart file (see &HOTFILE section for input file)
  LINID          = T                  ! Initial condition; F for default; use T if using WW3 as i.c. etc
  INITSTYLE      = 2                  ! 1 - Parametric Jonswap, 2 - Read from Global NETCDF files, work only if IBOUNDFORMAT=2
 /
@@ -509,7 +509,7 @@
 
 &HOTFILE
  LHOTF          = F                  ! Write hotfile
- FILEHOT_OUT    = 'wwm_hot_out'      !'.nc' suffix will be added 
+ FILEHOT_OUT    = 'hotfile_out.nc'      !name of output
  BEGTC          = '20030908.000000'  !Starting time of hotfile writing. With ihot!=0 in SCHISM,
                                      !this will be whatever the new hotstarted time is (even with ihot=2)
  DELTC          = 86400.             ! time between hotfile writes
@@ -526,11 +526,11 @@
                                      !  MPI_REDUCE is then used and thus you'd avoid too freq. output 
                                      ! 1: hotfiles in separate files, each associated
                                      ! with one process
- FILEHOT_IN     = 'wwm_hot_in.nc'    ! (Full) hot file name for input
+ FILEHOT_IN     = 'hotfile_in.nc'    ! (Full) hot file name for input (which can be FILEHOT_OUT above)
  HOTSTYLE_IN    = 2                  ! 1: binary hotfile of data as input
                                      ! 2: netcdf hotfile of data as input (default)
  IHOTPOS_IN     = 1                  ! Position in hotfile (only for netcdf)
-                                     ! for reading
+                                     ! for reading. If LCYCLEHOT=T, this can be 1 or 2
  MULTIPLEIN     = 0                  ! 0: read hotfile from one single file
                                      ! 1: read hotfile from multiple files (must use same # of CPU?)
 /

sample_inputs/wwminput.nml.spectra

@@ -15,10 +15,10 @@
  LMONO_OUT      = F
  LNAUTOUT       = T                  ! Nautical output of all quantities in degrees
  BEGTC          = '19980901.000000'  ! Time for start the simulation, ex:yyyymmdd. hhmmss
- DELTC          = 5                 ! Time step 
+ DELTC          = 5                 ! Time step (MUST match dt*nstep_wwm in SCHISM!)
  UNITC          = 'SEC'              ! Unity of time step
  ENDTC          = '19980901.010000'  ! Time for stop the simulation, ex:yyyymmdd. hhmmss
- DMIN           = 0.0001              ! Minimum water depth. This is not used in selfe; with selfe this is set automatically to h0 in param.in
+ DMIN           = 0.01              ! Minimum water depth. Use 'h0' as in SCHISM
 /
 
 &COUPL
@@ -52,7 +52,7 @@
 /
 
 &INIT
- LHOTR          = F                  ! Use hotstart file
+ LHOTR          = F                  ! Use hotstart file. If T, the name is specified in &HOTFILE
  LINID          = F                  ! Initial condition; F for default; use T if using WW3 as i.c. etc
  INITSTYLE      = 1                  ! 1 - Parametric Jonswap, 2 - Read from Global NETCDF files, work only if IBOUNDFORMAT=3
 /
@@ -429,14 +429,14 @@ ASPAR_LOCAL_LEVEL = 0               ! locality level 0 - a lot of memory 10 - no
                      !    MPI_REDUCE is then used.
                      ! 1: hotfiles in separate files, each associated
                      !    with one process
- FILEHOT_OUT    = 'hotfile_out.dat'
+ FILEHOT_OUT    = 'hotfile_out.nc' !name of hot outputs
  HOTSTYLE_IN    = 2  ! 1: binary hotfile of data as input
                      ! 2: netcdf hotfile of data as input (default)
  IHOTPOS_IN        = 1 ! Position in hotfile (only for netcdf)
-                       ! for reading
+                       ! for reading. If LCYCLEHOT=T, this can be 1 or 2
  MULTIPLEIN     = 0    ! 0: read hotfile from one single file
                        ! 1: read hotfile from multiple files.
- FILEHOT_IN     = 'hotfile_in.dat' ! Hot file name for input
+ FILEHOT_IN     = 'hotfile_in.nc' ! Hot file name for input if LHOTR=T in &INIT. This can be simply hotfile_out.nc above
 /
 
 &NESTING

src/Core/schism_glbl.F90

@@ -105,7 +105,7 @@ 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,h1_bcc,h2_bcc,vclose_surf_frac, &
-                      &hmin_airsea_ex
+                      &hmin_airsea_ex,shapiro0
 
   ! Misc. variables shared between routines
   integer,save :: nz_r,ieqstate,kr_co, &
@@ -250,9 +250,9 @@ module schism_glbl
   !where j is the axis id, i is the component id, ie is the local element id (aug.)
   !Undefined for ics=1
   real(rkind),save,allocatable :: eframe(:,:,:)
-  !x,y coordinates of each element node in the _element_ frame
-  !xel(4,nea)
-  real(rkind),save,allocatable :: xel(:,:),yel(:,:)
+  !x,y coordinates of each element node/side in the _element_ frame
+  !xel(4,nea), xs_el(4,nea)
+  real(rkind),save,allocatable :: xel(:,:),yel(:,:),xs_el(:,:),ys_el(:,:)
   !shape_c2(4,2,nea)- for quads only. The shape function values for the mid-pts of 2 diagnonals inside
   ! the quad formed by 4 sidecenters
   real(rkind),save,allocatable :: shape_c2(:,:,:)
@@ -325,8 +325,7 @@ module schism_glbl
   real(rkind),save,allocatable :: zs(:,:)         ! z-coord. (local frame - vertical up)
   !Transformation tensor for side frame: sframe(i,j,isd)
   ! where j is the axis id, i is the component id, isd is the local side id (aug.)
-  ! For ics=1, only sframe(1:2,1:2,isd) are used
-  ! x-axis is from adjacent elem 1 to 2.
+  ! For ics=1, only sframe(1:2,1:2,isd) are used. In all cases, x-axis is from adjacent elem 1 to 2.
   real(rkind),save,allocatable :: sframe(:,:,:)
   !cos/sin of side normals. If ics=1, these are same as sframe(1:2,1,isd)
   !If ics=2, these are product of sframe(:,1,:) and local lat/lon frame
@@ -409,7 +408,8 @@ module schism_glbl
   real(rkind),save,allocatable :: total_mass_error(:) !(ntracers) Total mass error after advection step for mass correction
   !x & y-component of velocity at side centers & whole levels
   !For ics=1, these are defined in the global frame
-  !For ics=2, these are defined in the ll frame
+  !For ics=2, these are defined in the ll frame (not local side frame). Note
+  !that these are not defined on sframe()
   real(rkind),save,allocatable :: su2(:,:),sv2(:,:) 
   !velocity at nodes in an element. In ll if ics=2
   !ufg(1:4,1:nvrt,1:nea)

src/Hydro/schism_init.F90

@@ -152,7 +152,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
                     &ubl1,ubl2,ubl3,ubl4,ubl5,ubl6,ubl7,ubl8,xn1, &
                     &xn2,yn1,yn2,xstal,ystal,ae,THAS,THAF,err_max,rr,suma, &
                     &te,sa,wx1,wx2,wy1,wy2,aux1,aux2,time,ttt, &
-                    &et,qq,tr,ft1,dep,wtratio,shapiro0,rmaxvel,sav_cd0
+                    &et,qq,tr,ft1,dep,wtratio,rmaxvel,sav_cd0
 
 
 #ifdef USE_FIB
@@ -571,26 +571,31 @@ subroutine schism_init(iorder,indir,iths,ntime)
 
 !...  Shapiro filter 
 !      call get_param('param.in','ishapiro',1,ishapiro,tmp,stringvalue)
-      if(iabs(ishapiro)>1) then
+      if(ishapiro<-1.or.ishapiro>2) then
         write(errmsg,*)'Illegal ishapiro:',ishapiro
         call parallel_abort(errmsg)
       endif
 
       if(ishapiro==1) then
-!        call get_param('param.in','shapiro',2,itmp,shapiro0,stringvalue)   
         if(shapiro0<0._rkind.or.shapiro0>0.5_rkind) then
           write(errmsg,*)'Illegal shapiro:',shapiro0
           call parallel_abort(errmsg)
         endif
       endif !ishapiro==1
 
+      if(ishapiro==2) then
+        if(shapiro0<0._rkind) then
+          write(errmsg,*)'Illegal shapiro(2):',shapiro0
+          call parallel_abort(errmsg)
+        endif
+      endif !ishapiro
+
       if(ishapiro/=0) then
-!        call get_param('param.in','niter_shap',1,niter_shap,tmp,stringvalue)
         if(niter_shap<0) then
           write(errmsg,*)'Illegal niter_shap:',niter_shap
           call parallel_abort(errmsg)
         endif
-      endif !ishapiro==1
+      endif !ishapiro
 
 !...  Horizontal diffusivity option; only works for upwind/TVD
 !     ihdif=0 means all hdif=0 and no hdif.gr3 is needed
@@ -1656,7 +1661,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &kbp00(npa),kbp_e(np),idry(npa),hmod(npa),znl(nvrt,npa), &
      &kbs(nsa),idry_s(nsa),isidenei2(4,ns),zs(nvrt,nsa), &
      &delj(ns),ibnd_ext_int(npa),pframe(3,3,npa),sigma_lcl(nvrt,npa),shape_c2(4,2,nea), &
-     &snx(nsa),sny(nsa),stat=istat)
+     &snx(nsa),sny(nsa),xs_el(4,nea),ys_el(4,nea),stat=istat)
       if(istat/=0) call parallel_abort('INIT: grid geometry arrays allocation failure')
 !'
 
@@ -1696,7 +1701,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
          &  pr2(npa),airt2(npa),shum2(npa),pr(npa),sflux(npa),srad(npa),tauxz(npa),tauyz(npa), &
          &  fluxsu(npa),fluxlu(npa),hradu(npa),hradd(npa),cori(nsa),Cd(nsa), &
          &  Cdp(npa),rmanning(npa),rough_p(npa),dfv(nvrt,npa),elev_nudge(npa),uv_nudge(npa), &
-         &  hdif(nvrt,npa),shapiro(nsa),fluxprc(npa),fluxevp(npa), & 
+         &  hdif(nvrt,npa),shapiro(ns),fluxprc(npa),fluxevp(npa), & 
          &  sparsem(0:mnei_p,np), & !sparsem for non-ghosts only
          &  tr_nudge(natrm,npa), & 
          &  fun_lat(0:2,npa),dav(2,npa),elevmax(npa),dav_max(2,npa),dav_maxmag(npa), &
@@ -2067,6 +2072,16 @@ subroutine schism_init(iorder,indir,iths,ntime)
 !      endif !ics
 
 !...  Finish off some remaining geometric calcualtions
+!     Sidecenter coord in the elem frame
+      do i=1,nea
+        do j=1,i34(i)
+          j1=nxq(1,j,i34(i))
+          j2=nxq(2,j,i34(i))
+          xs_el(j,i)=0.5d0*(xel(j1,i)+xel(j2,i))
+          ys_el(j,i)=0.5d0*(yel(j1,i)+yel(j2,i))
+        enddo !j
+      enddo !i
+
 !      !weno>
       if (itr_met==4) then !WENO
         call set_isbe !identify boundary elements
@@ -3051,7 +3066,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
       slam0=slam0*pi/180.d0
       sfea0=sfea0*pi/180.d0
 
-!...  Read in shaprio.gr3
+!...  Set shapiro(:). For ishapiro=2, this will be done in _step
       if(ishapiro==1) then
         shapiro(:)=shapiro0
       else if(ishapiro==-1) then
@@ -3065,7 +3080,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
           if(ipgl(i)%rank==myrank) swild(ipgl(i)%id)=tmp
         enddo !i
         close(32)
-        do i=1,nsa
+        do i=1,ns !a
           shapiro(i)=sum(swild(isidenode(1:2,i)))/2.d0
           !Check range
           if(shapiro(i)<0.d0.or.shapiro(i)>0.5d0) call parallel_abort('INIT: check shapiro')

src/Hydro/schism_step.F90

@@ -139,7 +139,7 @@ subroutine schism_step(it)
                  &ibnd,jfr,ncyc,iter,nlev,klev,kin,nqdim,limit,jmin, &
                  &ipsgb,iadvf,ifl_bnd,nnel,jlev,ndelt_min,ndelt_max,ii,id, &
                  &id2,id3,ip,ndelt,ibot_fl,ibelow,indx,nj,ind,ind2,lim,in1, &
-                 &in2,irank_s,itmp,itmp1,itmp2,node1,node2,ndim,mk,nd_lam, &
+                 &in2,in3,irank_s,itmp,itmp1,itmp2,node1,node2,ndim,mk,nd_lam, &
                  &iee,idel,irow,icol,ieq,ij,kbb,lwrite,lit,ihot_len,IHOTSTP, &
                  &itmpf,ibt,mmk,ndo,n,ind_tr,n_columns,ncid_hot,node_dim,elem_dim, &
                  &side_dim,nvrt_dim,ntracers_dim,three_dim,two_dim,one_dim, &
@@ -3257,24 +3257,8 @@ subroutine schism_step(it)
                   call vinter(1,nvrt,1,zs(k,j),kbs(jsj),nvrt,k,zs(:,jsj),sv2(:,jsj),swild(2),ibelow)
                 endif !isd/=i
                 swild10(i,1)=swild(1); swild10(i,2)=swild(2)
-
-!                !Project to side and elem. frame for ics=2
-!                if(ics==1) then
-!                  soln(i,1)=swild(1); soln(i,2)=swild(2)
-!                  !Elem frame
-!                  swild10(1,i)=swild(1); swild10(2,i)=swild(2)
-!                else
-!                  call project_hvec(swild(1),swild(2),sframe(:,:,jsj),sframe(:,:,j),soln(i,1),soln(i,2))
-!                  call project_hvec(swild(1),swild(2),sframe(:,:,jsj),eframe(:,:,ie),swild10(1,i),swild10(2,i))
-!                endif
               enddo !i=1,i34(ie)
 
-              !\nabla{u} const. within an elem. (elem. frame)
-!              dudx=dot_product(swild10(1,1:i34(ie)),dldxy_sd(1:i34(ie),1,ie))
-!              dudy=dot_product(swild10(1,1:i34(ie)),dldxy_sd(1:i34(ie),2,ie))
-!              dvdx=dot_product(swild10(2,1:i34(ie)),dldxy_sd(1:i34(ie),1,ie))
-!              dvdy=dot_product(swild10(2,1:i34(ie)),dldxy_sd(1:i34(ie),2,ie))
-
               !do i=1,2 !i34(ie) !2 sides per elem.
                 !jsj=elside(i,ie)
                 !if(isbs(jsj)==-1) then !deal with land bnd
@@ -3406,6 +3390,94 @@ subroutine schism_step(it)
         deallocate(swild98)
       endif !ihorcon/=0
 
+!...  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   workshare
+        shapiro=0.d0
+!$OMP   end workshare
+
+!$OMP   do
+        do j=1,ns !residents only
+!          if(isdel(2,j)==0.or.idry_s(j)==1) cycle
+!          if(idry_e(isdel(1,j))==1.or.idry_e(isdel(2,j))==1) cycle
+          if(idry_s(j)==1) cycle
+          if(ihydraulics/=0.and.nhtblocks>0) then
+            if(isblock_sd(1,j)/=0) cycle
+          endif
+
+          !wet side 
+          vmax=0.d0 !init max gradient
+          do k=kbs(j)+1,nvrt !strength= 0 at bottom
+            do l=1,2 !element
+              ie=isdel(l,j)
+              if(ie<=0) cycle
+              if(idry_e(ie)==1) cycle
+
+              !Wet elem
+              do i=1,i34(ie) !prep. side vel. via vertical interp
+                jsj=elside(i,ie)
+                if(jsj==j) then
+                  swild(1)=su2(k,j)
+                  swild(2)=sv2(k,j)
+                else
+                  call vinter(1,nvrt,1,zs(k,j),kbs(jsj),nvrt,k,zs(:,jsj),su2(:,jsj),swild(1),ibelow)
+                  call vinter(1,nvrt,1,zs(k,j),kbs(jsj),nvrt,k,zs(:,jsj),sv2(:,jsj),swild(2),ibelow)
+                endif !isd/=i
+                !in ll frame if ics=2
+                swild10(i,1)=swild(1); swild10(i,2)=swild(2)
+              enddo !i=1,i34(ie)
+
+              !Reconstruct local gradient
+              ll=lindex_s(j,ie)
+              if(ll==0) then
+                write(errmsg,*)'STEP: Cannot find a side(2)'
+                call parallel_abort(errmsg)
+              endif
+              in1=nxq(1,ll,i34(ie))
+              in2=nxq(i34(ie)-1,ll,i34(ie))
+
+              !2x2 matrix
+              if(i34(ie)==3) then
+                swild4(1,1)=xs_el(in1,ie)-xs_el(ll,ie)
+                swild4(1,2)=ys_el(in1,ie)-ys_el(ll,ie)
+                swild4(2,1)=xs_el(in2,ie)-xs_el(ll,ie)
+                swild4(2,2)=ys_el(in2,ie)-ys_el(ll,ie)
+                !RHS; 2nd index is (u,v)
+                swild2(1,1:2)=swild10(in1,1:2)-swild10(ll,1:2)
+                swild2(2,1:2)=swild10(in2,1:2)-swild10(ll,1:2)
+              else !quad
+                in3=nxq(2,ll,i34(ie))
+                swild4(1,1)=xs_el(in3,ie)-xs_el(ll,ie)
+                swild4(1,2)=ys_el(in3,ie)-ys_el(ll,ie)
+                swild4(2,1)=xs_el(in2,ie)-xs_el(in1,ie)
+                swild4(2,2)=ys_el(in2,ie)-ys_el(in1,ie)
+                swild2(1,1:2)=swild10(in3,1:2)-swild10(ll,1:2)
+                swild2(2,1:2)=swild10(in2,1:2)-swild10(in1,1:2)
+              endif !i34(ie)
+              delta_wc=swild4(1,1)*swild4(2,2)-swild4(1,2)*swild4(2,1)
+              if(delta_wc==0.d0) then
+                write(errmsg,*)'STEP: delta_wc=0:',delta_wc,ielg(ie)
+                call parallel_abort(errmsg)
+              endif
+
+              dudx=(swild2(1,1)*swild4(2,2)-swild2(2,1)*swild4(1,2))/delta_wc
+              dudy=(swild2(2,1)*swild4(1,1)-swild2(1,1)*swild4(2,1))/delta_wc
+              dvdx=(swild2(1,2)*swild4(2,2)-swild2(2,2)*swild4(1,2))/delta_wc
+              dvdy=(swild2(2,2)*swild4(1,1)-swild2(1,2)*swild4(2,1))/delta_wc
+
+              !Original Smag; Griffiths used a different one
+              vmax=max(vmax,sqrt(dudx*dudx+dvdy*dvdy+0.5d0*(dudy+dvdx)**2))
+            enddo !l=1,2
+          enddo !k=kbs(j)+1,nvrt 
+
+          shapiro(j)=0.5d0*tanh(vmax*shapiro0)
+        enddo !j=1,ns
+!$OMP   end do
+      endif !ishapiro==2
+
       if(myrank==0) write(16,*)'done hvis... '
 
 #ifdef TIMER2

src/Readme.beta_notes

@@ -2,6 +2,7 @@ Changes in param.[in,nml] and main parameter inputs for other modules:
 
 Github versions:
 
+(): added ishapiro=2 - Smagorinsky like filter. In this option, shapiro0 is the coefficient.
 b43afea (April 2, 2021): changed x,y to double in nc outputs (for newer visIT);
 9150d86 (Mar 31, 2021): added a new SAL option (iloadtide=2) using scaling;
 084e149 (Feb 25, 2021): Removed parameters: ibtrack_openbnd, iwindoff, dzb_decay (with hardwire); also