ww3_gse: testing some limiter stuff

erdc · Sep 22, 2023 · e2a9fa5 · e2a9fa5
1 parent ff00104
commit e2a9fa5
Show file tree

Hide file tree

Showing 7 changed files with 417 additions and 83 deletions.
diff --git a/model/src/PDLIB/yownodepool.F90 b/model/src/PDLIB/yownodepool.F90
@@ -77,7 +77,7 @@ module yowNodepool
 
   !> coordinates of the local +  ghost nodes. range [1:npa]
   real(rkind), public, target, allocatable :: x(:), y(:), z(:)
-  real(rkind), public, target, allocatable :: PDLIB_SI(:), PDLIB_TRIA(:), PDLIB_TRIA03(:), PDLIB_IEN(:,:)
+  real(rkind), public, target, allocatable :: PDLIB_SI(:), PDLIB_TRIA(:), PDLIB_TRIA03(:), PDLIB_IEN(:,:), PDLIB_IEND(:,:,:)
   integer, public, target, allocatable :: PDLIB_CCON(:), PDLIB_IA(:), PDLIB_JA(:)
   integer, public, target, allocatable :: PDLIB_IA_P(:), PDLIB_JA_P(:), PDLIB_JA_IE(:,:,:)
   integer, public, target, allocatable :: PDLIB_POS_CELL(:), PDLIB_IE_CELL(:)

diff --git a/model/src/PDLIB/yowpdlibmain.F90 b/model/src/PDLIB/yowpdlibmain.F90
@@ -1306,13 +1306,13 @@ subroutine ComputeTRIA_IEN_SI_CCON
     use yowerr,       only: parallel_abort
     use yowDatapool,    only: myrank
     use yowNodepool,    only: np_global, np, iplg, t_Node, ghostlg, ng, npa
-    use yowNodepool,    only: x, y, z, PDLIB_SI, PDLIB_IEN, PDLIB_TRIA, PDLIB_CCON, PDLIB_TRIA03
+    use yowNodepool,    only: x, y, z, PDLIB_SI, PDLIB_IEN, PDLIB_IEND, PDLIB_TRIA, PDLIB_CCON, PDLIB_TRIA03
 
     integer I1, I2, I3, stat, IE, NI(3)
     real  :: DXP1, DXP2, DXP3, DYP1, DYP2, DYP3, DBLTMP, TRIA03
     logical :: CROSSES_DATELINE
 
-    allocate(PDLIB_SI(npa), PDLIB_CCON(npa), PDLIB_IEN(6,ne), PDLIB_TRIA(ne), PDLIB_TRIA03(ne), stat=stat)
+    allocate(PDLIB_SI(npa), PDLIB_CCON(npa), PDLIB_IEN(6,ne), PDLIB_IEND(3,3,ne), PDLIB_TRIA(ne), PDLIB_TRIA03(ne), stat=stat)
     if(stat/=0) call parallel_abort('SI allocation failure')
 
     PDLIB_SI(:)   = 0.0d0 ! Median Dual Patch Area of each Node
@@ -1350,6 +1350,16 @@ subroutine ComputeTRIA_IEN_SI_CCON
         WRITE(*,*) 'AREA SMALLER ZERO IN PDLIB', IE, NE, PDLIB_TRIA(IE)
         STOP
       ENDIF
+      PDLIB_IEND(1,1,IE) = DOT_PRODUCT(PDLIB_IEN(1:2,IE),PDLIB_IEN(1:2,IE)) ! Tomaich. eq. 3.19 n1 * n1 etc.
+      PDLIB_IEND(1,2,IE) = DOT_PRODUCT(PDLIB_IEN(1:2,IE),PDLIB_IEN(3:4,IE))
+      PDLIB_IEND(1,3,IE) = DOT_PRODUCT(PDLIB_IEN(1:2,IE),PDLIB_IEN(5:6,IE))
+      PDLIB_IEND(2,1,IE) = DOT_PRODUCT(PDLIB_IEN(3:4,IE),PDLIB_IEN(1:2,IE))
+      PDLIB_IEND(2,2,IE) = DOT_PRODUCT(PDLIB_IEN(3:4,IE),PDLIB_IEN(3:4,IE))
+      PDLIB_IEND(2,3,IE) = DOT_PRODUCT(PDLIB_IEN(3:4,IE),PDLIB_IEN(5:6,IE))
+      PDLIB_IEND(3,1,IE) = DOT_PRODUCT(PDLIB_IEN(5:6,IE),PDLIB_IEN(1:2,IE))
+      PDLIB_IEND(3,2,IE) = DOT_PRODUCT(PDLIB_IEN(5:6,IE),PDLIB_IEN(3:4,IE))
+      PDLIB_IEND(3,3,IE) = DOT_PRODUCT(PDLIB_IEN(5:6,IE),PDLIB_IEN(5:6,IE))
+
 
       PDLIB_CCON(I1) = PDLIB_CCON(I1) + 1
       PDLIB_CCON(I2) = PDLIB_CCON(I2) + 1

diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90
@@ -802,9 +802,7 @@ MODULE W3GDATMD
 #ifdef W3_PR1
     REAL                  :: DUMMY
 #endif
-#ifdef W3_PR2
     REAL                  :: DTME, CLATMN
-#endif
 #ifdef W3_PR3
     REAL                  :: WDCG, WDTH
 #endif
@@ -1033,8 +1031,11 @@ MODULE W3GDATMD
     LOGICAL :: B_JGS_LIMITER
     LOGICAL :: B_JGS_USE_JACOBI
     LOGICAL :: B_JGS_BLOCK_GAUSS_SEIDEL
+    LOGICAL :: B_JGS_LGSE
+    INTEGER :: B_JGS_GSE_METHOD
     INTEGER :: B_JGS_MAXITER
     INTEGER :: B_JGS_LIMITER_FUNC
+    REAL*8  :: B_JGS_GSE_TS
     REAL*8  :: B_JGS_PMIN
     REAL*8  :: B_JGS_DIFF_THR
     REAL*8  :: B_JGS_NORM_THR
@@ -1231,9 +1232,7 @@ MODULE W3GDATMD
   !/
   !/ Data aliasses for structure PROP(S)
   !/
-#ifdef W3_PR2
   REAL, POINTER           :: DTME, CLATMN
-#endif
 #ifdef W3_PR3
   REAL, POINTER           :: WDCG, WDTH
 #endif
@@ -1388,6 +1387,9 @@ MODULE W3GDATMD
   LOGICAL, POINTER :: B_JGS_LIMITER
   LOGICAL, POINTER :: B_JGS_USE_JACOBI
   LOGICAL, POINTER :: B_JGS_BLOCK_GAUSS_SEIDEL
+  LOGICAL, POINTER :: B_JGS_LGSE
+  INTEGER, POINTER :: B_JGS_GSE_METHOD
+  REAL(8), POINTER :: B_JGS_GSE_TS
   INTEGER, POINTER :: B_JGS_MAXITER
   INTEGER, POINTER :: B_JGS_LIMITER_FUNC
   REAL(8), POINTER :: B_JGS_PMIN
@@ -2517,10 +2519,9 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     !
     !     Structure PROPS
     !
-#ifdef W3_PR2
     DTME   => MPARS(IMOD)%PROPS%DTME
     CLATMN => MPARS(IMOD)%PROPS%CLATMN
-#endif
+!AR: taken out i need this for gse ...
 #ifdef W3_PR3
     WDCG   => MPARS(IMOD)%PROPS%WDCG
     WDTH   => MPARS(IMOD)%PROPS%WDTH
@@ -2806,6 +2807,9 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     B_JGS_LIMITER => MPARS(IMOD)%SCHMS%B_JGS_LIMITER
     B_JGS_USE_JACOBI => MPARS(IMOD)%SCHMS%B_JGS_USE_JACOBI
     B_JGS_BLOCK_GAUSS_SEIDEL => MPARS(IMOD)%SCHMS%B_JGS_BLOCK_GAUSS_SEIDEL
+    B_JGS_LGSE => MPARS(IMOD)%SCHMS%B_JGS_LGSE
+    B_JGS_GSE_METHOD => MPARS(IMOD)%SCHMS%B_JGS_GSE_METHOD
+    B_JGS_GSE_TS => MPARS(IMOD)%SCHMS%B_JGS_GSE_TS
     B_JGS_MAXITER => MPARS(IMOD)%SCHMS%B_JGS_MAXITER
     B_JGS_LIMITER_FUNC => MPARS(IMOD)%SCHMS%B_JGS_LIMITER_FUNC
     B_JGS_PMIN => MPARS(IMOD)%SCHMS%B_JGS_PMIN

diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90
@@ -897,36 +897,39 @@ MODULE W3GRIDMD
 #ifdef W3_PR3
   REAL                    :: WDTHCG, WDTHTH
 #endif
-  LOGICAL :: JGS_TERMINATE_MAXITER
-  LOGICAL :: JGS_TERMINATE_DIFFERENCE
-  LOGICAL :: JGS_TERMINATE_NORM
-  LOGICAL :: JGS_LIMITER
-  INTEGER :: JGS_LIMITER_FUNC
-  LOGICAL :: JGS_BLOCK_GAUSS_SEIDEL
-  LOGICAL :: JGS_USE_JACOBI
-  LOGICAL :: JGS_SOURCE_NONLINEAR
-  LOGICAL :: UGOBCAUTO
-  LOGICAL :: UGBCCFL
-  LOGICAL :: EXPFSN
-  LOGICAL :: EXPFSPSI
-  LOGICAL :: EXPFSFCT
-  LOGICAL :: IMPFSN
-  LOGICAL :: EXPTOTAL
-  LOGICAL :: IMPTOTAL
-  LOGICAL :: IMPREFRACTION
-  LOGICAL :: IMPFREQSHIFT
-  LOGICAL :: IMPSOURCE
-  LOGICAL :: SETUP_APPLY_WLV
-  INTEGER :: JGS_MAXITER
+  LOGICAL :: JGS_TERMINATE_MAXITER = .TRUE.
+  LOGICAL :: JGS_TERMINATE_DIFFERENCE = .TRUE.
+  LOGICAL :: JGS_TERMINATE_NORM = .TRUE.
+  LOGICAL :: JGS_LIMITER = .FALSE.
+  INTEGER :: JGS_LIMITER_FUNC = 1
+  LOGICAL :: JGS_BLOCK_GAUSS_SEIDEL = .TRUE.
+  LOGICAL :: JGS_USE_JACOBI = .TRUE.
+  LOGICAL :: JGS_SOURCE_NONLINEAR = .FALSE.
+  LOGICAL :: JGS_LGSE = .FALSE. 
+  INTEGER :: JGS_GSE_METHOD = 1
+  REAL(8) :: JGS_GSE_TS = 1.d0
+  LOGICAL :: UGOBCAUTO = .FALSE.
+  LOGICAL :: UGBCCFL   = .FALSE.
+  LOGICAL :: EXPFSN    = .TRUE.
+  LOGICAL :: EXPFSPSI  = .FALSE.
+  LOGICAL :: EXPFSFCT  = .FALSE.
+  LOGICAL :: IMPFSN    = .FALSE.
+  LOGICAL :: EXPTOTAL  = .FALSE.
+  LOGICAL :: IMPTOTAL  = .FALSE.
+  LOGICAL :: IMPREFRACTION = .FALSE.
+  LOGICAL :: IMPFREQSHIFT = .FALSE.
+  LOGICAL :: IMPSOURCE = .FALSE.
+  LOGICAL :: SETUP_APPLY_WLV = .FALSE.
+  INTEGER :: JGS_MAXITER=100
   INTEGER :: nbSel
   INTEGER :: UNSTSCHEMES(6)
   INTEGER :: UNSTSCHEME
-  INTEGER :: JGS_NLEVEL
-  REAL*8  :: JGS_PMIN
-  REAL*8  :: JGS_DIFF_THR
-  REAL*8  :: JGS_NORM_THR
-  REAL*8  :: SOLVERTHR_SETUP
-  REAL*8  :: CRIT_DEP_SETUP
+  INTEGER :: JGS_NLEVEL = 0
+  REAL*8  :: JGS_PMIN = 0.
+  REAL*8  :: JGS_DIFF_THR = 1.E-10
+  REAL*8  :: JGS_NORM_THR = 1.E-20
+  REAL*8  :: SOLVERTHR_SETUP = 1.E-20
+  REAL*8  :: CRIT_DEP_SETUP = 0.
   !
   CHARACTER               :: UGOBCFILE*60
   REAL                    :: UGOBCDEPTH
