Changed horizontal diffusion method to filter. Tested.

Squashed commit of the following: commit 402928b Author: Joseph Zhang <[email protected]> Date: Tue Aug 6 19:46:20 2024 -0500 Add # of iteration for hdif; tested commit 7442629 Author: Joseph Zhang <[email protected]> Date: Mon Aug 5 14:54:02 2024 -0400 Change horizontal diffusion to a (conservative) filter
schism-dev · Aug 7, 2024 · 2f88fc5 · 2f88fc5
1 parent d04e0f4
commit 2f88fc5
Show file tree

Hide file tree

Showing 4 changed files with 88 additions and 34 deletions.
diff --git a/sample_inputs/param.nml b/sample_inputs/param.nml
@@ -328,9 +328,12 @@
   level_age = 9, -999 !default: -999 (all levels)
 
 !-----------------------------------------------------------------------
-! Horizontal diffusivity option. if ihdif=1, horizontal diffusivity is given in hdif.gr3
+! Horizontal diffusivity option. if ihdif=1, horizontal diffusivity in the 
+! form of a filter is given in hdif.gr3 (max strength should be ~1/8);
+! and # of iterations for filter is given by niter_hdif.
 !-----------------------------------------------------------------------
   ihdif = 0 
+  niter_hdif = 1 !used if ihdif/=0
 
 !-----------------------------------------------------------------------
 ! Bottom friction. 

diff --git a/src/Core/schism_glbl.F90 b/src/Core/schism_glbl.F90
@@ -94,7 +94,7 @@ module schism_glbl
                   &nstep_ice,niter_shap,iunder_deep,flag_fib,ielm_transport,max_subcyc, &
                   &itransport_only,iloadtide,nc_out,nu_sum_mult,iprecip_off_bnd, &
                   &iof_ugrid,model_type_pahm,iof_icm_sav,iof_icm_marsh,iof_icm_sfm,iof_icm_ba,&
-                  &iof_icm_clam,nbins_veg_vert,veg_lai,veg_cw
+                  &iof_icm_clam,nbins_veg_vert,veg_lai,veg_cw,niter_hdif
   integer,save :: ntrs(natrm),nnu_pts(natrm),mnu_pts,lev_tr_source(natrm)
   integer,save,dimension(:),allocatable :: iof_hydro,iof_wwm,iof_gen,iof_age,iof_sed,iof_eco, &
      &iof_icm,iof_icm_core,iof_icm_silica,iof_icm_zb,iof_icm_ph,iof_icm_srm,iof_cos,iof_fib, &
@@ -414,7 +414,7 @@ module schism_glbl
   real(rkind),save,allocatable :: bdy_frc(:,:,:) !body force at prism center Q_{i,k}
   real(rkind),save,allocatable :: flx_sf(:,:) !surface b.c. \kappa*dC/dz = flx_sf (at element center)
   real(rkind),save,allocatable :: flx_bt(:,:) !bottom b.c.
-  real(rkind),save,allocatable :: hdif(:,:) !horizontal diffusivity
+  real(rkind),save,allocatable :: hdif(:) !horizontal diffusivity
   real(rkind),save,allocatable :: tr_nd0(:,:,:) ! Initial tracer conc. at nodes
   real(rkind),save,allocatable :: rkai_num(:,:,:) !DVD (numerical mixing) [C^2]/sec
   real(rkind),save,allocatable :: eta1(:)   ! Elevation at nodes at previous timestep

diff --git a/src/Hydro/schism_init.F90 b/src/Hydro/schism_init.F90
@@ -209,7 +209,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &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,stemp_stc,stemp_dz, &
      &veg_vert_z,veg_vert_scale_cd,veg_vert_scale_N,veg_vert_scale_D,veg_lai,veg_cw, &
-     &RADFLAG
+     &RADFLAG,niter_hdif
 
      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_srm,iof_icm_sav,iof_icm_marsh,iof_icm_sfm, &
@@ -499,6 +499,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
       model_type_pahm=10
       stemp_stc=0; stemp_dz=1.0 !heat exchange between sediment and bottom water
       RADFLAG='LON' !if WWM is used, this will be overwritten
