Bug fix in misc_subs c/o of Fei: take care of odd cases where extrap

node is on interface
schism-dev · Apr 21, 2021 · 0cec024 · 0cec024
1 parent b96bc94
commit 0cec024
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 28 deletions.
diff --git a/src/Hydro/misc_subs.F90 b/src/Hydro/misc_subs.F90
@@ -1032,37 +1032,36 @@ subroutine levels1(iths,it)
 
 !         Interface (shoreline) sides
 !$OMP     workshare
-          icolor=0 !nodes on the interface sides
+!          icolor=0 !nodes on the interface sides (not needed)
           icolor2=0 !interface sides
 !$OMP     end workshare
 
 !$OMP     do
           do i=1,ns
             if(isdel(2,i)/=0) then; if(idry_e2(isdel(1,i))+idry_e2(isdel(2,i))==1) then
-!              icolor(isidenode(1:2,i))=1
               icolor2(i)=1
             endif; endif
           enddo !i
 !$OMP     end do
 
-!$OMP     do
-          loopinun: do i=1,np
-            do j=1,nne(i)
-              ie=indel(j,i)
-              id=iself(j,i)
-              do m=1,2 !2 neighboring sides
-                isd=elside(nxq(m+i34(ie)-3,id,i34(ie)),ie)
-                if(icolor2(isd)==1) then
-                  icolor(i)=1
-                  cycle loopinun
-                endif
-              enddo !m
-            enddo !j
-          end do loopinun !i
-!$OMP     end do
+!!$OMP     do
+!          loopinun: do i=1,np
+!            do j=1,nne(i)
+!              ie=indel(j,i)
+!              id=iself(j,i)
+!              do m=1,2 !2 neighboring sides
+!                isd=elside(nxq(m+i34(ie)-3,id,i34(ie)),ie)
+!                if(icolor2(isd)==1) then
+!                  icolor(i)=1
+!                  cycle loopinun
+!                endif
+!              enddo !m
+!            enddo !j
+!          end do loopinun !i
+!!$OMP     end do
 !$OMP end parallel
 
-          call exchange_p2di(icolor)
+!          call exchange_p2di(icolor)
           call exchange_s2di(icolor2)
 
 !         Aug. shoreline sides (must be internal sides)
@@ -1087,25 +1086,35 @@ subroutine levels1(iths,it)
             inew=0 !for initializing and counting su2 sv2
             do i=1,nsdf !aug.
               isd=isdf(i)
-              if(isdel(1,isd)<0.or.isdel(2,isd)<0) cycle
+!              if(isdel(1,isd)<0.or.isdel(2,isd)<0) cycle
               if(isdel(1,isd)==0.or.isdel(2,isd)==0) then
                 write(errmsg,*)'LEVELS1: bnd side (2):',isdel(:,isd),iplg(isidenode(1:2,isd))
                 call parallel_abort(errmsg)
               endif
-              if(idry_e2(isdel(1,isd))+idry_e2(isdel(2,isd))/=1) cycle
+!              if(idry_e2(isdel(1,isd))+idry_e2(isdel(2,isd))/=1) cycle
+
+              !Try to find a dry elem (to take care of some odd cases where
+              !nodeA is interface btw sub-domains)
+!              if(idry_e2(isdel(1,isd))==1) then
+!                ie=isdel(1,isd)
+!              else 
+!                ie=isdel(2,isd)
+!              endif
+              ie=0
+              do m=1,2
+                if(isdel(m,isd)>0) then; if(idry_e2(isdel(m,isd))==1) then
+                  ie=isdel(m,isd); exit
+                endif; endif
+              enddo !m
+              if(ie==0) cycle
 
-              if(idry_e2(isdel(1,isd))==1) then
-                ie=isdel(1,isd)
-              else 
-                ie=isdel(2,isd)
-              endif
               n1=isidenode(1,isd)
               n2=isidenode(2,isd)
               nodeA=elnode(1,ie)+elnode(2,ie)+elnode(3,ie)-n1-n2
 
               if(icolor(nodeA)==1) cycle !this node is done
 
-              icolor(nodeA)=1 !this node is done
+              icolor(nodeA)=1 !this node will be done
               if(nodeA>np) cycle
 !             nodeA is resident
 
@@ -1191,8 +1200,8 @@ subroutine levels1(iths,it)
               ltmp=ltmp.or.inew(i)/=0
               if(inew(i)/=0) then
 !                srwt_xchng(1)=.true.
-                su2(1:nvrt,i)=su2(1:nvrt,i)/inew(i)
-                sv2(1:nvrt,i)=sv2(1:nvrt,i)/inew(i)
+                su2(1:nvrt,i)=su2(1:nvrt,i)/dble(inew(i))
+                sv2(1:nvrt,i)=sv2(1:nvrt,i)/dble(inew(i))
               endif
             enddo !i
 !$OMP end parallel do 

diff --git a/src/Readme.beta_notes b/src/Readme.beta_notes
@@ -480,6 +480,7 @@ git versions:
 (118) 28ff9d1 (Dec 16, 2020): added optional self-attraction loading tides. The option shares some constants
                               with tidal potential: freq names and cut-off depth;
 (119) b1bcaa0 (April 20, 2021): removed most of 'goto'. Remaining ones: harm.F90, lap.F90, WWM
+(120) (April 21, 2021): bug fix in misc_subs c/o of Fei (nodeA in final extrap stage may be interface node)
 ---------------------------------------------------------------
 Auto-test history:
 R100 (Nov 2011); R168 (minus WWM); R224 (basic); R240 (full); R306 (basic); R370 (basic); R408 (basic); R430 (basic); R601 (basic); R730 (basic); R747 (all hydro); R1120 (hydro+ SF BayDelta); R1305: (all hydro; ICM); R1532 (branch/selfe_opt1): all; R1641: hydro; R2018: hydro; R2232 (all);