@@ -1081,6 +1084,9 @@ MODULE W3GRIDMD
        JGS_LIMITER,                               &
        JGS_LIMITER_FUNC,                          &
        JGS_USE_JACOBI,                            &
+       JGS_LGSE,                                  &
+       JGS_GSE_METHOD,                            & 
+       JGS_GSE_TS,                                &
        JGS_BLOCK_GAUSS_SEIDEL,                    &
        JGS_MAXITER,                               &
        JGS_PMIN,                                  &
@@ -2400,8 +2406,11 @@ SUBROUTINE W3GRID()
     JGS_TERMINATE_NORM = .FALSE.
     JGS_LIMITER = .FALSE.
     JGS_LIMITER_FUNC = 1
+    JGS_GSE_TS = 350000
     JGS_BLOCK_GAUSS_SEIDEL = .TRUE.
     JGS_USE_JACOBI = .TRUE.
+    JGS_LGSE = .FALSE. 
+    JGS_GSE_METHOD = 1
     JGS_MAXITER=100
     JGS_PMIN = 1
     JGS_DIFF_THR = 1.E-10
@@ -2424,6 +2433,9 @@ SUBROUTINE W3GRID()
     B_JGS_NORM_THR = JGS_NORM_THR
     B_JGS_NLEVEL = JGS_NLEVEL
     B_JGS_SOURCE_NONLINEAR = JGS_SOURCE_NONLINEAR
+    B_JGS_LGSE = JGS_LGSE
+    B_JGS_GSE_TS = JGS_GSE_TS
+    B_JGS_GSE_METHOD = JGS_GSE_METHOD 
 
     nbSel=0
 

diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90
@@ -745,8 +745,11 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT )
              B_JGS_DIFF_THR,                                  &
              B_JGS_NORM_THR,                                  &
              B_JGS_NLEVEL,                                    &
-             B_JGS_SOURCE_NONLINEAR
-        !Init COUNTCON and IOBDP to zero, it needs to be set somewhere or
+             B_JGS_SOURCE_NONLINEAR,                          & 
+             B_JGS_LGSE,                                      &
+             B_JGS_GSE_TS,                                    & 
+             B_JGS_GSE_METHOD
+        !Init COUNTCON to zero, it needs to be set somewhere or
         !removed
         COUNTCON=0
         IOBDP=0  
@@ -849,7 +852,11 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT )
              B_JGS_DIFF_THR,                                  &
              B_JGS_NORM_THR,                                  &
              B_JGS_NLEVEL,                                    &
-             B_JGS_SOURCE_NONLINEAR
+             B_JGS_SOURCE_NONLINEAR,                          & 
+             B_JGS_LGSE,                                      &
+             B_JGS_GSE_TS,                                    &
+             B_JGS_GSE_METHOD
+
         IF (.NOT. GUGINIT) THEN
           CALL W3DIMUG ( IGRD, NTRI, NX, COUNTOT, NNZ, NDSE, NDST )
         END IF

diff --git a/model/src/w3iorsmd.F90 b/model/src/w3iorsmd.F90
@@ -653,7 +653,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
 #ifdef W3_MPI
         ELSE
           !
-          IF (LPDLIB) THEN
+          IF (LPDLIB .and. (GTYPE.eq.UNGTYPE)) THEN
 #endif
 #ifdef W3_TIMINGS
             CALL PRINT_MY_TIME("Before UNST_PDLIB_WRITE_TO_FILE")
@@ -735,7 +735,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
         WRITE (NDST,9020) TYPE
 #endif
       ELSE
-        IF (LPDLIB) THEN
+        IF (LPDLIB .and. (GTYPE.eq.UNGTYPE)) THEN
 #ifdef W3_TIMINGS
           CALL PRINT_MY_TIME("Before UNST_PDLIB_READ_FROM_FILE")
 #endif
@@ -796,7 +796,8 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
           ELSE
 #endif
             VA = 0.
-            DO JSEA=1, NSEA
+!AR: 2do here is most likely a bug in develop there NSEA was used, which is global
+            DO JSEA=1, NSEAL
               CALL INIT_GET_ISEA(ISEA, JSEA)
               NREC   = ISEA + 2
               RPOS   = 1_8 + LRECL*(NREC-1_8)