+      niter_hdif=1
       veg_vert_z=(/((i-1)*0.4d0,i=1,nbins_veg_vert+1)/) ![m]
       veg_vert_scale_cd=(/(1.0d0,i=1,nbins_veg_vert+1)/) !scaling [-]
       veg_vert_scale_N=(/(1.0d0,i=1,nbins_veg_vert+1)/)
@@ -1389,7 +1390,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),shapiro_smag(nsa),fluxprc(npa),fluxevp(npa),prec_snow(npa),prec_rain(npa), & 
+         & hdif(npa),shapiro(nsa),shapiro_smag(nsa),fluxprc(npa),fluxevp(npa),prec_snow(npa),prec_rain(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), &
@@ -3016,7 +3017,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &call parallel_abort('Check hdif.gr3')
           do i=1,np_global
             read(32,*)j,xtmp,ytmp,tmp 
-            if(tmp<0.d0) then
+            if(tmp<0.d0.or.tmp>0.2d0) then
               write(errmsg,*)'hdif out of bound:',tmp,i
               call parallel_abort(errmsg)
             endif
@@ -3027,7 +3028,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
         call mpi_bcast(buf3,ns_global,rtype,0,comm,istat)
 
         do i=1,np_global
-          if(ipgl(i)%rank==myrank) hdif(:,ipgl(i)%id)=buf3(i) !tmp
+          if(ipgl(i)%rank==myrank) hdif(ipgl(i)%id)=buf3(i) !tmp
         enddo !i
       endif !ihdif/=0
 

diff --git a/src/Hydro/transport_TVD_imp.F90 b/src/Hydro/transport_TVD_imp.F90
@@ -1892,7 +1892,7 @@ subroutine do_transport_tvd_imp(it,ntr,difnum_max_l) !,nvrt1,npa1,dfh1)
       difnum_max_l=0.d0 !max. diffusion number reached by this process (check stability)
 !$OMP end single
 
-!$OMP do reduction(max: difnum_max_l) 
+!$OMP do 
       do i=1,ne
         if(idry_e(i)==1) cycle
 
@@ -1983,29 +1983,28 @@ subroutine do_transport_tvd_imp(it,ntr,difnum_max_l) !,nvrt1,npa1,dfh1)
             !Body source
             rrhs(1,kin)=rrhs(1,kin)+dt*bdy_frc(m,k,i)
 
-            !Horizontal diffusion
-            if(ihdif/=0) then
-              do j=1,i34(i) !sides
-                jsj=elside(j,i) !residents
-                iel=ic3(j,i)
-                if(iel==0.or.idry_e(max(1,iel))==1) cycle
-
-                nd1=isidenode(1,jsj)
-                nd2=isidenode(2,jsj)
-                hdif_tmp=(hdif(k,nd1)+hdif(k,nd2)+hdif(k-1,nd1)+hdif(k-1,nd2))/4.d0
-                if(k>=kbs(jsj)+1) then
-                  !av_h=(znl(k,nd1)-znl(k-1,nd1)+znl(k,nd2)-znl(k-1,nd2))/2.d0
-                  !!average height
-                  av_h=zs(k,jsj)-zs(k-1,jsj)
-                  if(av_h<=0.d0) call parallel_abort('TRAN_IMP: av_h<=0')
-                  !Check diffusion number; write warning message
-                  difnum=dt_by_bigv*hdif_tmp/delj(jsj)*av_h*distj(jsj)
-!                  if(difnum>difnum_max_l) difnum_max_l=difnum
-                  difnum_max_l=max(difnum_max_l,difnum)
-                  rrhs(1,kin)=rrhs(1,kin)+difnum*(trel_tmp(m,k,iel)-trel_tmp(m,k,i))
-                endif !k>=
-              enddo !j
-            endif !ihdif/=0
+!            !Horizontal diffusion
+!            if(ihdif/=0) then
+!              do j=1,i34(i) !sides
+!                jsj=elside(j,i) !residents
+!                iel=ic3(j,i)
+!                if(iel==0.or.idry_e(max(1,iel))==1) cycle
+! 
+!                nd1=isidenode(1,jsj)
+!                nd2=isidenode(2,jsj)
+!                hdif_tmp=(hdif(k,nd1)+hdif(k,nd2)+hdif(k-1,nd1)+hdif(k-1,nd2))/4.d0
+!                if(k>=kbs(jsj)+1) then
+!                  !av_h=(znl(k,nd1)-znl(k-1,nd1)+znl(k,nd2)-znl(k-1,nd2))/2.d0
+!                  !!average height
+!                  av_h=zs(k,jsj)-zs(k-1,jsj)
+!                  if(av_h<=0.d0) call parallel_abort('TRAN_IMP: av_h<=0')
+!                  !Check diffusion number; write warning message
+!                  difnum=dt_by_bigv*hdif_tmp/delj(jsj)*av_h*distj(jsj)
+!                  difnum_max_l=max(difnum_max_l,difnum)
+!                  rrhs(1,kin)=rrhs(1,kin)+difnum*(trel_tmp(m,k,iel)-trel_tmp(m,k,i))
+!                endif !k>=
+!              enddo !j
+!            endif !ihdif/=0
 
           enddo !k=kbe(i)+1,nvrt
 
