From c6aeea3e283bdcde087446b9af657327fd1bbd08 Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Fri, 7 Apr 2023 12:05:48 -0400 Subject: [PATCH 1/3] fix for block_explicit * use npa when filling va at the end of pdlib_explict_block * re-indent w3wavemd after removal of if-block (from wise pr) * trailing whitespace cleanup (from wise pr) --- model/src/w3iorsmd.F90 | 1 - model/src/w3profsmd_pdlib.F90 | 278 +++++++++++++++++----------------- model/src/w3wavemd.F90 | 104 ++++++------- 3 files changed, 191 insertions(+), 192 deletions(-) diff --git a/model/src/w3iorsmd.F90 b/model/src/w3iorsmd.F90 index 76c2a5cb8..a969aa9af 100644 --- a/model/src/w3iorsmd.F90 +++ b/model/src/w3iorsmd.F90 @@ -769,7 +769,6 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT ) ! Include remainder values (switch to record format) ---- * JSEA = NSEAL_MIN + 1 IF ( JSEA.EQ.NSEAL ) THEN - !ISEA = IAPROC + (JSEA - 1) * NAPROC CALL INIT_GET_ISEA(ISEA, JSEA) NREC = ISEA + 2 RPOS = 1_8 + LRECL*(NREC-1_8) diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 144ace494..d29cb3f15 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -125,7 +125,7 @@ MODULE PDLIB_W3PROFSMD INTEGER :: FreqShiftMethod = 2 LOGICAL :: FSGEOADVECT LOGICAL, SAVE :: LINIT_OUTPUT = .TRUE. - REAL, SAVE :: RTIME = 0.d0 + REAL, SAVE :: RTIME = 0.d0 INTEGER :: POS_TRICK(3,2) #ifdef W3_DEBUGSRC @@ -245,7 +245,7 @@ SUBROUTINE PDLIB_INIT(IMOD) WRITE(740+IAPROC,*) 'NTPROC=', NTPROC FLUSH(740+IAPROC) #endif - + PDLIB_NSEAL = 0 IF (IAPROC .le. NAPROC) THEN @@ -263,7 +263,7 @@ SUBROUTINE PDLIB_INIT(IMOD) CALL initFromGridDim(NX,NTRI,TRIGP,NTH,MPI_COMM_WCMP) ELSE CALL initFromGridDim(NX,NTRI,TRIGP,NSPEC,MPI_COMM_WCMP) - ENDIF + ENDIF ! #ifdef W3_DEBUGSOLVER WRITE(740+IAPROC,*) 'After initFromGridDim' @@ -936,19 +936,19 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) USE W3PARALL, only : ONESIXTH, ZERO, THR USE yowRankModule, only : IPGL_npa - INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, + INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, ! actual Wave Direction - REAL, INTENT(IN) :: DT ! Time intervall for which the + REAL, INTENT(IN) :: DT ! Time intervall for which the ! advection should be computed ! for the given velocity field - REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's + 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(INOUT) :: AC(npa) ! Wave Action before and ! after advection - REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation + REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation ! coefficients for boundary ! conditions - LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of + LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of ! the max. Global Time step #ifdef W3_S INTEGER, SAVE :: IENT = 0 @@ -1068,7 +1068,7 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) ITER(IK,ITH) = ABS(NINT(CFLXY)) END IF END IF ! LCALC - + #ifdef W3_DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 4' FLUSH(740+IAPROC) @@ -1089,7 +1089,7 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) FLUSH(740+IAPROC) #endif - + DO IT = 1, ITER(IK,ITH) #ifdef W3_DEBUGSOLVER WRITE(740+IAPROC,*) 'IK=', IK, ' ITH=', ITH @@ -1264,19 +1264,19 @@ SUBROUTINE PDLIB_W3XYPFSPSI2 ( ISP, C, LCALC, RD10, RD20, DT, AC) USE W3PARALL, only : ONESIXTH, ZERO, THR USE yowRankModule, only : IPGL_npa IMPLICIT NONE - INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, + INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, ! actual Wave Direction - REAL, INTENT(IN) :: DT ! Time intervall for which the + REAL, INTENT(IN) :: DT ! Time intervall for which the ! advection should be computed ! for the given velocity field - REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's + 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(INOUT) :: AC(npa) ! Wave Action before and ! after advection - REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation + REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation ! coefficients for boundary ! conditions - LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of + LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of ! the max. Global Time step #ifdef W3_S INTEGER, SAVE :: IENT = 0 @@ -1561,19 +1561,19 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) USE yowRankModule, only : IPGL_npa IMPLICIT NONE - INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, + INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, ! actual Wave Direction - REAL, INTENT(IN) :: DT ! Time intervall for which the + REAL, INTENT(IN) :: DT ! Time intervall for which the ! advection should be computed ! for the given velocity field - REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's + 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(INOUT) :: AC(npa) ! Wave Action before and ! after advection - REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation + REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation ! coefficients for boundary ! conditions - LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of + LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of ! the max. Global Time step #ifdef W3_S INTEGER, SAVE :: IENT = 0 @@ -1595,7 +1595,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) REAL*8 :: FT, UTILDE REAL*8 :: FL11, FL12, FL21, FL22, FL31, FL32 REAL*8 :: FL111, FL112, FL211, FL212, FL311, FL312 - REAL :: DTSI(npa), U(npa), UL(npa) + REAL :: DTSI(npa), U(npa), UL(npa) REAL :: DTMAX_GL, DTMAX, DTMAXEXP, REST REAL*8 :: LAMBDA(2), KTMP(3) REAL*8 :: KELEM(3,NE), FLALL(3,NE) @@ -1724,7 +1724,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) END IF ST(NI) = ST(NI) + THETA_L(:,IE) ! the 2nd term are the theta values of each node ... THETA_H = (1./3.+DT/(2.*PDLIB_TRIA(IE)) * KELEM(:,IE) ) * FT ! LAX -! THETA_H = (1./3.+2./3.*KELEM(:,IE)/SUM(MAX(ZERO,KELEM(:,IE))))*FT ! CENTRAL ... can be tested as well a bit more dispersive then LAX + ! THETA_H = (1./3.+2./3.*KELEM(:,IE)/SUM(MAX(ZERO,KELEM(:,IE))))*FT ! CENTRAL ... can be tested as well a bit more dispersive then LAX THETA_ACE(:,IE) = THETA_H-THETA_L(:,IE) PP(NI) = PP(NI) + MAX(ZERO, -THETA_ACE(:,IE)) * DTSI(NI) PM(NI) = PM(NI) + MIN(ZERO, -THETA_ACE(:,IE)) * DTSI(NI) @@ -1745,7 +1745,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) #endif USTARI(1,:) = MAX(UL,U) - USTARI(2,:) = MIN(UL,U) + USTARI(2,:) = MIN(UL,U) UIP = 0. UIM = 0. @@ -1861,7 +1861,7 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) #endif END SUBROUTINE PDLIB_W3XYPFSFCT2 -!/ ------------------------------------------------------------------- / + !/ ------------------------------------------------------------------- / SUBROUTINE TEST_MPI_STATUS(string) !/ @@ -1939,7 +1939,7 @@ SUBROUTINE TEST_MPI_STATUS(string) WRITE(740+IAPROC,*) 'Leaving the TEST_MPI_STATUS' FLUSH(740+IAPROC) END SUBROUTINE TEST_MPI_STATUS -!/ ------------------------------------------------------------------- / + !/ ------------------------------------------------------------------- / SUBROUTINE SCAL_INTEGRAL_PRINT_GENERAL(V, string, maxidx, CheckUncovered, PrintFullValue) !/ @@ -1955,7 +1955,7 @@ SUBROUTINE SCAL_INTEGRAL_PRINT_GENERAL(V, string, maxidx, CheckUncovered, PrintF !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ - ! 1. Purpose : Source code for parallel debugging + ! 1. Purpose : Source code for parallel debugging ! 2. Method : maxidx = npa or np for arrays that have been synchronized or not ! CheckUncovered is because some the triangulation may not cover all nodes ! 3. Parameters : @@ -2785,7 +2785,7 @@ SUBROUTINE PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC ) USE W3ODATMD, only: IAPROC USE W3GDATMD, only: B_JGS_USE_JACOBI - LOGICAL, INTENT(IN) :: LCALC + LOGICAL, INTENT(IN) :: LCALC INTEGER, INTENT(IN) :: IMOD REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY #ifdef W3_DEBUGSOLVER @@ -2855,7 +2855,7 @@ SUBROUTINE PDLIB_W3XYPUG_BLOCK_EXPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) USE W3ODATMD, only: IAPROC USE W3GDATMD, only: B_JGS_USE_JACOBI - LOGICAL, INTENT(IN) :: LCALC + LOGICAL, INTENT(IN) :: LCALC INTEGER, INTENT(IN) :: IMOD REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY @@ -3664,7 +3664,7 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) CXY(1,IP) = CCOS * CG1/CLATS(IP_GLOB) CXY(2,IP) = CSIN * CG1 IF (FLCUR) THEN - CXY(1,IP) = CXY(1,IP) + FACX * CX(IP_GLOB)/CLATS(IP_GLOB)*IOBDP_LOC(IP) + CXY(1,IP) = CXY(1,IP) + FACX * CX(IP_GLOB)/CLATS(IP_GLOB)*IOBDP_LOC(IP) CXY(2,IP) = CXY(2,IP) + FACY * CY(IP_GLOB)*IOBDP_LOC(IP) ENDIF #ifdef W3_MGP @@ -3700,7 +3700,7 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) IB1 = (1-IOBPA_LOC(IP)) * IOBPD_LOC(ITH,IP) IB2 = IOBPD_LOC(ITH,IP) #ifdef W3_REF1 - IBR = (1-IOBP_LOC(IP)) * (1-IOBPD_LOC(ITH,IP)) * (1-IOBPA_LOC(IP)) + IBR = (1-IOBP_LOC(IP)) * (1-IOBPD_LOC(ITH,IP)) * (1-IOBPA_LOC(IP)) #endif IF (IOBDP_LOC(IP) .eq. 1) THEN DO I = 1, PDLIB_CCON(IP) @@ -3710,11 +3710,11 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) #ifdef W3_DEBUGSRC WRITE(740+IAPROC,*) 'I1=', I1, ' PDLIB_I_DIAG=', PDLIB_I_DIAG(IP) #endif - -#ifdef W3_REF1 + +#ifdef W3_REF1 IF (IBR == 1) THEN - DTK = KP(POS,IE) * DTG - B_JAC(ISP,IP) = B_JAC(ISP,IP) + PDLIB_TRIA03(IE) * VA(ISP,IP) + DTK = KP(POS,IE) * DTG + B_JAC(ISP,IP) = B_JAC(ISP,IP) + PDLIB_TRIA03(IE) * VA(ISP,IP) ELSE DTK = KP(POS,IE) * DTG * IB1 B_JAC(ISP,IP) = B_JAC(ISP,IP) + PDLIB_TRIA03(IE) * VA(ISP,IP) * IB2 @@ -4467,7 +4467,7 @@ SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1(DTG) INTEGER :: ITH0 LOGICAL :: LSIG = .FALSE. - !AR: TODO: check&report if needed ... + !AR: TODO: check&report if needed ... LSIG = FLCUR .OR. FLLEV DO IP = 1, np @@ -4495,8 +4495,8 @@ SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1(DTG) END IF CAS_SIG(:,IP) = CAS ELSE IF (FreqShiftMethod .eq. 2) THEN - IF (IOBP_LOC(IP).eq.1.and.IOBDP_LOC(IP).eq.1.and.IOBPA_LOC(IP).eq.0) THEN - CALL PROP_FREQ_SHIFT_M2(IP, ISEA, CWNB_M2, DWNI_M2, DTG) + IF (IOBP_LOC(IP).eq.1.and.IOBDP_LOC(IP).eq.1.and.IOBPA_LOC(IP).eq.0) THEN + CALL PROP_FREQ_SHIFT_M2(IP, ISEA, CWNB_M2, DWNI_M2, DTG) #ifdef W3_DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'sum(CWNB_M2)=', sum(CWNB_M2) #endif @@ -5534,10 +5534,10 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL USE W3SRC4MD, only: W3SPR4 #endif #ifdef W3_REF1 - USE W3GDATMD, only: REFPARS + USE W3GDATMD, only: REFPARS #endif implicit none - LOGICAL, INTENT(IN) :: LCALC + LOGICAL, INTENT(IN) :: LCALC INTEGER, INTENT(IN) :: IMOD REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY ! @@ -5709,7 +5709,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL ! ! geographical advection ! - IF (IMEM == 1) THEN + IF (IMEM == 1) THEN call calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) ENDIF @@ -5974,10 +5974,10 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL DO IK=1,NK DO ITH=1,NTH ISP = ITH + (IK-1)*NTH - IF (REFPARS(3) .LT. 0.5 .AND. IOBPD_LOC(ITH,IP) .EQ. 0 .AND. IOBPA_LOC(IP) .EQ. 0) THEN + IF (REFPARS(3) .LT. 0.5 .AND. IOBPD_LOC(ITH,IP) .EQ. 0 .AND. IOBPA_LOC(IP) .EQ. 0) THEN VA(ISP,IP) = VAOLD(ISP,IP) * IOBDP_LOC(IP) ! Restores reflected action spectra ... ENDIF - ENDDO + ENDDO ENDDO #endif ELSE @@ -6045,7 +6045,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL ! ! Terminate via differences ! - IF (B_JGS_TERMINATE_DIFFERENCE .and. INT(MOD(NBITER,10)) == 0) THEN ! Every 10th step check conv. + IF (B_JGS_TERMINATE_DIFFERENCE .and. INT(MOD(NBITER,10)) == 0) THEN ! Every 10th step check conv. CALL MPI_ALLREDUCE(is_converged, itmp, 1, MPI_INT, MPI_SUM, MPI_COMM_WCMP, ierr) is_converged = itmp prop_conv = (DBLE(NX) - DBLE(is_converged))/DBLE(NX) * 100. @@ -6160,7 +6160,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL VA(ISP,IP)=MAX(ZERO, VA(ISP,IP))*IOBDP_LOC(IP)*DBLE(IOBPD_LOC(ITH,IP)) #ifdef W3_REF1 IF (REFPARS(3).LT.0.5.AND.IOBPD_LOC(ITH,IP).EQ.0.AND.IOBPA_LOC(IP).EQ.0) THEN - VA(ISP,IP) = VAOLD(ISP,IP) ! restores reflected boundary values + VA(ISP,IP) = VAOLD(ISP,IP) ! restores reflected boundary values ENDIF #endif END DO @@ -6238,56 +6238,56 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL #endif IF (FLSOU) THEN - IF (B_JGS_LIMITER) THEN + IF (B_JGS_LIMITER) THEN - DO ISP=1,NSPEC - IK = 1 + (ISP-1)/NTH - SPEC(ISP) = VAOLD(ISP,JSEA) - ENDDO + 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), & + 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) - END DO -! - DO IK=1, NK - IS0 = (IK-1)*NTH - DO ITH=2, NTH - DAM(ITH+IS0) = DAM(1+IS0) - 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) - 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 + 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) + END DO + ! + DO IK=1, NK + IS0 = (IK-1)*NTH + DO ITH=2, NTH + DAM(ITH+IS0) = DAM(1+IS0) + 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) + 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 ENDIF ! B_JGS_LIMITER ENDIF ! FLSOU END DO ! JSEA @@ -6343,9 +6343,9 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : Explicit block solver - ! 2. Method : It uses the n-scheme and the idea is to reduce latency due - ! to DD communication and increase vectorization level on the - ! single core + ! 2. Method : It uses the n-scheme and the idea is to reduce latency due + ! to DD communication and increase vectorization level on the + ! single core ! 3. Parameters : ! ! Parameter list @@ -6393,14 +6393,14 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) USE MPI, only : MPI_MIN #endif #ifdef W3_REF1 - USE W3GDATMD, only: REFPARS + USE W3GDATMD, only: REFPARS #endif - IMPLICIT NONE - - LOGICAL, INTENT(IN) :: LCALC + IMPLICIT NONE + + LOGICAL, INTENT(IN) :: LCALC - INTEGER, INTENT(IN) :: IMOD + INTEGER, INTENT(IN) :: IMOD REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY @@ -6412,7 +6412,7 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) REAL :: LAMBDAX(NTH), LAMBDAY(NTH) REAL :: DTMAX(NTH), DTMAXEXP(NTH), DTMAXOUT, DTMAXGL REAL :: FIN(1), FOUT(1), REST, CFLXY, RD1, RD2, RD10, RD20 - REAL :: UOLD(NTH,NPA), U(NTH,NPA) + REAL :: UOLD(NTH,NPA), U(NTH,NPA) REAL, PARAMETER :: ONESIXTH = 1.0/6.0 REAL, PARAMETER :: ZERO = 0.0 @@ -6425,32 +6425,32 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) ! 2. Calculate velocities ---------------- * ! ! 2a. Vectorized for all points looping over each wave number (maybe do a dirty save will be nice!) - ! + ! DO IK = 1, NK IF (LCALC) THEN - DO IP = 1, NPA - CALL WAVNU3 (SIG(IK), DW(iplg(IP)), KSIG(IP), CGSIG(IP)) - ENDDO + DO IP = 1, NPA + CALL WAVNU3 (SIG(IK), DW(iplg(IP)), KSIG(IP), CGSIG(IP)) + ENDDO DO ITH = 1, NTH DO IP = 1, NPA ISEA = IPLG(IP) CXX(ITH,IP) = CGSIG(IP) * FACX * ECOS(ITH) / CLATS(ISEA) CYY(ITH,IP) = CGSIG(IP) * FACY * ESIN(ITH) - ENDDO + ENDDO IF (FLCUR) THEN - DO IP = 1, NPA + DO IP = 1, NPA ISEA = IPLG(IP) IF (IOBP_LOC(IP) .GT. 0) THEN CXX(ITH,IP) = CXX(ITH,IP) + FACX * CX(ISEA)/CLATS(ISEA) CYY(ITH,IP) = CYY(ITH,IP) + FACY * CY(ISEA) ENDIF - ENDDO + ENDDO ENDIF - ENDDO + ENDDO DO IE = 1, NE @@ -6466,14 +6466,14 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) KELEM1(ITH,IE,IK) = LAMBDAX(ITH) * PDLIB_IEN(1,IE) + LAMBDAY(ITH) * PDLIB_IEN(2,IE) ! K-Values - so called Flux Jacobians KELEM2(ITH,IE,IK) = LAMBDAX(ITH) * PDLIB_IEN(3,IE) + LAMBDAY(ITH) * PDLIB_IEN(4,IE) KELEM3(ITH,IE,IK) = LAMBDAX(ITH) * PDLIB_IEN(5,IE) + LAMBDAY(ITH) * PDLIB_IEN(6,IE) - KTMP(1) = KELEM1(ITH,IE,IK) ! Extract + KTMP(1) = KELEM1(ITH,IE,IK) ! Extract KTMP(2) = KELEM2(ITH,IE,IK) KTMP(3) = KELEM3(ITH,IE,IK) NM(ITH,IE,IK) = - 1.D0/MIN(-THR,SUM(MIN(ZERO,KTMP))) ! N-Values KELEM1(ITH,IE,IK) = MAX(ZERO,KTMP(1)) KELEM2(ITH,IE,IK) = MAX(ZERO,KTMP(2)) KELEM3(ITH,IE,IK) = MAX(ZERO,KTMP(3)) - ENDDO + ENDDO FL11 = CXX(:,I2) * PDLIB_IEN(1,IE) + CYY(:,I2) * PDLIB_IEN(2,IE) ! Weights for Simpson Integration FL12 = CXX(:,I3) * PDLIB_IEN(1,IE) + CYY(:,I3) * PDLIB_IEN(2,IE) @@ -6498,14 +6498,14 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) KKSUM = ZERO DO IE = 1, NE NI = INE(:,IE) - DO ITH = 1, NTH + DO ITH = 1, NTH KKSUM(ITH,NI(1)) = KKSUM(ITH,NI(1)) + KELEM1(ITH,IE,IK) KKSUM(ITH,NI(2)) = KKSUM(ITH,NI(2)) + KELEM2(ITH,IE,IK) KKSUM(ITH,NI(3)) = KKSUM(ITH,NI(3)) + KELEM3(ITH,IE,IK) ENDDO END DO - DTMAXEXP = 1.E10 + DTMAXEXP = 1.E10 DTMAX = 1.E10 DO IP = 1, np IF (IOBP_LOC(IP) .EQ. 1 .OR. FSBCCFL) THEN @@ -6537,25 +6537,25 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) END IF ! LCALC -! Exact and convert Wave Action - should be some subroutine function or whatever + ! Exact and convert Wave Action - should be some subroutine function or whatever DO ITH = 1, NTH ISP = ITH + (IK-1) * NTH DO IP = 1, NPA U(ITH,IP) = VA(ISP,IP) / CGSIG(IP) * CLATS(IPLG(IP)) - ENDDO - ENDDO + ENDDO + ENDDO UOLD = U DO IT = 1, ITER(IK) ST = ZERO DO IE = 1, NE NI = INE(:,IE) - DO ITH = 1, NTH + DO ITH = 1, NTH UTILDE(ITH) = NM(ITH,IE,IK) * (FLALL1(ITH,IE,IK)*U(ITH,NI(1)) + FLALL2(ITH,IE,IK)*U(ITH,NI(2)) + FLALL3(ITH,IE,IK)*U(ITH,NI(3))) ST(ITH,NI(1)) = ST(ITH,NI(1)) + KELEM1(ITH,IE,IK) * (U(ITH,NI(1)) - UTILDE(ITH)) ! the 2nd term are the theta values of each node ... ST(ITH,NI(2)) = ST(ITH,NI(2)) + KELEM2(ITH,IE,IK) * (U(ITH,NI(2)) - UTILDE(ITH)) ! the 2nd term are the theta values of each node ... ST(ITH,NI(3)) = ST(ITH,NI(3)) + KELEM3(ITH,IE,IK) * (U(ITH,NI(3)) - UTILDE(ITH)) ! the 2nd term are the theta values of each node ... - ENDDO + ENDDO END DO ! IE DO IP = 1, NP DO ITH = 1, NTH @@ -6564,7 +6564,7 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) #ifdef W3_REF1 IF (REFPARS(3).LT.0.5.AND.IOBPD_LOC(ITH,IP).EQ.0.AND.IOBPA_LOC(IP).EQ.0) U(ITH,IP) = UOLD(ITH,IP) ! restores reflected boundary values #endif - ENDDO + ENDDO ENDDO ! IE IF ( FLBPI ) THEN @@ -6586,17 +6586,17 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) U(ITH,JX) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) ) / CGSIG(ISBPI(IBI)) * CLATS(ISBPI(IBI)) END IF END DO - ENDDO + ENDDO ENDIF ! FLBPI CALL PDLIB_exchange2DREAL(U) ENDDO ! IT -! Exact and convert Wave Action + ! Exact and convert Wave Action DO ITH = 1, NTH ISP = ITH + (IK-1) * NTH - DO IP = 1, NP + DO IP = 1, NPA VA(ISP,IP) = U(ITH,IP) * CGSIG(IP) / CLATS(IPLG(IP)) ENDDO ENDDO @@ -6663,9 +6663,9 @@ SUBROUTINE BLOCK_SOLVER_EXPLICIT_INIT() ALLOCATE(FLALL1(NTH,NE,NK), FLALL2(NTH,NE,NK), FLALL3(NTH,NE,NK)) ALLOCATE(KELEM1(NTH,NE,NK), KELEM2(NTH,NE,NK), KELEM3(NTH,NE,NK)) - ALLOCATE(NM(NTH,NE,NK), DTSI(NP)) + ALLOCATE(NM(NTH,NE,NK), DTSI(NP)) ALLOCATE(ITER(NK)) - + !/ ------------------------------------------------------------------- / !/ END SUBROUTINE BLOCK_SOLVER_EXPLICIT_INIT @@ -7346,7 +7346,7 @@ SUBROUTINE DEALLOCATE_PDLIB_GLOBAL(IMOD) INTEGER, INTENT(IN) :: IMOD DEALLOCATE ( & - ! GRIDS(IMOD)%TRIGP, & + ! GRIDS(IMOD)%TRIGP, & GRIDS(IMOD)%SI, & GRIDS(IMOD)%TRIA, & GRIDS(IMOD)%CROSSDIFF, & @@ -7364,8 +7364,8 @@ SUBROUTINE DEALLOCATE_PDLIB_GLOBAL(IMOD) GRIDS(IMOD)%POSI, & GRIDS(IMOD)%I_DIAG, & GRIDS(IMOD)%JA_IE, & - !GRIDS(IMOD)%IOBP, & - !GRIDS(IMOD)%IOBPD, & + !GRIDS(IMOD)%IOBP, & + !GRIDS(IMOD)%IOBPD, & GRIDS(IMOD)%IOBDP, & GRIDS(IMOD)%IOBPA ) !/ @@ -7423,27 +7423,27 @@ SUBROUTINE ERGOUT(FHNDL, ERGNAME) USE W3GDATMD, only: NSPEC, NTH, NK, NSEAL USE W3WDATMD, only: VA, VAOLD IMPLICIT NONE - - INTEGER, INTENT(IN) :: FHNDL + + INTEGER, INTENT(IN) :: FHNDL CHARACTER(LEN=*), INTENT(IN) :: ERGNAME REAL :: SUMVA(NSEAL) - INTEGER :: JSEA + INTEGER :: JSEA IF (LINIT_OUTPUT) THEN OPEN(FHNDL, FILE = TRIM(ERGNAME), FORM = 'UNFORMATTED') - LINIT_OUTPUT = .false. - ENDIF + LINIT_OUTPUT = .false. + ENDIF RTIME = RTIME + 1. - DO JSEA = 1, NSEAL + DO JSEA = 1, NSEAL SUMVA(JSEA) = SUM(VA(:,JSEA)) - ENDDO + ENDDO - WRITE(FHNDL) RTIME + WRITE(FHNDL) RTIME WRITE(FHNDL) (SUMVA(JSEA), SUMVA(JSEA), SUMVA(JSEA), JSEA = 1, NSEAL) - - END SUBROUTINE + + END SUBROUTINE ERGOUT !/ ------------------------------------------------------------------- / SUBROUTINE JACOBI_INIT(IMOD) !/ diff --git a/model/src/w3wavemd.F90 b/model/src/w3wavemd.F90 index fb29f9152..6cbc7e74f 100644 --- a/model/src/w3wavemd.F90 +++ b/model/src/w3wavemd.F90 @@ -1156,19 +1156,19 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & END IF #ifdef W3_DEBUGRUN - DO JSEA = 1, NSEAL - DO IS = 1, NSPEC - IF (VA(IS, JSEA) .LT. 0.) THEN - WRITE(740+IAPROC,*) 'TEST W3WAVE 5', VA(IS,JSEA) - CALL FLUSH(740+IAPROC) - ENDIF + DO JSEA = 1, NSEAL + DO IS = 1, NSPEC + IF (VA(IS, JSEA) .LT. 0.) THEN + WRITE(740+IAPROC,*) 'TEST W3WAVE 5', VA(IS,JSEA) + CALL FLUSH(740+IAPROC) + ENDIF + ENDDO ENDDO - ENDDO - IF (SUM(VA) .NE. SUM(VA)) THEN - WRITE(740+IAPROC,*) 'NAN in ACTION 5', IX, IY, SUM(VA) - CALL FLUSH(740+IAPROC) - STOP - ENDIF + IF (SUM(VA) .NE. SUM(VA)) THEN + WRITE(740+IAPROC,*) 'NAN in ACTION 5', IX, IY, SUM(VA) + CALL FLUSH(740+IAPROC) + STOP + ENDIF #endif call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 6') @@ -1513,48 +1513,48 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & #endif ! #ifdef W3_DEBUGSRC - IF (IX .eq. DEBUG_NODE) THEN - WRITE(740+IAPROC,*) 'NODE_SRCE_IMP_PRE : IX=', IX, ' JSEA=', JSEA - END IF - WRITE(740+IAPROC,*) 'IT/IX/IY/IMOD=', IT, IX, IY, IMOD - WRITE(740+IAPROC,*) 'ISEA/JSEA=', ISEA, JSEA - WRITE(740+IAPROC,*) 'Before sum(VA)=', sum(VA(:,JSEA)) - FLUSH(740+IAPROC) -#endif - CALL W3SRCE(srce_imp_pre, IT, ISEA, JSEA, IX, IY, IMOD, & - VAold(:,JSEA), VA(:,JSEA), & - VSioDummy, VDioDummy, SHAVETOT(JSEA), & - ALPHA(1:NK,JSEA), WN(1:NK,ISEA), & - CG(1:NK,ISEA), CLATS(ISEA), DW(ISEA), U10(ISEA), & - U10D(ISEA), & + IF (IX .eq. DEBUG_NODE) THEN + WRITE(740+IAPROC,*) 'NODE_SRCE_IMP_PRE : IX=', IX, ' JSEA=', JSEA + END IF + WRITE(740+IAPROC,*) 'IT/IX/IY/IMOD=', IT, IX, IY, IMOD + WRITE(740+IAPROC,*) 'ISEA/JSEA=', ISEA, JSEA + WRITE(740+IAPROC,*) 'Before sum(VA)=', sum(VA(:,JSEA)) + FLUSH(740+IAPROC) +#endif + CALL W3SRCE(srce_imp_pre, IT, ISEA, JSEA, IX, IY, IMOD, & + VAold(:,JSEA), VA(:,JSEA), & + VSioDummy, VDioDummy, SHAVETOT(JSEA), & + ALPHA(1:NK,JSEA), WN(1:NK,ISEA), & + CG(1:NK,ISEA), CLATS(ISEA), DW(ISEA), U10(ISEA), & + U10D(ISEA), & #ifdef W3_FLX5 - TAUA(ISEA), TAUADIR(ISEA), & -#endif - AS(ISEA), UST(ISEA), & - USTDIR(ISEA), CX(ISEA), CY(ISEA), & - ICE(ISEA), ICEH(ISEA), ICEF(ISEA), & - ICEDMAX(ISEA), & - REFLEC, REFLED, DELX, DELY, DELA, & - TRNX(IY,IX), TRNY(IY,IX), BERG(ISEA), & - FPIS(ISEA), DTDYN(JSEA), & - FCUT(JSEA), DTGpre, TAUWX(JSEA), TAUWY(JSEA), & - TAUOX(JSEA), TAUOY(JSEA), TAUWIX(JSEA), & - TAUWIY(JSEA), TAUWNX(JSEA), & - TAUWNY(JSEA), PHIAW(JSEA), CHARN(JSEA), & - TWS(JSEA), PHIOC(JSEA), TMP1, D50, PSIC, TMP2, & - PHIBBL(JSEA), TMP3, TMP4, PHICE(JSEA), & - TAUOCX(JSEA), TAUOCY(JSEA), WNMEAN(JSEA), & - RHOAIR(ISEA), ASF(ISEA)) - IF (.not. LSLOC) THEN - VSTOT(:,JSEA) = VSioDummy - VDTOT(:,JSEA) = VDioDummy - ENDIF + TAUA(ISEA), TAUADIR(ISEA), & +#endif + AS(ISEA), UST(ISEA), & + USTDIR(ISEA), CX(ISEA), CY(ISEA), & + ICE(ISEA), ICEH(ISEA), ICEF(ISEA), & + ICEDMAX(ISEA), & + REFLEC, REFLED, DELX, DELY, DELA, & + TRNX(IY,IX), TRNY(IY,IX), BERG(ISEA), & + FPIS(ISEA), DTDYN(JSEA), & + FCUT(JSEA), DTGpre, TAUWX(JSEA), TAUWY(JSEA), & + TAUOX(JSEA), TAUOY(JSEA), TAUWIX(JSEA), & + TAUWIY(JSEA), TAUWNX(JSEA), & + TAUWNY(JSEA), PHIAW(JSEA), CHARN(JSEA), & + TWS(JSEA), PHIOC(JSEA), TMP1, D50, PSIC, TMP2, & + PHIBBL(JSEA), TMP3, TMP4, PHICE(JSEA), & + TAUOCX(JSEA), TAUOCY(JSEA), WNMEAN(JSEA), & + RHOAIR(ISEA), ASF(ISEA)) + IF (.not. LSLOC) THEN + VSTOT(:,JSEA) = VSioDummy + VDTOT(:,JSEA) = VDioDummy + ENDIF #ifdef W3_DEBUGSRC - WRITE(740+IAPROC,*) 'After sum(VA)=', sum(VA(:,JSEA)) - WRITE(740+IAPROC,*) ' sum(VSTOT)=', sum(VSTOT(:,JSEA)) - WRITE(740+IAPROC,*) ' sum(VDTOT)=', sum(VDTOT(:,JSEA)) - WRITE(740+IAPROC,*) ' SHAVETOT=', SHAVETOT(JSEA) - FLUSH(740+IAPROC) + WRITE(740+IAPROC,*) 'After sum(VA)=', sum(VA(:,JSEA)) + WRITE(740+IAPROC,*) ' sum(VSTOT)=', sum(VSTOT(:,JSEA)) + WRITE(740+IAPROC,*) ' sum(VDTOT)=', sum(VDTOT(:,JSEA)) + WRITE(740+IAPROC,*) ' SHAVETOT=', SHAVETOT(JSEA) + FLUSH(740+IAPROC) #endif END DO ! JSEA END IF ! PDLIB From ba748aa43e23a274338637cb5d2bf8d9d3743b41 Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Fri, 7 Apr 2023 12:10:42 -0400 Subject: [PATCH 2/3] re-align comments --- model/src/w3profsmd_pdlib.F90 | 48 +++++++++++++++++------------------ 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index d29cb3f15..787afcf63 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -937,19 +937,19 @@ SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC) USE yowRankModule, only : IPGL_npa INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, - ! actual Wave Direction + ! actual Wave Direction REAL, INTENT(IN) :: DT ! Time intervall for which the - ! advection should be computed - ! for the given velocity field + ! advection should be computed + ! for the given velocity field REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's - ! X- and Y- Components, + ! X- and Y- Components, REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and - ! after advection + ! after advection REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation - ! coefficients for boundary - ! conditions + ! coefficients for boundary + ! conditions LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of - ! the max. Global Time step + ! the max. Global Time step #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif @@ -1265,19 +1265,19 @@ SUBROUTINE PDLIB_W3XYPFSPSI2 ( ISP, C, LCALC, RD10, RD20, DT, AC) USE yowRankModule, only : IPGL_npa IMPLICIT NONE INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, - ! actual Wave Direction + ! actual Wave Direction REAL, INTENT(IN) :: DT ! Time intervall for which the - ! advection should be computed - ! for the given velocity field + ! advection should be computed + ! for the given velocity field REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's - ! X- and Y- Components, + ! X- and Y- Components, REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and - ! after advection + ! after advection REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation - ! coefficients for boundary - ! conditions + ! coefficients for boundary + ! conditions LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of - ! the max. Global Time step + ! the max. Global Time step #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif @@ -1562,19 +1562,19 @@ SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC) IMPLICIT NONE INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber, - ! actual Wave Direction + ! actual Wave Direction REAL, INTENT(IN) :: DT ! Time intervall for which the - ! advection should be computed - ! for the given velocity field + ! advection should be computed + ! for the given velocity field REAL, INTENT(IN) :: C(npa,2) ! Velocity field in it's - ! X- and Y- Components, + ! X- and Y- Components, REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and - ! after advection + ! after advection REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation - ! coefficients for boundary - ! conditions + ! coefficients for boundary + ! conditions LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of - ! the max. Global Time step + ! the max. Global Time step #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif From fa04bec2b33a5a41658d8c1b0aefbc05c0d25b40 Mon Sep 17 00:00:00 2001 From: "denise.worthen" Date: Mon, 10 Apr 2023 08:35:25 -0400 Subject: [PATCH 3/3] switch loop ordering --- model/src/w3profsmd_pdlib.F90 | 37 ++++++++++++++++++----------------- 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index 787afcf63..9564fc49d 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -6435,20 +6435,20 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) CALL WAVNU3 (SIG(IK), DW(iplg(IP)), KSIG(IP), CGSIG(IP)) ENDDO - DO ITH = 1, NTH - DO IP = 1, NPA + DO IP = 1, NPA + DO ITH = 1, NTH ISEA = IPLG(IP) CXX(ITH,IP) = CGSIG(IP) * FACX * ECOS(ITH) / CLATS(ISEA) CYY(ITH,IP) = CGSIG(IP) * FACY * ESIN(ITH) - ENDDO + ENDDO ! ith IF (FLCUR) THEN - DO IP = 1, NPA + DO ITH = 1, NTH ISEA = IPLG(IP) IF (IOBP_LOC(IP) .GT. 0) THEN CXX(ITH,IP) = CXX(ITH,IP) + FACX * CX(ISEA)/CLATS(ISEA) CYY(ITH,IP) = CYY(ITH,IP) + FACY * CY(ISEA) ENDIF - ENDDO + ENDDO !ith ENDIF ENDDO @@ -6538,12 +6538,13 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) END IF ! LCALC ! Exact and convert Wave Action - should be some subroutine function or whatever - DO ITH = 1, NTH - ISP = ITH + (IK-1) * NTH - DO IP = 1, NPA - U(ITH,IP) = VA(ISP,IP) / CGSIG(IP) * CLATS(IPLG(IP)) - ENDDO - ENDDO + do ip = 1,npa + isp = 0 + do ith = 1,nth + isp = ith + (ik-1)*nth + u(ith,ip) = va(isp,ip) / cgsig(ip) * clats(iplg(ip)) + enddo + enddo UOLD = U DO IT = 1, ITER(IK) @@ -6559,7 +6560,6 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) END DO ! IE DO IP = 1, NP DO ITH = 1, NTH - ISP = ITH + (IK-1) * NTH U(ITH,IP) = MAX(ZERO,U(ITH,IP)-DTSI(IP)*ST(ITH,IP)*(1-IOBPA_LOC(IP)))*IOBPD_LOC(ITH,IP)*IOBDP_LOC(IP) #ifdef W3_REF1 IF (REFPARS(3).LT.0.5.AND.IOBPD_LOC(ITH,IP).EQ.0.AND.IOBPA_LOC(IP).EQ.0) U(ITH,IP) = UOLD(ITH,IP) ! restores reflected boundary values @@ -6594,12 +6594,13 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) ENDDO ! IT ! Exact and convert Wave Action - DO ITH = 1, NTH - ISP = ITH + (IK-1) * NTH - DO IP = 1, NPA - VA(ISP,IP) = U(ITH,IP) * CGSIG(IP) / CLATS(IPLG(IP)) - ENDDO - ENDDO + do ip = 1,npa + isp = 0 + do ith = 1,nth + isp = ith + (ik-1)*nth + va(isp,ip) = u(ith,ip) * cgsig(ip) / clats(iplg(ip)) + end do + end do ENDDO ! IK