From 0aa1b5ad397b383f99a2d4edacdafb5d7194badb Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Thu, 26 Oct 2023 21:00:25 +0200 Subject: [PATCH 01/10] ww3_gse2: merge gse stuff ... --- model/src/w3gdatmd.F90 | 15 +- model/src/w3gridmd.F90 | 66 +++--- model/src/w3profsmd_pdlib.F90 | 378 ++++++++++++++++++++++++++++++---- 3 files changed, 387 insertions(+), 72 deletions(-) diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90 index 7bc5e2f30..449c9b81d 100644 --- a/model/src/w3gdatmd.F90 +++ b/model/src/w3gdatmd.F90 @@ -813,9 +813,7 @@ MODULE W3GDATMD #ifdef W3_PR1 REAL :: DUMMY #endif -#ifdef W3_PR2 REAL :: DTME, CLATMN -#endif #ifdef W3_PR3 REAL :: WDCG, WDTH #endif @@ -1046,8 +1044,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 @@ -1244,9 +1245,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 @@ -1403,6 +1402,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 @@ -2532,10 +2534,8 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST ) ! ! Structure PROPS ! -#ifdef W3_PR2 DTME => MPARS(IMOD)%PROPS%DTME CLATMN => MPARS(IMOD)%PROPS%CLATMN -#endif #ifdef W3_PR3 WDCG => MPARS(IMOD)%PROPS%WDCG WDTH => MPARS(IMOD)%PROPS%WDTH @@ -2829,6 +2829,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 index fa8128afb..f2fef2bdd 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -903,36 +903,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 @@ -1089,6 +1092,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, & @@ -2432,8 +2438,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 @@ -2456,6 +2465,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/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 6759fb53e..fc8723a9d 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -973,7 +973,7 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) REAL :: FL111, FL112, FL211, FL212, FL311, FL312 REAL :: DTSI(npa), U(npa) REAL :: DTMAX_GL, DTMAX, DTMAXEXP, REST - REAL :: LAMBDA(2), KTMP(3) + REAL :: LAMBDA(2), KTMP(3), DFAK(NE) REAL :: KELEM(3,NE), FLALL(3,NE) REAL :: KKSUM(npa), ST(npa) REAL :: NM(NE) @@ -1105,7 +1105,7 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) DO IE = 1, NE NI = INE(:,IE) UTILDE = NM(IE) * (DOT_PRODUCT(FLALL(:,IE),U(NI))) - ST(NI) = ST(NI) + KELEM(:,IE) * (U(NI) - UTILDE) ! the 2nd term are the theta values of each node ... + ST(NI) = ST(NI) + KELEM(:,IE) * (U(NI) - UTILDE) END DO ! IE #ifdef W3_DEBUGSOLVER IF (testWrite) THEN @@ -5495,7 +5495,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL USE CONSTANTS, only : TPI, TPIINV, GRAV USE W3GDATMD, only: MAPSTA USE W3GDATMD, only: FSREFRACTION, FSFREQSHIFT, FSSOURCE, NX, DSIP - USE W3GDATMD, only: B_JGS_NORM_THR, B_JGS_TERMINATE_NORM, B_JGS_PMIN + USE W3GDATMD, only: B_JGS_NORM_THR, B_JGS_TERMINATE_NORM, B_JGS_PMIN, B_JGS_LIMITER_FUNC USE W3GDATMD, only: B_JGS_TERMINATE_DIFFERENCE, B_JGS_MAXITER, B_JGS_LIMITER USE W3GDATMD, only: B_JGS_TERMINATE_MAXITER, B_JGS_BLOCK_GAUSS_SEIDEL, B_JGS_DIFF_THR USE W3GDATMD, only: MAPWN @@ -6239,21 +6239,6 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL IF (FLSOU) THEN IF (B_JGS_LIMITER) THEN - - DO ISP=1,NSPEC - IK = 1 + (ISP-1)/NTH - SPEC(ISP) = VAOLD(ISP,JSEA) - ENDDO -#ifdef W3_ST4 - CALL W3SPR4 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEAN1, WNMEAN, & - AMAX, U10(ISEA), U10D(ISEA), & -#ifdef W3_FLX5 - TAUA, TAUADIR, DAIR, & -#endif - USTAR, USTDIR, & - TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS, DLWMEAN) -#endif - DAM = 0. DO IK=1, NK DAM(1+(IK-1)*NTH) = 0.0081*0.1 / ( 2 * SIG(IK) * WN(IK,ISEA)**3 * CG(IK,ISEA)) * CG1(IK) / CLATS(ISEA) @@ -6266,28 +6251,52 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL END DO END DO - DAM2 = 0. - DO IK=1, NK - JAC2 = 1./TPI/SIG(IK) - FRLOCAL = SIG(IK)*TPIINV - DAM2(1+(IK-1)*NTH) = 1E-06 * GRAV/FRLOCAL**4 * USTAR * MAX(FMEANWS,FMEAN) * DTG * JAC2 * CG1(IK) / CLATS(ISEA) - END DO - DO IK=1, NK - IS0 = (IK-1)*NTH - DO ITH=2, NTH - DAM2(ITH+IS0) = DAM2(1+IS0) + IF (B_JGS_LIMITER_FUNC == 2) THEN + DO ISP=1,NSPEC + IK = 1 + (ISP-1)/NTH + SPEC(ISP) = VAOLD(ISP,JSEA) + ENDDO +#ifdef W3_ST4 + CALL W3SPR4 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEAN1, WNMEAN, & + AMAX, U10(ISEA), U10D(ISEA), & +#ifdef W3_FLX5 + TAUA, TAUADIR, DAIR, & +#endif + USTAR, USTDIR, & + TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS, DLWMEAN) +#endif + DAM2 = 0. + DO IK=1, NK + JAC2 = 1./TPI/SIG(IK) + FRLOCAL = SIG(IK)*TPIINV + DAM2(1+(IK-1)*NTH) = 5E-07 * GRAV/FRLOCAL**4 * UST(ISEA) * FMEAN * DTG * JAC2 * CG1(IK) / CLATS(ISEA) END DO - END DO - - DO IK = 1, NK - DO ITH = 1, NTH - ISP = ITH + (IK-1)*NTH - newdac = VA(ISP,IP) - VAOLD(ISP,JSEA) - maxdac = max(DAM(ISP),DAM2(ISP)) - NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC) - VA(ISP,IP) = max(0., VAOLD(ISP,IP) + NEWDAC) + DO IK=1, NK + IS0 = (IK-1)*NTH + DO ITH=2, NTH + DAM2(ITH+IS0) = DAM2(1+IS0) + END DO + END DO + DO IK = 1, NK + DO ITH = 1, NTH + ISP = ITH + (IK-1)*NTH + newdac = VA(ISP,IP) - VAOLD(ISP,JSEA) + maxdac = max(DAM(ISP),DAM2(ISP)) + NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC) + VA(ISP,IP) = max(0., VAOLD(ISP,IP) + NEWDAC) + ENDDO ENDDO - ENDDO + ELSE + DO IK = 1, NK + DO ITH = 1, NTH + ISP = ITH + (IK-1)*NTH + newdac = VA(ISP,IP) - VAOLD(ISP,JSEA) + maxdac = DAM(ISP) + NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC) + VA(ISP,IP) = max(0., VAOLD(ISP,IP) + NEWDAC) + ENDDO + ENDDO + ENDIF ENDIF ! B_JGS_LIMITER ENDIF ! FLSOU END DO ! JSEA @@ -7630,4 +7639,295 @@ SUBROUTINE JACOBI_FINALIZE !/ END SUBROUTINE JACOBI_FINALIZE !/ ------------------------------------------------------------------- / -END MODULE PDLIB_W3PROFSMD + SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) + !/ + !/ +-----------------------------------+ + !/ | WAVEWATCH III NOAA/NCEP | + !/ | | + !/ | Aron Roland (Roland & Partner, | + !/ | BGS IT&E GmbH) | + !/ | | + !/ | FORTRAN 90 | + !/ | Last update : 01-June-2023 | + !/ +-----------------------------------+ + !/ + !/ 01-June-2018 : Origination. ( version 6.04 ) + !/ + ! 1. Purpose : Solves the anisotropic diffusion equation + ! 2. Method :i Galerkin + ! 3. Parameters : a lot + ! + ! Parameter list + ! ---------------------------------------------------------------- + ! ---------------------------------------------------------------- + ! + ! 4. Subroutines used : + ! + ! Name Type Module Description + ! ---------------------------------------------------------------- + ! STRACE Subr. W3SERVMD Subroutine tracing. + ! ---------------------------------------------------------------- + ! + ! 5. Called by : + ! + ! Name Type Module Description + ! ---------------------------------------------------------------- + ! ---------------------------------------------------------------- + ! + ! 6. Error messages : + ! 7. Remarks + ! 8. Structure : + ! 9. Switches : + ! + ! !/S Enable subroutine tracing. + ! + ! 10. Source code : + ! + !/ ------------------------------------------------------------------- / +#ifdef W3_S + USE W3SERVMD, only: STRACE +#endif + ! + USE CONSTANTS, only : LPDLIB, TPI, TPIINV, GRAV, DERA, RADIUS + USE W3GDATMD, only: MAPSF, NSEAL, DMIN, IOBDP, MAPSTA, IOBP, MAPFS, NX, CLATS, CLATMN, SIG + USE W3GDATMD, only: ESIN, ECOS, XFR, DTH, NSPEC, B_JGS_GSE_TS, XGRD, YGRD + USE W3ADATMD, only: DW, MPI_COMM_WCMP, U10, WN, CG + USE W3WDATMD, ONLY: VA + USE W3PARALL, only: INIT_GET_ISEA + USE W3GDATMD, only: IOBP_LOC, IOBPD_LOC, IOBPA_LOC, IOBDP_LOC + USE YOWNODEPOOL, only: iplg, np, pdlib_tria, pdlib_ien, pdlib_si + USE yowfunction, only: pdlib_abort + use YOWNODEPOOL, only: npa + use yowElementpool, only : INE, NE + use yowDatapool, only: rtype + use yowExchangeModule, only : PDLIB_exchange1DREAL + USE W3GDATMD, only: B_JGS_USE_JACOBI, FLAGLL + USE W3GDATMD, only: NSPEC, NTH, NK + USE W3GDATMD, only: FSTOTALIMP + USE W3ODATMD, only: IAPROC + USE MPI, only : MPI_MIN + !/ + REAL, INTENT(IN) :: DTG + ! + !/ ------------------------------------------------------------------- / + !/ + INTEGER ITH, IK, IE, IS, IERR, MYRANK + INTEGER NewISP, JTH, istati, JSEA, ISEA + REAL*8 TFAC, DNN, DSS, DSSD, DNND, CGD, DFRR, DTME + REAL PHI_V(NPA) + REAL*8 eDet, DEDX(3), DEDY(3), DIFFV(2) + INTEGER NI(3), ITR, IDX, IP, ISP, IT + REAL*8 XSEL(3), DVDXIE, DVDYIE, FACX + REAL*8 GRAD(2), V(2), eScal, DT_DIFF, DIFFVEC(3,NPA) + INTEGER NB_ITER, iIter, ip_global + REAL*8 DeltaTmax, eDeltaT, CLATSMN, DFAC, RFAC, eDiffNorm + REAL*8 eNorm, DTquot, diffc, dcell, XWIND, DIFFTOT + REAL*8 DVDX(NPA), DVDY(NPA), DV2DXX(NPA), DV2DXY(NPA), KH, SWFAC + REAL eRealA(1), eRealB(1) + + + DVDX = 0.d0 + DVDY = 0.d0 + DV2DXX = 0.d0 + DV2DXY = 0.d0 + + DTME = B_JGS_GSE_TS + FACX = 1./(DERA * RADIUS) + !WRITE(*,*) 'START BLOCK SOLVER' + CALL MPI_COMM_RANK(MPI_COMM_WCMP, myrank, ierr) + + !WRITE (3000+myrank,*) 'Entering Diffusion' + !CALL MPI_BARRIER (MPI_COMM_WCMP,IERR) + + IF ( FLAGLL ) THEN + RFAC = DERA * RADIUS + DFAC = 1. / RFAC**2 + ELSE + RFAC = 1. + DFAC = 1. + END IF + CLATMN = COS ( 70 * DERA ) + + DO ITH = 1, NTH + DO IK = 1, NK + ISP = ITH + (IK-1)*NTH + DO JSEA = 1, NPA + CALL INIT_GET_ISEA(ISEA, JSEA) + TFAC = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 ) + CGD = 0.5 * GRAV / SIG(IK) * IOBDP_LOC(JSEA) + DSS = ( CGD * (XFR-1.))**2 * DTME / 12.! * TFAC + DNN = ( CGD * DTH )**2 * DTME / 12.! * TFAC + DCELL = CGD / 10.0 ! -> CELLP needs probably redifinition ... + KH = WN(IK,ISEA)*DW(JSEA) + IF (DW(JSEA) .gt. 1000.) THEN + SWFAC = 1. + ELSE + SWFAC = 0.!(1.-0.5*(1+(2*KH/SINH(MIN(20.d0,2*KH)))))**8 + ENDIF + XWIND = 3.3 * U10(ISEA)*WN(IK,ISEA)/SIG(IK) - 2.3 + XWIND = MAX ( 0. , MIN ( 1. , XWIND ) ) +#ifdef W3_XW0 + XWIND = 0. +#endif +#ifdef W3_XW1 + XWIND = 1. +#endif + TFAC = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 ) + DSS = SWFAC * (XWIND * DCELL + (1.-XWIND) * DSS * TFAC) +#ifdef W3_DSS0 + DSS = 0. +#endif + DNN = SWFAC * (XWIND * DCELL + (1.-XWIND) * DNN * TFAC) + + DIFFVEC(1,JSEA) = (DSS*ECOS(ITH)**2+DNN*ESIN(ITH)**2) + DIFFVEC(2,JSEA) = (DSS*ESIN(ITH)**2+DNN*ECOS(ITH)**2) / CLATS(ISEA)**2 + DIFFVEC(3,JSEA) = ((DSS-DNN) * ESIN(ITH)*ECOS(ITH)) / CLATS(ISEA) + !write(3000+myrank,'(I10,20F20.10)') IK, SIG(IK), CGD, CLATS(ISEA), DSS, DNN, DIFFVEC(1,JSEA), DIFFVEC(2,JSEA) + END DO + + CALL DIFFERENTIATE_XYDIR(DBLE(VA(ISP,:)),DVDX,DVDY) ! AR: this two lines can be optimized, however seems fast enough by now. + CALL DIFFERENTIATE_XYDIR(DVDX,DV2DXX,DV2DXY) + + DeltaTmax = 1./TINY(1.) + DO IE=1,NE + NI = INE(:,IE) + eDet = 2. * PDLIB_TRIA(IE) + DO IDX=1,3 + V(1) = PDLIB_IEN(2*IDX-1,IE) + V(2) = PDLIB_IEN(2*IDX ,IE) + eNorm = DOT_PRODUCT(V,V) + IP = INE(IDX,IE) + eDiffNorm = HYPOT(DIFFVEC(1,IP), DIFFVEC(2,IP)) + if (eDiffNorm .gt. tiny(1.)) then + eDeltaT = min(DTG, PDLIB_SI(IP) / (eDiffNorm * dfac * eNorm /(4. * eDet))) ! modified for diffusion vector + else + eDeltaT = DTG + endif + !write(4000+myrank,*) IE, iplg(IP), eDeltaT, eDiffNorm, eNorm, 4 * eDet + IF (eDeltaT .lt. DeltaTmax) THEN + DeltaTmax = eDeltaT + END IF + END DO + END DO + + eRealA(1)=DeltaTmax + CALL MPI_ALLREDUCE(eRealA, eRealB, 1, rtype, MPI_MIN, MPI_COMM_WCMP, ierr) + DeltaTmax=eRealB(1) + + DTquot = DTG / DeltaTmax + NB_ITER = NINT(DTquot) + 1 + DT_DIFF = DTG/NB_ITER + PHI_V = 0. + + !WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax + !CALL FLUSH(5000+myrank) + + DO IT = 1, NB_ITER + DO IE = 1, NE + NI = INE(:,IE) + eDet = 2. * PDLIB_TRIA(IE) + DEDX(1) = PDLIB_IEN(1,IE) + DEDX(2) = PDLIB_IEN(3,IE) + DEDX(3) = PDLIB_IEN(5,IE) + DEDY(1) = PDLIB_IEN(2,IE) + DEDY(2) = PDLIB_IEN(4,IE) + DEDY(3) = PDLIB_IEN(6,IE) + XSEL = VA(ISP,NI) + DVDXIE = DOT_PRODUCT(XSEL,DEDX) + DVDYIE = DOT_PRODUCT(XSEL,DEDY) + GRAD(1) = DVDXIE / eDet * 1./3. * SUM(DIFFVEC(1,NI)) + GRAD(2) = DVDYIE / eDet * 1./3. * SUM(DIFFVEC(2,NI)) + DO IDX=1,3 + V(1) = 0.5 * PDLIB_IEN(2*IDX-1,IE) + V(2) = 0.5 * PDLIB_IEN(2*IDX ,IE) + eScal = DOT_PRODUCT(V, GRAD(1:2)) + IP = INE(IDX,IE) + PHI_V(IP) = PHI_V(IP) + eScal + END DO + END DO + CALL PDLIB_exchange1DREAL(PHI_V) + DO JSEA =1, NSEAL + VA(ISP,JSEA) = MAX(0.,VA(ISP,JSEA) - (DT_DIFF * PHI_V(JSEA) / PDLIB_SI(JSEA) + 2 * DV2DXY(JSEA) * DIFFVEC(3,JSEA)) * DFAC * IOBDP_LOC(JSEA)) + END DO + END DO + END DO + END DO + + !WRITE(3000+myrank,*) 'FINISHED DIFFUSION' + !CALL FLUSH(3000+myrank) + !CALL MPI_BARRIER (MPI_COMM_WCMP,IERR) + + END SUBROUTINE BLOCK_SOLVER_DIFFUSION + + SUBROUTINE DIFFERENTIATE_XYDIR(VAR, DVDX, DVDY) + USE W3GDATMD, only: ECOS, ESIN, DMIN, NTH, SIG, NK + USE W3ADATMD, only: CG, CX, CY, DW + USE yowExchangeModule, only : PDLIB_exchange1DREAL + USE yowNodepool, only : PDLIB_IEN, PDLIB_TRIA, NP, NPA + USE yowElementpool, only : NE, INE + IMPLICIT NONE + REAL*8, INTENT(IN) :: VAR(NPA) + REAL*8, INTENT(INOUT) :: DVDX(NPA), DVDY(NPA) + REAL*4 :: DVDX4(NPA), DVDY4(NPA) + REAL*8 :: DEDY(3),DEDX(3) + REAL*8 :: DVDXIE, DVDYIE + REAL*8 :: WEI(NPA) + INTEGER :: NI(3) + INTEGER :: IE, JSEA + + WEI(:) = 0.d0 + DVDX(:) = 0.d0 + DVDY(:) = 0.d0 + +#ifdef DEBUG + WRITE(STAT%FHNDL,*) 'DIFFERENTIATE_XYDIR' + WRITE(STAT%FHNDL,*) 'sum(VAR ) = ', sum(VAR) + WRITE(STAT%FHNDL,*) 'sum(IEN ) = ', sum(IEN) + WRITE(STAT%FHNDL,*) 'sum(TRIA) = ', sum(TRIA) +#endif + + DO IE = 1, NE + NI = INE(:,IE) + WEI(NI) = WEI(NI) + 2.*PDLIB_TRIA(IE) + DEDX(1) = PDLIB_IEN(1,IE) + DEDX(2) = PDLIB_IEN(3,IE) + DEDX(3) = PDLIB_IEN(5,IE) + DEDY(1) = PDLIB_IEN(2,IE) + DEDY(2) = PDLIB_IEN(4,IE) + DEDY(3) = PDLIB_IEN(6,IE) + DVDXIE = DOT_PRODUCT( VAR(NI),DEDX) + DVDYIE = DOT_PRODUCT( VAR(NI),DEDY) + DVDX(NI) = DVDX(NI) + DVDXIE + DVDY(NI) = DVDY(NI) + DVDYIE + END DO + + DVDX = DVDX/WEI + DVDY = DVDY/WEI + + DO JSEA = 1, NP + IF (DW(JSEA) .LT. DMIN) THEN + DVDX(JSEA) = 0. + DVDY(JSEA) = 0. + END IF + END DO + + DVDX4 = DVDX + DVDY4 = DVDY + CALL PDLIB_exchange1DREAL(DVDX4) ! AR: todo, checck pipes. + CALL PDLIB_exchange1DREAL(DVDY4) + DVDX = DVDX4 + DVDY = DVDY4 + +#ifdef DEBUG + WRITE(STAT%FHNDL,*) 'sum(DVDX) = ', sum(DVDX) + WRITE(STAT%FHNDL,*) 'sum(DVDY) = ', sum(DVDY) +#endif + + IF (.FALSE.) THEN + OPEN(2305, FILE = 'erggrad.bin' , FORM = 'UNFORMATTED') + WRITE(2305) 1. + WRITE(2305) (DVDX(JSEA), DVDY(JSEA), SQRT(DVDY(JSEA)**2+DVDY(JSEA)**2), JSEA = 1, NP) + ENDIF + END SUBROUTINE + + END MODULE PDLIB_W3PROFSMD From beee8618f98b469ca0772b454fbe90d705fbbac0 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Sun, 5 Nov 2023 17:31:58 +0100 Subject: [PATCH 02/10] ww3_gse2: merge some missing parts ... --- model/src/PDLIB/yownodepool.F90 | 2 +- model/src/PDLIB/yowpdlibmain.F90 | 14 ++++++++++++-- model/src/w3profsmd_pdlib.F90 | 4 ++-- model/src/w3wavemd.F90 | 11 ++++++++++- 4 files changed, 25 insertions(+), 6 deletions(-) diff --git a/model/src/PDLIB/yownodepool.F90 b/model/src/PDLIB/yownodepool.F90 index 8e98c3fc5..e005ebea2 100644 --- 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 index cc72b97fd..76dfb616e 100644 --- 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/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index fc8723a9d..29b1770ea 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -915,7 +915,7 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) USE W3SERVMD, only: STRACE #endif ! - USE W3GDATMD, only: NK, NTH, NX, IEN, CLATS, MAPSF + USE W3GDATMD, only: NK, NTH, NX, IEN, CLATS, MAPSF, DTH USE W3GDATMD, only: IOBPD_LOC, IOBP_LOC, IOBDP_LOC, IOBPA_LOC, FSBCCFL USE W3WDATMD, only: TIME USE W3ADATMD, only: CG, ITER, DW , CFLXYMAX, NSEALM @@ -926,7 +926,7 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) #ifdef W3_REF1 USE W3GDATMD, only: REFPARS #endif - USE YOWNODEPOOL, only: PDLIB_SI, PDLIB_IEN, PDLIB_TRIA, ipgl, iplg, npa, np + USE YOWNODEPOOL, only: PDLIB_SI, PDLIB_IEN, PDLIB_IEND, PDLIB_TRIA, ipgl, iplg, npa, np use yowElementpool, only: ne, INE use yowDatapool, only: rtype use yowExchangeModule, only : PDLIB_exchange1DREAL diff --git a/model/src/w3wavemd.F90 b/model/src/w3wavemd.F90 index 44c80964d..0545eb396 100644 --- a/model/src/w3wavemd.F90 +++ b/model/src/w3wavemd.F90 @@ -448,7 +448,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & USE W3IOSFMD #ifdef W3_PDLIB USE PDLIB_W3PROFSMD, only : APPLY_BOUNDARY_CONDITION_VA - USE PDLIB_W3PROFSMD, only : PDLIB_W3XYPUG, PDLIB_W3XYPUG_BLOCK_IMPLICIT, PDLIB_W3XYPUG_BLOCK_EXPLICIT + USE PDLIB_W3PROFSMD, only : PDLIB_W3XYPUG, PDLIB_W3XYPUG_BLOCK_IMPLICIT, PDLIB_W3XYPUG_BLOCK_EXPLICIT, block_solver_diffusion USE PDLIB_W3PROFSMD, only : ALL_VA_INTEGRAL_PRINT, ALL_VAOLD_INTEGRAL_PRINT, ALL_FIELD_INTEGRAL_PRINT USE W3PARALL, only : PDLIB_NSEAL, PDLIB_NSEALM USE yowNodepool, only: npa, iplg, np @@ -1819,6 +1819,9 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & DO ISPEC=1,NSPEC CALL PDLIB_W3XYPUG ( ISPEC, FACX, FACX, DTG, VGX, VGY, UGDTUPDATE ) END DO + IF (B_JGS_LGSE) THEN + CALL BLOCK_SOLVER_DIFFUSION(DTG) + ENDIF END IF END IF #endif @@ -1831,12 +1834,18 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & #endif #ifdef W3_PDLIB CALL PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACX, DTG, VGX, VGY, UGDTUPDATE ) + IF (B_JGS_LGSE) THEN + CALL BLOCK_SOLVER_DIFFUSION(DTG) + ENDIF #endif #ifdef W3_PDLIB ELSE IF(FSTOTALEXP .and. (IT .ne. 0)) THEN #endif #ifdef W3_PDLIB CALL PDLIB_W3XYPUG_BLOCK_EXPLICIT(IMOD, FACX, FACX, DTG, VGX, VGY, UGDTUPDATE ) + IF (B_JGS_LGSE) THEN + CALL BLOCK_SOLVER_DIFFUSION(DTG) + ENDIF #endif #ifdef W3_PDLIB ENDIF From 15696ac79ae4a0474232eec25decbec338152628 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Tue, 16 Jan 2024 21:55:36 +0100 Subject: [PATCH 03/10] ww3_gse2: describe new switches in the namelist file --- model/nml/namelists.nml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/model/nml/namelists.nml b/model/nml/namelists.nml index 390fdb874..f4172ddec 100644 --- a/model/nml/namelists.nml +++ b/model/nml/namelists.nml @@ -318,6 +318,8 @@ $ JGS_DIFF_THR : implicit solver threshold for JGS_TE $ JGS_NORM_THR : terminate based on the norm of the solution $ JGS_LIMITER : use total (quasi-steady: limits whole equation) instead of local limiter (un-steady: limits only source terms) $ JGS_LIMITER_FUNC : 1 - old limiter; 2 - alternatnive limiter +$ JGS_LGSE : T/F - turn on/off GSE correction on unstructured grids +$ JGS_GSE_TS : Time constant of the diffusion tensor, see Manual, same as for PR2 $ SETUP_APPLY_WLV : Compute wave setup (experimental) $ SOLVERTHR_SETUP : Solver threshold for setup computations $ CRIT_DEP_SETUP : Critical depths for setup computations From ac5b1045235db4c0ea44f82fe43dc84c803e66e3 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Tue, 16 Jan 2024 21:56:36 +0100 Subject: [PATCH 04/10] ww3_gse2: update inp examples --- model/inp/ww3_grid.inp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/model/inp/ww3_grid.inp b/model/inp/ww3_grid.inp index 655a10493..a3223b9c0 100644 --- a/model/inp/ww3_grid.inp +++ b/model/inp/ww3_grid.inp @@ -352,6 +352,8 @@ $ JGS_LIMITER : TRUE: Use total (quasi-steady: limits $ FALSE: default $ JGS_LIMITER_FUNC : 1 - old limiter (default) $ 2 - alternatnive limiter +$ JGS_LGSE : T/F - turn on/off GSE correction on unstructured grids +$ JGS_GSE_TS : Time constant of the diffusion tensor, see Manual, same as for PR2 $ SETUP_APPLY_WLV : Compute wave setup (TRUE/FALSE, default TRUE) $ SOLVERTHR_SETUP : Solver threshold for setup computations (default 1E-6) $ CRIT_DEP_SETUP : Critical depth for setup computations (default 0.1) From 9e07bbe112c5643e560a71b1c712fe7eb1b08707 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Tue, 16 Jan 2024 22:03:33 +0100 Subject: [PATCH 05/10] ww3_gse2: revert init part in w3grid --- model/src/w3gridmd.F90 | 60 +++++++++++++++++++++--------------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90 index f2fef2bdd..e1151d29d 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -903,39 +903,39 @@ MODULE W3GRIDMD #ifdef W3_PR3 REAL :: WDTHCG, WDTHTH #endif - 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 + 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_LGSE + INTEGER :: JGS_GSE_METHOD + REAL(8) :: JGS_GSE_TS INTEGER :: nbSel INTEGER :: UNSTSCHEMES(6) INTEGER :: UNSTSCHEME - 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. + 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 ! CHARACTER :: UGOBCFILE*60 REAL :: UGOBCDEPTH From afd8ae3615faf8e162dcfd1da748877f7fa8fc47 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Sat, 3 Aug 2024 09:29:04 +0000 Subject: [PATCH 06/10] ww3_gse2: fix severall bugs and improve scheme for filtered points --- model/src/w3gridmd.F90 | 13 +++++++-- model/src/w3iogrmd.F90 | 9 ++++-- model/src/w3profsmd_pdlib.F90 | 53 ++++++++++++++++++++++------------- model/src/w3wavemd.F90 | 4 +++ 4 files changed, 54 insertions(+), 25 deletions(-) diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90 index 5e6fb78bc..f5c6aa6b6 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -2467,8 +2467,9 @@ SUBROUTINE W3GRID() JGS_NLEVEL = 0 JGS_SOURCE_NONLINEAR = .FALSE. ! read data from the unstructured devoted namelist - CALL READNL ( NDSS, 'UNST', STATUS ) + CALL READNL ( NDSS, 'UNST', STATUS ) + B_JGS_USE_JACOBI = JGS_USE_JACOBI B_JGS_TERMINATE_MAXITER = JGS_TERMINATE_MAXITER B_JGS_TERMINATE_DIFFERENCE = JGS_TERMINATE_DIFFERENCE @@ -3346,7 +3347,10 @@ SUBROUTINE W3GRID() JGS_DIFF_THR, & JGS_NORM_THR, & JGS_NLEVEL, & - JGS_SOURCE_NONLINEAR + JGS_SOURCE_NONLINEAR, & + JGS_LGSE, & + JGS_GSE_METHOD, & + JGS_GSE_TS ! WRITE (NDSO,2976) P2SF, I1P2SF, I2P2SF, & US3D, I1US3D, I2US3D, & @@ -6672,7 +6676,10 @@ SUBROUTINE W3GRID() ', JGS_DIFF_THR=', F8.3, & ', JGS_NORM_THR=', F8.3, & ', JGS_NLEVEL=', I3, & - ', JGS_SOURCE_NONLINEAR=', L3 / ) + ', JGS_SOURCE_NONLINEAR=', L3, & + ', JGS_LGSE=', L3, & + ', JGS_GSE_METHOD=', I3, & + ', JGS_GSE_TS=', F15.3/) ! 960 FORMAT (/' Miscellaneous ',A/ & ' --------------------------------------------------') diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90 index ce4403ba3..0c3e3ec94 100644 --- a/model/src/w3iogrmd.F90 +++ b/model/src/w3iogrmd.F90 @@ -796,7 +796,8 @@ 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 #ifdef W3_ASCII WRITE (NDSA,*) & 'FSN, FSPSI,FSFCT,FSNIMP,FSTOTALIMP,FSTOTALEXP, & @@ -830,7 +831,8 @@ 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 #endif !Init COUNTCON and IOBDP to zero, it needs to be set somewhere or !removed @@ -995,7 +997,8 @@ 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/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 29b1770ea..03b796db7 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -7722,7 +7722,7 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) INTEGER NB_ITER, iIter, ip_global REAL*8 DeltaTmax, eDeltaT, CLATSMN, DFAC, RFAC, eDiffNorm REAL*8 eNorm, DTquot, diffc, dcell, XWIND, DIFFTOT - REAL*8 DVDX(NPA), DVDY(NPA), DV2DXX(NPA), DV2DXY(NPA), KH, SWFAC + REAL*8 DVDX(NPA), DVDY(NPA), DV2DXX(NPA), DV2DXY(NPA), KH, SWFAC(NPA) REAL eRealA(1), eRealB(1) @@ -7747,23 +7747,26 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) DFAC = 1. END IF CLATMN = COS ( 70 * DERA ) - + DO ITH = 1, NTH DO IK = 1, NK ISP = ITH + (IK-1)*NTH DO JSEA = 1, NPA CALL INIT_GET_ISEA(ISEA, JSEA) - TFAC = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 ) - CGD = 0.5 * GRAV / SIG(IK) * IOBDP_LOC(JSEA) - DSS = ( CGD * (XFR-1.))**2 * DTME / 12.! * TFAC - DNN = ( CGD * DTH )**2 * DTME / 12.! * TFAC - DCELL = CGD / 10.0 ! -> CELLP needs probably redifinition ... - KH = WN(IK,ISEA)*DW(JSEA) - IF (DW(JSEA) .gt. 1000.) THEN - SWFAC = 1. + IF (DW(ISEA) .gt. 1000.) THEN + SWFAC(JSEA) = 1. ELSE - SWFAC = 0.!(1.-0.5*(1+(2*KH/SINH(MIN(20.d0,2*KH)))))**8 + SWFAC(JSEA) = 0.!(1.-0.5*(1+(2*KH/SINH(MIN(20.d0,2*KH)))))**8 ENDIF + IF (SWFAC(JSEA) .LT. TINY(1.)) THEN + DIFFVEC(:,JSEA) = 0. + CYCLE + ENDIF + CGD = 0.5 * GRAV / SIG(IK) * IOBDP_LOC(JSEA) + DSS = ( CGD * (XFR-1.))**2 * DTME / 12. + DNN = ( CGD * DTH )**2 * DTME / 12. + DCELL = CGD / 10.0 ! -> CELLP needs probably redifinition ... + KH = WN(IK,ISEA)*DW(ISEA) XWIND = 3.3 * U10(ISEA)*WN(IK,ISEA)/SIG(IK) - 2.3 XWIND = MAX ( 0. , MIN ( 1. , XWIND ) ) #ifdef W3_XW0 @@ -7773,16 +7776,16 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) XWIND = 1. #endif TFAC = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 ) - DSS = SWFAC * (XWIND * DCELL + (1.-XWIND) * DSS * TFAC) + DSS = SWFAC(JSEA) * (XWIND * DCELL + (1.-XWIND) * DSS * TFAC) #ifdef W3_DSS0 DSS = 0. #endif - DNN = SWFAC * (XWIND * DCELL + (1.-XWIND) * DNN * TFAC) + DNN = SWFAC(JSEA) * (XWIND * DCELL + (1.-XWIND) * DNN * TFAC) DIFFVEC(1,JSEA) = (DSS*ECOS(ITH)**2+DNN*ESIN(ITH)**2) DIFFVEC(2,JSEA) = (DSS*ESIN(ITH)**2+DNN*ECOS(ITH)**2) / CLATS(ISEA)**2 DIFFVEC(3,JSEA) = ((DSS-DNN) * ESIN(ITH)*ECOS(ITH)) / CLATS(ISEA) - !write(3000+myrank,'(I10,20F20.10)') IK, SIG(IK), CGD, CLATS(ISEA), DSS, DNN, DIFFVEC(1,JSEA), DIFFVEC(2,JSEA) + !IF (DW(ISEA) .gt. 1000.) write(3000+myrank,'(I10,20F20.10)') IK, SIG(IK), CGD, CLATS(ISEA), DSS, DNN, DIFFVEC(1,JSEA), DIFFVEC(2,JSEA), DCELL, XWIND END DO CALL DIFFERENTIATE_XYDIR(DBLE(VA(ISP,:)),DVDX,DVDY) ! AR: this two lines can be optimized, however seems fast enough by now. @@ -7791,7 +7794,10 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) DeltaTmax = 1./TINY(1.) DO IE=1,NE NI = INE(:,IE) - eDet = 2. * PDLIB_TRIA(IE) + IF (SUM(SWFAC(NI)) .EQ. 0) THEN + CYCLE + ENDIF + eDet = 2. * PDLIB_TRIA(IE) DO IDX=1,3 V(1) = PDLIB_IEN(2*IDX-1,IE) V(2) = PDLIB_IEN(2*IDX ,IE) @@ -7819,12 +7825,15 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) DT_DIFF = DTG/NB_ITER PHI_V = 0. - !WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax + WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax !CALL FLUSH(5000+myrank) DO IT = 1, NB_ITER DO IE = 1, NE NI = INE(:,IE) + IF (SUM(SWFAC(NI)) .EQ. 0) THEN + CYCLE + ENDIF eDet = 2. * PDLIB_TRIA(IE) DEDX(1) = PDLIB_IEN(1,IE) DEDX(2) = PDLIB_IEN(3,IE) @@ -7847,7 +7856,11 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) END DO CALL PDLIB_exchange1DREAL(PHI_V) DO JSEA =1, NSEAL - VA(ISP,JSEA) = MAX(0.,VA(ISP,JSEA) - (DT_DIFF * PHI_V(JSEA) / PDLIB_SI(JSEA) + 2 * DV2DXY(JSEA) * DIFFVEC(3,JSEA)) * DFAC * IOBDP_LOC(JSEA)) + DIFFTOT = - (DT_DIFF * PHI_V(JSEA) / PDLIB_SI(JSEA) + DT_DIFF * 2 * DV2DXY(JSEA) * DIFFVEC(3,JSEA)) * DFAC * IOBDP_LOC(JSEA) + VA(ISP,JSEA) = MAX(0.,VA(ISP,JSEA) + DIFFTOT ) + !IF (ABS(DIFFTOT) .GT. 0.) THEN + ! WRITE(50000+myrank,'(2I10,10F20.10)') JSEA, ISP, VA(ISP,JSEA), DT_DIFF, PHI_V(JSEA), PDLIB_SI(JSEA), DV2DXY(JSEA), DIFFVEC(3,JSEA), DFAC, IOBDP_LOC(JSEA) + !ENDIF END DO END DO END DO @@ -7865,6 +7878,7 @@ SUBROUTINE DIFFERENTIATE_XYDIR(VAR, DVDX, DVDY) USE yowExchangeModule, only : PDLIB_exchange1DREAL USE yowNodepool, only : PDLIB_IEN, PDLIB_TRIA, NP, NPA USE yowElementpool, only : NE, INE + USE W3PARALL, only: INIT_GET_ISEA IMPLICIT NONE REAL*8, INTENT(IN) :: VAR(NPA) REAL*8, INTENT(INOUT) :: DVDX(NPA), DVDY(NPA) @@ -7873,7 +7887,7 @@ SUBROUTINE DIFFERENTIATE_XYDIR(VAR, DVDX, DVDY) REAL*8 :: DVDXIE, DVDYIE REAL*8 :: WEI(NPA) INTEGER :: NI(3) - INTEGER :: IE, JSEA + INTEGER :: IE, JSEA, ISEA WEI(:) = 0.d0 DVDX(:) = 0.d0 @@ -7905,7 +7919,8 @@ SUBROUTINE DIFFERENTIATE_XYDIR(VAR, DVDX, DVDY) DVDY = DVDY/WEI DO JSEA = 1, NP - IF (DW(JSEA) .LT. DMIN) THEN + CALL INIT_GET_ISEA(ISEA, JSEA) + IF (DW(ISEA) .LT. DMIN) THEN DVDX(JSEA) = 0. DVDY(JSEA) = 0. END IF diff --git a/model/src/w3wavemd.F90 b/model/src/w3wavemd.F90 index e6e106c39..46577c960 100644 --- a/model/src/w3wavemd.F90 +++ b/model/src/w3wavemd.F90 @@ -1812,6 +1812,10 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & END IF IF (LPDLIB) THEN + + !WRITE(*,*) B_JGS_LGSE, FSTOTALIMP, FSTOTALEXP + !STOP 'TEST TEST TEST' + ! #ifdef W3_PDLIB IF (FLCX .or. FLCY) THEN From 5aae07d1ff8a5b8aea724563ad2d6877ae19bc9d Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Sat, 3 Aug 2024 14:33:48 +0000 Subject: [PATCH 07/10] ww3_gse2: code cleaning and consolidation --- model/src/w3profsmd_pdlib.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 03b796db7..99a2cd037 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -7805,7 +7805,7 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) IP = INE(IDX,IE) eDiffNorm = HYPOT(DIFFVEC(1,IP), DIFFVEC(2,IP)) if (eDiffNorm .gt. tiny(1.)) then - eDeltaT = min(DTG, PDLIB_SI(IP) / (eDiffNorm * dfac * eNorm /(4. * eDet))) ! modified for diffusion vector + eDeltaT = min(DTG, 0.8 * PDLIB_SI(IP) / (eDiffNorm * dfac * eNorm /(4. * eDet))) ! modified for diffusion vector else eDeltaT = DTG endif @@ -7825,7 +7825,7 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) DT_DIFF = DTG/NB_ITER PHI_V = 0. - WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax + !WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax !CALL FLUSH(5000+myrank) DO IT = 1, NB_ITER @@ -7851,12 +7851,16 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) V(2) = 0.5 * PDLIB_IEN(2*IDX ,IE) eScal = DOT_PRODUCT(V, GRAD(1:2)) IP = INE(IDX,IE) - PHI_V(IP) = PHI_V(IP) + eScal + PHI_V(IP) = PHI_V(IP) + eScal + 1./3. * 2 * DV2DXY(IP) * DIFFVEC(3,IP) END DO END DO CALL PDLIB_exchange1DREAL(PHI_V) DO JSEA =1, NSEAL - DIFFTOT = - (DT_DIFF * PHI_V(JSEA) / PDLIB_SI(JSEA) + DT_DIFF * 2 * DV2DXY(JSEA) * DIFFVEC(3,JSEA)) * DFAC * IOBDP_LOC(JSEA) + IF (IOBP_LOC(JSEA) .EQ. 1) THEN + DIFFTOT = - DT_DIFF * PHI_V(JSEA) / PDLIB_SI(JSEA) * DFAC + ELSE + DIFFTOT = 0 + ENDIF VA(ISP,JSEA) = MAX(0.,VA(ISP,JSEA) + DIFFTOT ) !IF (ABS(DIFFTOT) .GT. 0.) THEN ! WRITE(50000+myrank,'(2I10,10F20.10)') JSEA, ISP, VA(ISP,JSEA), DT_DIFF, PHI_V(JSEA), PDLIB_SI(JSEA), DV2DXY(JSEA), DIFFVEC(3,JSEA), DFAC, IOBDP_LOC(JSEA) From 27ef1494e6111cd70930a1cb2fcf70945fa622e2 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Fri, 9 Aug 2024 08:48:34 +0000 Subject: [PATCH 08/10] ww3_gse2: add my latest work on gse ... --- model/src/w3profsmd_pdlib.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 99a2cd037..d0332864b 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -7825,8 +7825,8 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) DT_DIFF = DTG/NB_ITER PHI_V = 0. - !WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax - !CALL FLUSH(5000+myrank) + WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax + CALL FLUSH(5000+myrank) DO IT = 1, NB_ITER DO IE = 1, NE @@ -7851,13 +7851,13 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) V(2) = 0.5 * PDLIB_IEN(2*IDX ,IE) eScal = DOT_PRODUCT(V, GRAD(1:2)) IP = INE(IDX,IE) - PHI_V(IP) = PHI_V(IP) + eScal + 1./3. * 2 * DV2DXY(IP) * DIFFVEC(3,IP) + PHI_V(IP) = PHI_V(IP) + eScal !+ 1./3. * 2 * DV2DXY(IP) * DIFFVEC(3,IP) END DO END DO CALL PDLIB_exchange1DREAL(PHI_V) DO JSEA =1, NSEAL IF (IOBP_LOC(JSEA) .EQ. 1) THEN - DIFFTOT = - DT_DIFF * PHI_V(JSEA) / PDLIB_SI(JSEA) * DFAC + DIFFTOT = - DT_DIFF * DFAC * ( PHI_V(JSEA) / PDLIB_SI(JSEA) + 2 * DV2DXY(IP) * DIFFVEC(3,IP) ) ELSE DIFFTOT = 0 ENDIF From d28cbc306feccd4657a343120fbdc99e2dc245b5 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Fri, 9 Aug 2024 10:49:47 +0200 Subject: [PATCH 09/10] ww3_gse2: try something ... --- model/src/w3pro2md.F90 | 16 +++++++++ model/src/w3profsmd_pdlib.F90 | 68 ++++++++++++++++++++++++++++------- 2 files changed, 71 insertions(+), 13 deletions(-) diff --git a/model/src/w3pro2md.F90 b/model/src/w3pro2md.F90 index a23f893ef..4ab6a3e69 100644 --- a/model/src/w3pro2md.F90 +++ b/model/src/w3pro2md.F90 @@ -967,6 +967,22 @@ SUBROUTINE W3XYP2 ( ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY ) #endif DNN = XWIND * DCELL + (1.-XWIND) * DNND * TFAC + + CGD = 0.5 * GRAV / SIG(IK) * IOBDP_LOC(JSEA) + DSS = ( CGD * (XFR-1.))**2 * DTME / 12. + DNN = ( CGD * DTH )**2 * DTME / 12. + DCELL = CGD / 10.0 ! -> CELLP needs probably redifinition ... + KH = WN(IK,ISEA)*DW(ISEA) + XWIND = 3.3 * U10(ISEA)*WN(IK,ISEA)/SIG(IK) - 2.3 + XWIND = MAX ( 0. , MIN ( 1. , XWIND ) ) + TFAC = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 ) + DSS = SWFAC(JSEA) * (XWIND * DCELL + (1.-XWIND) * DSS * TFAC) + DNN = SWFAC(JSEA) * (XWIND * DCELL + (1.-XWIND) * DNN * TFAC) + DIFFVEC(1,JSEA) = (DSS*ECOS(ITH)**2+DNN*ESIN(ITH)**2) + DIFFVEC(2,JSEA) = (DSS*ESIN(ITH)**2+DNN*ECOS(ITH)**2) / CLATS(ISEA)**2 + DIFFVEC(3,JSEA) = ((DSS-DNN) * ESIN(ITH)*ECOS(ITH)) / CLATS(ISEA) + + VDXX(IXY) = DTLOC * (DSS*ECOS(ITH)**2+DNN*ESIN(ITH)**2) VDYY(IXY) = DTLOC * (DSS*ESIN(ITH)**2+DNN*ECOS(ITH)**2) & / CLATS(ISEA)**2 diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 03b796db7..6fcd3bba0 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -676,7 +676,7 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) ! make the interface between the WAVEWATCH and the WWM code. ! ! 8. Structure : - ! + !PDLIB_W3XYPUG ! ! 9. Switches : ! @@ -718,17 +718,20 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ - INTEGER :: ITH, IK, ISEA + INTEGER :: ITH, IK, ISEA, ITHL, ITHR INTEGER :: I, J, IE, IBND_MAP INTEGER :: IP_glob REAL :: CCOS, CSIN, CCURX, CCURY, WN1, CG1 - REAL :: C(npa,2) + REAL :: CCOSL, CSINL, CCOSR, CSINR + REAL :: C(npa,2), CL(npa,2), CR(npa,2) REAL :: RD1, RD2 !/ !/ Automatic work arrays !/ REAL :: VLCFLX(npa), VLCFLY(npa) - REAL :: AC(npa) + REAL :: VLCFLXL(npa), VLCFLYL(npa) + REAL :: VLCFLXR(npa), VLCFLYR(npa) + REAL :: AC(npa), ACL(npa), ACM(npa), ACR(npa) REAL :: AC_MAP(NBND_MAP) INTEGER :: JSEA, IP !/ ------------------------------------------------------------------- / @@ -745,8 +748,29 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) #endif ITH = 1 + MOD(ISP-1,NTH) IK = 1 + (ISP-1)/NTH + + IF (ITH == 1) THEN ! Compute indices for computation of ith,isp for gse + ITHL = NTH + ELSE IF (ITH == NTH) THEN + ITHR = 1 + ELSE + ITHL = ITH - 1 + ITHR = ITH + 1 + ENDIF + CCOS = FACX * ECOS(ITH) CSIN = FACY * ESIN(ITH) + + CCOSL = FACX * ECOS(ITHL) + CSINL = FACY * ESIN(ITHL) + CCOSR = FACX * ECOS(ITHR) + CSINR = FACY * ESIN(ITHR) + + CCOSL = 0.5 * (CCOSL + CCOS) + CSINL = 0.5 * (CSINL + CSIN) + CCOSR = 0.5 * (CCOSR + CCOS) + CSINR = 0.5 * (CSINR + CSIN) + CCURX = FACX CCURY = FACY ! @@ -771,6 +795,11 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) AC(IP) = VA(ISP,JSEA) / CG(IK,ISEA) * CLATS(ISEA) VLCFLX(IP) = CCOS * CG(IK,ISEA) / CLATS(ISEA) VLCFLY(IP) = CSIN * CG(IK,ISEA) + + VLCFLXL(IP) = CCOSL * CG(IK,ISEA) / CLATS(ISEA) + VLCFLYL(IP) = CSINR * CG(IK,ISEA) + VLCFLXR(IP) = CCOSL * CG(IK,ISEA) / CLATS(ISEA) + VLCFLYR(IP) = CSINR * CG(IK,ISEA) #endif #ifdef W3_MGP VLCFLX(IP) = VLCFLX(IP) - CCURX*VGX/CLATS(ISEA) @@ -802,6 +831,12 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) C(:,1) = VLCFLX(:) * IOBDP_LOC C(:,2) = VLCFLY(:) * IOBDP_LOC + + CL(:,1) = VLCFLXL(:) * IOBDP_LOC + CL(:,2) = VLCFLYL(:) * IOBDP_LOC + + CR(:,1) = VLCFLXR(:) * IOBDP_LOC + CR(:,2) = VLCFLYR(:) * IOBDP_LOC ! ! 4. Prepares boundary update ! @@ -833,7 +868,10 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) ELSE IF (FSPSI) THEN CALL PDLIB_W3XYPFSPSI2(ISP, C, LCALC, RD1, RD2, DTG, AC) ELSE IF (FSFCT) THEN - CALL PDLIB_W3XYPFSFCT2(ISP, C, LCALC, RD1, RD2, DTG, AC) + CALL PDLIB_W3XYPFSFCT2(ISP, CL, LCALC, RD1, RD2, DTG, 0.25 * AC, ACL) + CALL PDLIB_W3XYPFSFCT2(ISP, C, LCALC, RD1, RD2, DTG, 0.5 * AC, ACM) + CALL PDLIB_W3XYPFSFCT2(ISP, CR, LCALC, RD1, RD2, DTG, 0.25 * AC, ACR) + AC = ACL + ACM + ACR ELSE IF (FSNIMP) THEN STOP 'For PDLIB and FSNIMP, no function has been programmed yet' ENDIF @@ -1492,7 +1530,7 @@ SUBROUTINE PDLIB_W3XYPFSPSI2 ( ISP, C, LCALC, RD10, RD20, DT, AC) END SUBROUTINE PDLIB_W3XYPFSPSI2 !/ ------------------------------------------------------------------- / - SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) + SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, ACIN, ACOUT) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | @@ -1568,7 +1606,9 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) ! for the given velocity field REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's ! X- and Y- Components, - REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and + REAL, INTENT(IN) :: ACIN(npa) ! Wave Action before and + ! after advection + REAL, INTENT(OUT) :: ACOUT(npa) ! Wave Action before and ! after advection REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation ! coefficients for boundary @@ -1616,7 +1656,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) #ifdef W3_DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 1' FLUSH(740+IAPROC) - CALL SCAL_INTEGRAL_PRINT_R4(AC, "AC in input") + CALL SCAL_INTEGRAL_PRINT_R4(ACIN, "AC in input") #endif ITH = 1 + MOD(ISP-1,NTH) @@ -1699,9 +1739,11 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) DTSI(IP) = DBLE(DT)/DBLE(ITER(IK,ITH))/PDLIB_SI(IP) ! Some precalculations for the time integration. END DO + ACOUT = ACIN + DO IT = 1, ITER(IK,ITH) - U = DBLE(AC) + U = DBLE(ACOUT) ST = ZERO PM = ZERO PP = ZERO @@ -1791,7 +1833,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) #endif END DO - AC = REAL(U) + ACOUT = REAL(U) #ifdef W3_DEBUGSOLVER IF (testWrite) THEN @@ -1823,10 +1865,10 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) IP_glob = MAPSF(ISBPI(IBI),1) JX=IPGL_npa(IP_glob) IF (JX .gt. 0) THEN - AC(JX) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) ) & + ACOUT(JX) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) ) & / CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI)) #ifdef W3_DEBUGSOLVER - sumAC=sumAC + AC(JX) + sumAC=sumAC + ACOUT(JX) sumBPI0=sumBPI0 + BBPI0(ISP,IBI) sumBPIN=sumBPIN + BBPIN(ISP,IBI) sumCG=sumCG + CG(IK,ISBPI(IBI)) @@ -1846,7 +1888,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumCLATS=', sumCLATS FLUSH(740+IAPROC) #endif - CALL PDLIB_exchange1DREAL(AC) + CALL PDLIB_exchange1DREAL(ACOUT) #ifdef W3_DEBUGSOLVER IF (testWrite) THEN From 3aff39b0c1e535de34c5799241d155976eee96d8 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Wed, 28 Aug 2024 17:12:13 +0000 Subject: [PATCH 10/10] ww3_gse2: some checks due to grid related issues ... to be continued --- model/src/w3gridmd.F90 | 4 ++-- model/src/w3profsmd_pdlib.F90 | 26 ++++++++++++++------------ model/src/w3triamd.F90 | 3 +++ 3 files changed, 19 insertions(+), 14 deletions(-) diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90 index f5c6aa6b6..b72ed8e8e 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -7153,7 +7153,7 @@ SUBROUTINE W3GRID() ' Number of longitudes :',I10/ & ' Number of latitudes :',I10/ & ' Number of grid points :',I10/ & - ' Number of sea points :',I10,' (',F4.1,'%)'/& + ' Number of sea points :',I10,' (',F15.4,'%)'/& ' Number of input b. points :',I10/ & ' Number of land points :',I10/ & ' Number of excluded points :',I10/) @@ -7162,7 +7162,7 @@ SUBROUTINE W3GRID() ' Number of longitudes :',I10/ & ' Number of latitudes :',I10/ & ' Number of grid points :',I10/ & - ' Number of sea points :',I10,' (100%)'/ & + ' Number of sea points :',I10,' (test)'/ & ' Number of input b. points :',I10/ & ' Number of land points :',I10/ & ' Number of excluded points :',I10/) diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 91b0f09e2..0af0e95b5 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -751,7 +751,9 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) IF (ITH == 1) THEN ! Compute indices for computation of ith,isp for gse ITHL = NTH + ITHR = 2 ELSE IF (ITH == NTH) THEN + ITHL = NTH - 1 ITHR = 1 ELSE ITHL = ITH - 1 @@ -766,11 +768,6 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) CCOSR = FACX * ECOS(ITHR) CSINR = FACY * ESIN(ITHR) - CCOSL = 0.5 * (CCOSL + CCOS) - CSINL = 0.5 * (CSINL + CSIN) - CCOSR = 0.5 * (CCOSR + CCOS) - CSINR = 0.5 * (CSINR + CSIN) - CCURX = FACX CCURY = FACY ! @@ -837,6 +834,10 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) CR(:,1) = VLCFLXR(:) * IOBDP_LOC CR(:,2) = VLCFLYR(:) * IOBDP_LOC + + CL = 0.5 * (C + CL) + CR = 0.5 * (C + CR) + ! ! 4. Prepares boundary update ! @@ -868,10 +869,11 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC ) ELSE IF (FSPSI) THEN CALL PDLIB_W3XYPFSPSI2(ISP, C, LCALC, RD1, RD2, DTG, AC) ELSE IF (FSFCT) THEN - CALL PDLIB_W3XYPFSFCT2(ISP, CL, LCALC, RD1, RD2, DTG, 0.25 * AC, ACL) - CALL PDLIB_W3XYPFSFCT2(ISP, C, LCALC, RD1, RD2, DTG, 0.5 * AC, ACM) - CALL PDLIB_W3XYPFSFCT2(ISP, CR, LCALC, RD1, RD2, DTG, 0.25 * AC, ACR) - AC = ACL + ACM + ACR + !CALL PDLIB_W3XYPFSFCT2(ISP, CL, LCALC, RD1, RD2, DTG, 0.25 * AC, ACL) + !CALL PDLIB_W3XYPFSFCT2(ISP, C, LCALC, RD1, RD2, DTG, 0.5 * AC, ACM) + !CALL PDLIB_W3XYPFSFCT2(ISP, CR, LCALC, RD1, RD2, DTG, 0.25 * AC, ACR) + !AC = ACL + ACM + ACR + CALL PDLIB_W3XYPFSFCT2(ISP, C, LCALC, RD1, RD2, DTG, AC, AC) ELSE IF (FSNIMP) THEN STOP 'For PDLIB and FSNIMP, no function has been programmed yet' ENDIF @@ -7867,8 +7869,8 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) DT_DIFF = DTG/NB_ITER PHI_V = 0. - WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax - CALL FLUSH(5000+myrank) + !WRITE(5000+myrank,*) 'NUMBER OF SUB ITERATIONS', ITH, IK, NB_ITER, DT_DIFF, DeltaTmax + !CALL FLUSH(5000+myrank) DO IT = 1, NB_ITER DO IE = 1, NE @@ -7899,7 +7901,7 @@ SUBROUTINE BLOCK_SOLVER_DIFFUSION(DTG) CALL PDLIB_exchange1DREAL(PHI_V) DO JSEA =1, NSEAL IF (IOBP_LOC(JSEA) .EQ. 1) THEN - DIFFTOT = - DT_DIFF * DFAC * ( PHI_V(JSEA) / PDLIB_SI(JSEA) + 2 * DV2DXY(IP) * DIFFVEC(3,IP) ) + DIFFTOT = - DT_DIFF * DFAC * ( PHI_V(JSEA) / PDLIB_SI(JSEA) + 2 * DV2DXY(JSEA) * DIFFVEC(3,JSEA) ) ELSE DIFFTOT = 0 ENDIF diff --git a/model/src/w3triamd.F90 b/model/src/w3triamd.F90 index 9fac503b6..ff7e5d86c 100644 --- a/model/src/w3triamd.F90 +++ b/model/src/w3triamd.F90 @@ -249,6 +249,7 @@ SUBROUTINE READMSH(NDS,FNAME) DO I= 1, NODES READ(NDS,*) j, XYBTMP1(1,I), XYBTMP1(2,I), XYBTMP1(3,I) END DO + WRITE(*,*) 'TEST NUMBER OF NODES', NODES ! ! read number of elements and elements from Gmsh files ! @@ -321,6 +322,8 @@ SUBROUTINE READMSH(NDS,FNAME) ! ! Number of nodes after clean up ! + WRITE(*,*) 'TEST NUMBER OF NODES AFTER CLEAN UP', J + NX = J ! DO I = 1, NTRI