@@ -2061,9 +2060,60 @@ subroutine do_transport_tvd_imp(it,ntr,difnum_max_l) !,nvrt1,npa1,dfh1)
       wtimer(9,2)=wtimer(9,2)+cwtmp2-cwtmp
 #endif
 
-!     Output warning for diffusion number
-      if(difnum_max_l>0.5d0) write(12,*)'TRAN_IMP, diffusion # exceeds 0.5:',it,difnum_max_l
-!'
+!     Horizontal diffusion as filter
+!     Save the final array from previous step as trel_tmp
+      if(ihdif/=0) then
+        do kk=1,niter_hdif
+!$OMP     parallel default(shared) !private(i,j,k,m,sum1,jsj,iel,nd1,nd2,hdif_tmp)
+!$OMP     workshare
+          trel_tmp(1:ntr,:,1:nea)=tr_el(1:ntr,:,1:nea)
+!$OMP     end workshare
+
+!$OMP     do
+          do i=1,ne
+            if(idry_e(i)==1) cycle
+
+            !Wet elements with wet nodes
+            do m=1,ntr !cycle through tracers
+              do k=kbe(i),nvrt
+                sum1=0.d0
+                do j=1,i34(i) !sides
+                  jsj=elside(j,i) !residents
+                  iel=ic3(j,i)
+                  if(iel==0.or.idry_e(max(1,iel))==1) cycle
+
+                  nd1=isidenode(1,jsj)
+                  nd2=isidenode(2,jsj)
+                  hdif_tmp=(hdif(nd1)+hdif(nd2))/2.d0
+                  if(k>=kbs(jsj)+1) then
+                    sum1=sum1+hdif_tmp*(trel_tmp(m,k,iel)-trel_tmp(m,k,i))
+                  endif !k>=
+                enddo !j
+
+                tr_el(m,k,i)=trel_tmp(m,k,i)+sum1
+              enddo !k
+
+              !Extend
+!              do k=1,kbe(i)
+!                tr_el(m,k,i)=tr_el(m,kbe(i)+1,i)
+!              enddo !k
+            enddo !m: tracers
+          enddo !i=1,ne
+!$OMP     end do
+!$OMP     end parallel
+
+!         Update ghosts
+#ifdef INCLUDE_TIMING
+          cwtmp=mpi_wtime()
+!          timer_ns(1)=timer_ns(1)+cwtmp-cwtmp2
+#endif
+          call exchange_e3d_2t_tr(tr_el)
+#ifdef INCLUDE_TIMING
+          cwtmp2=mpi_wtime()
+          wtimer(9,2)=wtimer(9,2)+cwtmp2-cwtmp
+#endif
+        enddo !kk=1,niter_hdif
+      endif !ihdif/=0
 
 #ifdef DEBUG
 !debug conservation