#include "w3macros.h"
!/
!/ ------------------------------------------------------------------- /
MODULE PDLIB_W3PROFSMD
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 95 |
!/ | Last update : 1-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2016 : Origination ( version 6.04 )
!/
! 1. Purpose : PDLIB version of UGTYPE including fully implicit
! discretization. This works is based on the thesis
! of Roland, 2008 and represents the continues
! development of the solution of the WAE on unstructured
! grids. Following the quest since one decade we
! continuesly improve the aplicability and robustness of
! the source code and the methods. The development and
! implementation of the involved schemes was funded over
! the past decade by IFREMER, SHOM, USACE, NCEP/NOAA,
! BGS IT&E GmbH, Zanke & Partner and Roland & Partner.
! The PDLIB (Parallel Decomposition Library) library,
! which is used here is courtesy to BGS IT&E GmbH and
! has it's own license, which is the same as WW3. As of
! the origin of the methods, ideas and source code. This
! code was 1st developed in the WWM-III (Roland, 2008) and
! the ported to WW3. This is true for all source code
! related to UGTYPE.
!
!
! 2. Method : We apply here the framework of Residual Distributions
! schemes for hyperbolic problems for nonlinear propagation
! laws based on the work of Richiuotto et al. 2005.
! We supply the N-scheme, PSI-scheme and Lax-FCT-scheme
! as explicit methods ranging from 1st order time space
! to most optimal PSI method up to 2nd order Lax-FCT-scheme.
! For the implicit implementation we used up to now ONLY
! the N-Scheme. Higher order schemes are up to now rather
! a research feature than for practical application. The
! reason is given in Cavalleri et al. 2018, we do not
! resolve enough physics in order to be able to run
! 2nd or even higher order schemes.
! Use the numerical schemes with the needed care and
! do proper convergence analysis on the time step and
! grid size as well as solver threshold depedency.
! Think about the time and spatial scales of your
! u r intending to resolve! Multiscale modelling needs
! much work on the side of the modeler and much more
! time for the grid generation and the validation of the
! final model
!
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/ ------------------------------------------------------------------- /
!/
PUBLIC
!/
!/ Public variables
!/
LOGICAL :: MAPSTA_HACK = .FALSE.
REAL, ALLOCATABLE :: ASPAR_JAC(:,:), ASPAR_DIAG_SOURCES(:,:)
REAL, ALLOCATABLE :: B_JAC(:,:), CAD_THE(:,:), CAS_SIG(:,:)
REAL, ALLOCATABLE :: CWNB_SIG_M2(:,:)
REAL, ALLOCATABLE :: U_JAC(:,:)
REAL, ALLOCATABLE :: COFRM4(:)
INTEGER, ALLOCATABLE :: IS0_pdlib(:)
INTEGER :: FreqShiftMethod = 2
LOGICAL :: FSGEOADVECT
!/DEBUGSRC INTEGER :: TESTNODE = 1
!
!/ ------------------------------------------------------------------- /
!
CONTAINS
!
!/ ------------------------------------------------------------------- /
!
SUBROUTINE VA_SETUP_IOBPD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-Mai-2018 |
!/ +-----------------------------------+
!/
!/ 01-Mai-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Setup boundary pointer
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: IOBPD, GTYPE, UNGTYPE
USE W3GDATMD, ONLY: NSPEC, NTH, NSEAL
USE W3WDATMD, ONLY: VA
USE YOWNODEPOOL, ONLY: iplg
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!
INTEGER JSEA, IP, IP_glob, ITH, ISP
!/S CALL STRACE (IENT, 'VA_SETUP_IOBPD')
IF (GTYPE .eq. UNGTYPE) THEN
DO JSEA=1,NSEAL
IP = JSEA
IP_glob = iplg(IP)
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
VA(ISP,JSEA) = VA(ISP,JSEA)*IOBPD(ITH,IP_glob)
END DO
END DO
END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_STYLE_INIT(IMOD)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Init pdlib part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: FLCX, FLCY
USE CONSTANTS, ONLY : GRAV, TPI
USE W3GDATMD, ONLY: XYB, XGRD, YGRD, NX, NSEA, NTRI, TRIGP, NSPEC
USE W3GDATMD, ONLY: MAPSTA, MAPFS, GRIDS, NTH
USE W3GDATMD, ONLY: IOBP, IOBPD, IOBP_loc, IOBPD_loc, SIG, NK
USE W3GDATMD, ONLY: TRIA, IEN, LEN, ANGLE, ANGLE0
USE W3GDATMD, ONLY: CCON, COUNTCON, INDEX_CELL, IE_CELL
USE W3GDATMD, ONLY: POS_CELL, SI, IAA, JAA, POSI, I_DIAG, JA_IE
USE W3ADATMD, ONLY: MPI_COMM_WCMP, MPI_COMM_WAVE
USE W3ODATMD, ONLY: IAPROC, NAPROC, NTPROC
USE yowDatapool, ONLY: istatus
USE yowpdlibMain, ONLY: initFromGridDim
USE YOWNODEPOOL, ONLY: npa, iplg
USE W3PARALL, ONLY : PDLIB_NSEAL, PDLIB_NSEALM
USE W3PARALL, ONLY : JX_TO_JSEA, ISEA_TO_JSEA
USE yowfunction, ONLY : ComputeListNP_ListNPA_ListIPLG, pdlib_abort
!/
IMPLICIT NONE
INCLUDE "mpif.h"
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!! INCLUDE "mpif.h"
REAL, ALLOCATABLE :: XP_IN(:), YP_IN(:), DEP_IN(:)
INTEGER, ALLOCATABLE :: INE_IN(:,:)
INTEGER :: istat
INTEGER :: I, J, IBND_MAP, ISEA, IP, IX, JSEA, nb
INTEGER :: IP_glob
INTEGER :: myrank, ierr, iproc
INTEGER, ALLOCATABLE :: NSEAL_arr(:)
INTEGER :: IERR_MPI
INTEGER :: IScal(1)
INTEGER, INTENT(in) :: IMOD
INTEGER :: IK, ISP
INTEGER IK0, ISP0, ITH
REAL :: eSIG, eFR
REAL, PARAMETER :: COEF4 = 5.0E-7
!/S CALL STRACE (IENT, 'PDLIB_STYLE_INIT')
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_STYLE_INIT, IMOD (no print)'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NAPROC=', NAPROC
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NTPROC=', NTPROC
!/DEBUGSOLVER FLUSH(740+IAPROC)
PDLIB_NSEAL=0
IF (IAPROC .le. NAPROC) THEN
ALLOCATE(XP_IN(NX), YP_IN(NX), DEP_IN(NX), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(1)
DO I=1,NX
XP_IN(I)=XYB(I,1)
YP_IN(I)=XYB(I,2)
DEP_IN(I)=XYB(I,3)
END DO
ALLOCATE(INE_IN(3,NTRI), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(2)
DO I=1,NTRI
DO J=1,3
INE_IN(J,I)=TRIGP(I,J)
END DO
END DO
CALL MPI_COMM_RANK(MPI_COMM_WCMP, myrank, ierr)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_STYLE_INIT, IAPROC=', IAPROC
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_STYLE_INIT, NAPROC=', NAPROC
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_STYLE_INIT, myrank=', myrank
!/DEBUGSOLVER FLUSH(740+IAPROC)
!
CALL initFromGridDim(NX,XP_IN,YP_IN,DEP_IN,NTRI,INE_IN,NSPEC,MPI_COMM_WCMP)
!
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'After initFromGridDim'
!/DEBUGSOLVER FLUSH(740+IAPROC)
DEALLOCATE(XP_IN, YP_IN, DEP_IN, INE_IN)
!
! Now the computation of NSEAL
!
DO IP = 1, npa
IX = iplg(IP)
ISEA = MAPFS(1,IX)
IF (ISEA .gt. 0) PDLIB_NSEAL = PDLIB_NSEAL + 1
END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'npa is augmented domain over NX'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_NSEAL is basicall npa but ONLY over the wet points'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NSEAL is set to PDLIB_NSEAL'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_NSEAL=', PDLIB_NSEAL
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'npa=', npa
!/DEBUGSOLVER FLUSH(740+IAPROC)
ALLOCATE(JX_TO_JSEA(npa), ISEA_TO_JSEA(NSEA), stat=istat)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISEA_TO_JSEA ALLOCATEd'
!/DEBUGSOLVER FLUSH(740+IAPROC)
if(istat /= 0) CALL PDLIB_ABORT(3)
JSEA = 0
JX_TO_JSEA = 0
ISEA_TO_JSEA = 0
DO IP = 1, npa
IX = iplg(IP)
ISEA = MAPFS(1,IX)
IF (ISEA .gt. 0) THEN
JSEA=JSEA+1
JX_TO_JSEA(IP)=JSEA
ISEA_TO_JSEA(ISEA)=JSEA
END IF
END DO
!
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'After JX_TO_JSEA, ISEA_TO_JSEA and friend computation'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!
! Map a point in (1:PDLIB_NSEAL) to a point in (1:NSEA)
!
nb=0
DO IX=1,NX
IF (MAPFS(1,IX) .gt. 0) nb = nb + 1
END DO
IF (nb .ne. NSEA) THEN
WRITE(*,*) 'Logical error in computation of NSEA / nb'
WRITE(*,*) 'nb=', nb, ' NSEA=', NSEA
STOP
END IF
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'nb / NSEA consistency check'
!/DEBUGSOLVER FLUSH(740+IAPROC)
END IF
FSGEOADVECT=.FALSE.
IF ((FLCX .eqv. .TRUE.).and.(FLCY .eqv. .TRUE.)) THEN
FSGEOADVECT=.TRUE.
END IF
!
! Compute NSEALM
!
IF (IAPROC .le. NAPROC) THEN
IF (IAPROC .eq. 1) THEN
ALLOCATE(NSEAL_arr(NAPROC))
NSEAL_arr(1)=PDLIB_NSEAL
DO IPROC=2,NAPROC
CALL MPI_RECV(IScal,1,MPI_INT, IPROC-1, 23, MPI_COMM_WAVE, istatus, IERR_MPI)
NSEAL_arr(IPROC)=IScal(1)
END DO
PDLIB_NSEALM=maxval(NSEAL_arr)
deALLOCATE(NSEAL_arr)
ELSE
IScal(1)=PDLIB_NSEAL
CALL MPI_SEND(IScal,1,MPI_INT, 0, 23, MPI_COMM_WAVE, IERR_MPI)
END IF
END IF
!
IF (IAPROC .eq. 1) THEN
IScal(1)=PDLIB_NSEALM
DO IPROC = 2 , NTPROC
CALL MPI_SEND(IScal,1,MPI_INT, IPROC-1, 24, MPI_COMM_WAVE, IERR_MPI)
END DO
ELSE
CALL MPI_RECV(IScal,1,MPI_INT, 0, 24, MPI_COMM_WAVE, istatus, IERR_MPI)
PDLIB_NSEALM=IScal(1)
END IF
!
WRITE(*,'(A40,3I20)') 'IAPROC, PDLIB_NSEAL, PDLIB_NSEALM', IAPROC, PDLIB_NSEAL, PDLIB_NSEALM
!
! Computation of IOBP_loc, IOBPD_loc
!
ALLOCATE(GRIDS(IMOD)%IOBP_loc(npa), GRIDS(IMOD)%IOBPD_loc(NTH,npa), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(4)
DO IP=1,npa
IP_glob=iplg(IP)
GRIDS(IMOD)%IOBP_loc(IP) = IOBP(IP_glob)
GRIDS(IMOD)%IOBPD_loc(:,IP) = IOBPD(:,IP_glob)
END DO
!/DEBUGINIT WRITE(740+IAPROC,*) 'ALLOCATEd(ISEA_TO_JSEA)=', allocated(ISEA_TO_JSEA)
!/DEBUGINIT WRITE(740+IAPROC,*) 'PDLIB_NSEALM=', PDLIB_NSEALM
!/DEBUGINIT FLUSH(740+IAPROC)
! WRITE(*,*) 'Begin, ComputeListNP_ListNPA_ListIPLG'
CALL ComputeListNP_ListNPA_ListIPLG
ALLOCATE(COFRM4(NK))
DO IK=1,NK
eSIG=SIG(IK)
eFR=eSIG/TPI
COFRM4(IK)=COEF4*GRAV/(eFR**4)
END DO
ALLOCATE(IS0_pdlib(NSPEC))
DO ISP=1, NSPEC
IS0_pdlib(ISP) = ISP - 1
END DO
DO ISP=1, NSPEC, NTH
IS0_pdlib(ISP) = IS0_pdlib(ISP) + NTH
END DO
!
! DEALLOCATE USELESS STUFF
!
DEALLOCATE(LEN, ANGLE, ANGLE0, TRIA)
DEALLOCATE(CCON, COUNTCON, INDEX_CELL, IE_CELL)
DEALLOCATE(POS_CELL, SI, IAA, JAA, POSI, I_DIAG, JA_IE)
!/
!/ End of PDLIB_STYLE_INIT ------------------------------------------- /
!/
END SUBROUTINE PDLIB_STYLE_INIT
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_MAPSTA_INIT(IMOD)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Init mapsta part for pdlib
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : INDEX_MAP, NBND_MAP, NSEA, NSEAL, MAPSTA, GRIDS, NX
USE W3ODATMD, ONLY : IAPROC, NAPROC
USE YOWNODEPOOL, ONLY: iplg, npa
USE yowfunction, ONLY: pdlib_abort
USE W3ODATMD, ONLY: IAPROC
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER :: IBND_MAP, ISEA, JSEA, IX, IP, IP_glob
INTEGER, INTENT(in) :: IMOD
INTEGER :: Status(NX), istat
!/S CALL STRACE (IENT, 'PDLIB_MAPSTA_INIT')
!/DEBUGINIT WRITE(*,*) 'Passing by PDLIB_MAPSTA_INIT IAPROC=', IAPROC
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
ALLOCATE(GRIDS(IMOD)%MAPSTA_LOC(npa), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(5)
Status=0
DO IP=1,npa
IP_glob=iplg(IP)
Status(IP_glob)=IP
GRIDS(IMOD)%MAPSTA_LOC(IP)=MAPSTA(1,IP_glob)
END DO
NBND_MAP=0
DO IX=1,NX
IF ((MAPSTA(1,IX) .lt. 1).and.(Status(IX).gt.0)) THEN
NBND_MAP = NBND_MAP + 1
END IF
END DO
ALLOCATE(GRIDS(IMOD)%INDEX_MAP(NBND_MAP), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(6)
IBND_MAP=0
DO IX=1,NX
IF ((MAPSTA(1,IX) .lt. 1).and.(Status(IX).gt.0)) THEN
IBND_MAP = IBND_MAP + 1
GRIDS(IMOD)%INDEX_MAP(IBND_MAP) = Status(IX)
END IF
END DO
!/
!/ End of W3SPR4 ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 10-Jan-2011 |
!/ +-----------------------------------+
!/
!/ 10-Jan-2008 : Origination. ( version 3.13 )
!/ 10-Jan-2011 : Addition of implicit scheme ( version 3.14.4 )
!/ 06-Feb-2014 : PDLIB parallelization
!/
! 1. Purpose : Explicit advection schemes driver
!
! Propagation in physical space for a given spectral component.
! Gives the choice of scheme on unstructured grid
! Use the geographical parall algorithms for further speed.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! Local variables.
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! 5. Called by :
!
! W3WAVE Wave model routine.
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
! make the interface between the WAVEWATCH and the WWM code.
!
! 8. Structure :
!
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
!/
!
USE CONSTANTS
!
USE W3TIMEMD, ONLY: DSEC21
!
USE W3GDATMD, ONLY: NX, NY, MAPFS, CLATS, &
FLCX, FLCY, NK, NTH, DTH, XFR, &
ECOS, ESIN, SIG, PFMOVE, &
IOBP, IOBPD, &
FSN, FSPSI, FSFCT, FSNIMP, &
GTYPE, UNGTYPE, NBND_MAP, INDEX_MAP
USE YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA
USE YOWNODEPOOL, ONLY: iplg, npa
USE W3WDATMD, ONLY: TIME, VA
USE W3ODATMD, ONLY: TBPI0, TBPIN, FLBPI
USE W3ADATMD, ONLY: CG, CX, CY, ITIME
USE W3IDATMD, ONLY: FLCUR
USE W3GDATMD, ONLY: NSEAL
USE W3ODATMD, ONLY: IAPROC
IMPLICIT NONE
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: ISP
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
LOGICAL, INTENT(IN) :: LCALC
LOGICAL :: SCHEME
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
INTEGER :: ITH, IK, ISEA
INTEGER :: I, J, IE, IBND_MAP
INTEGER :: IP_glob
REAL :: CCOS, CSIN, CCURX, CCURY
REAL :: C(npa,2)
REAL :: RD1, RD2
!/
!/ Automatic work arrays
!/
REAL :: VLCFLX(npa), VLCFLY(npa)
REAL :: AC(npa)
REAL :: AC_MAP(NBND_MAP)
INTEGER :: JSEA, IP
!/ ------------------------------------------------------------------- /
!
! 1. Preparations --------------------------------------------------- *
! 1.a Set constants
!
!/S CALL STRACE (IENT, 'W3XYPUG')
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Begin of PDLIB_W3XYPUG'
!/DEBUGSOLVER FLUSH(740+IAPROC)
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
CCURX = FACX
CCURY = FACY
!
! 1.b Initialize arrays
!
VLCFLX = 0.
VLCFLY = 0.
AC = 0.
CALL SETDEPTH_PDLIB
!
! 2. Calculate velocities ---------------- *
!
DO JSEA = 1, NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA = MAPFS(1,IP_glob)
!write(*,*) 'IP TEST', JSEA, ISEA, IP, IP_glob
AC(IP) = VA(ISP,JSEA) / CG(IK,ISEA) * CLATS(ISEA)
VLCFLX(IP) = CCOS * CG(IK,ISEA) / CLATS(ISEA)
VLCFLY(IP) = CSIN * CG(IK,ISEA)
!/MGP VLCFLX(IP) = VLCFLX(IP) - CCURX*VGX/CLATS(ISEA)
!/MGP VLCFLY(IP) = VLCFLY(IP) - CCURY*VGY
END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, ' ITH=', ITH, ' IK=', IK
!/DEBUGSOLVER WRITE(740+IAPROC,*) '1: maxval(VLCFLX)=', maxval(VLCFLX)
!/DEBUGSOLVER WRITE(740+IAPROC,*) '1: maxval(VLCFLY)=', maxval(VLCFLY)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FLCUR=', FLCUR
!/DEBUGSOLVER FLUSH(740+IAPROC)
IF ( FLCUR ) THEN
DO JSEA=1, NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA = MAPFS(1,IP_glob)
!
! Currents are not included on coastal boundaries (COUNTSEACON(IXY) .NE. PDLIB_CCON(IXY))
!
IF (IOBP(IP_glob) .GT. 0) THEN
VLCFLX(IP) = VLCFLX(IP) + CCURX*CX(ISEA)/CLATS(ISEA)
VLCFLY(IP) = VLCFLY(IP) + CCURY*CY(ISEA)
END IF
END DO
END IF
C(:,1) = VLCFLX(:)
C(:,2) = VLCFLY(:)
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'CCURXY=', CCURX, CCURY
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'max(CX)=', maxval(CX)
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'max(CY)=', maxval(CY)
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'min(CLATS)=', minval(CLATS)
!!/DEBUGSOLVER WRITE(740+IAPROC,*) '2: maxval(VLCFLX)=', maxval(VLCFLX)
!!/DEBUGSOLVER WRITE(740+IAPROC,*) '2: maxval(VLCFLY)=', maxval(VLCFLY)
!!/DEBUGSOLVER FLUSH(740+IAPROC)
!
! 4. Prepares boundary update
!
IF ( FLBPI ) THEN
RD1 = DSEC21 ( TBPI0, TIME )
RD2 = DSEC21 ( TBPI0, TBPIN )
ELSE
RD1=1.
RD2=0.
END IF
!
! Saving data for MAPSTA business
!
IF (MAPSTA_HACK) THEN
DO IBND_MAP=1,NBND_MAP
IP=INDEX_MAP(IBND_MAP)
AC_MAP(IBND_MAP) = AC(IP)
END DO
END IF
!
! 4. propagate using the selected scheme
!
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'maxval(C)=', maxval(C)
!/DEBUGSOLVER FLUSH(740+IAPROC)
IF (FSN) THEN
CALL PDLIB_W3XYPFSN2(ISP, C, LCALC, RD1, RD2, DTG, AC)
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)
ELSE IF (FSNIMP) THEN
STOP 'For PDLIB and FSNIMP, no function has been programmed yet'
ENDIF
!
IF (MAPSTA_HACK) THEN
DO IBND_MAP=1,NBND_MAP
IP=INDEX_MAP(IBND_MAP)
AC(IP) = AC_MAP(IBND_MAP)
END DO
END IF
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'After solutioning'
!/DEBUGSOLVER FLUSH(740+IAPROC)
! 6. Store results in VQ in proper format --------------------------- *
!
!/C90/!DIR$ IVDEP
!
DO JSEA=1, NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA=MAPFS(1,IP_glob)
VA(ISP,JSEA) = MAX ( 0. , CG(IK,ISEA)/CLATS(ISEA)*AC(IP) )
END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Leaving PDLIB_W3XYPUG'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/
!/ End of W3SPR4 ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Explicit N-Scheme
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NK, NTH, NX, &
IEN, CLATS, MAPSF, IOBPD, IOBP, IOBDP, IOBPA
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: CG, ITER, DW , CFLXYMAX, NSEALM
USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPIN, ISBPI, BBPI0, BBPIN
USE W3TIMEMD, ONLY: DSEC21
USE W3ADATMD, ONLY: MPI_COMM_WCMP
USE W3GDATMD, ONLY: NSEAL, DMIN, NSEA
!/REF1 USE W3GDATMD, ONLY: REFPARS
USE YOWNODEPOOL, ONLY: PDLIB_SI, PDLIB_IEN, PDLIB_TRIA, ipgl, iplg, npa, np
use yowElementpool, ONLY: ne, INE
use yowDatapool, ONLY: rtype
use yowExchangeModule, ONLY : PDLIB_exchange1DREAL
USE W3ODATMD, ONLY : IAPROC
USE MPI, ONLY : MPI_MIN
USE W3PARALL, ONLY : INIT_GET_JSEA_ISPROC
USE W3PARALL, ONLY : ONESIXTH, ZERO, THR
USE yowRankModule, ONLY : IPGL_npa
IMPLICIT NONE
INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber,
! actual Wave Direction
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
! X- and Y- Components,
REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and
! after advection
REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation
! coefficients for boundary
! conditions
LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of
! the max. Global Time step
!/S INTEGER, SAVE :: IENT = 0
!/REF1 INTEGER(KIND=1) :: IOBPDR(NX)
INTEGER :: IP, IE, POS, IT, I1, I2, I3, I, J, ITH, IK
INTEGER :: IBI, NI(3)
INTEGER :: JX
!
! local REAL
!
REAL :: RD1, RD2
!:
! local double
!
REAL :: UTILDE
REAL :: SUMTHETA
REAL :: FT, CFLXY
REAL :: FL11, FL12, FL21, FL22, FL31, FL32
REAL :: FL111, FL112, FL211, FL212, FL311, FL312
REAL :: DTSI(npa), U(npa)
REAL :: DTMAX_GL, DTMAX, DTMAXEXP, REST
REAL :: LAMBDA(2), KTMP(3)
REAL :: KELEM(3,NE), FLALL(3,NE)
REAL :: KKSUM(npa), ST(npa)
REAL :: NM(NE)
INTEGER :: ISPROC, JSEA, IP_glob, ierr, IX
REAL :: eSumAC, sumAC, sumBPI0, sumBPIN, sumCG, sumCLATS
LOGICAL :: testWrite
REAL :: FIN(1), FOUT(1)
!/S CALL STRACE (IENT, 'W3XYPFSN')
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 1'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVER CALL GET_SCAL_INTEGRAL_COH_R4(AC, eSumAC)
!/DEBUGSOLVER testWrite=.FALSE.
!/DEBUGSOLVER IF (eSumAC .gt. 0) THEN
!/DEBUGSOLVER testWrite=.TRUE.
!/DEBUGSOLVER END IF
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R4(AC, "AC in input")
!/DEBUGSOLVER END IF
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
DTMAX = DBLE(10.E10)
! DO IP = 1, npa
! IP_glob = iplg(IP)
! IF (IOBPD(ITH,IP_glob) .EQ. 0) THEN
! if ( AC(IP) .gt. 0.) write(*,*) 'TEST w3profs_pdlib 0', ip, ith, ik, AC(IP)
! endif
! END DO
!
!/REF1 IOBPDR(:)=(1-IOBP(:))*(1-IOBPD(ITH,:))
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NX=', NX
!/DEBUGSOLVER DO IX=1,NX
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'IX/IOBP=', IX, IOBP(IX)
!/DEBUGSOLVER END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'sum(IOBP)=', sum(IOBP)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 2'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!
!2 Propagation
!2.a Calculate K-Values and contour based quantities ...
!
DO IE = 1, ne
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
LAMBDA(1) = ONESIXTH *(C(I1,1)+C(I2,1)+C(I3,1)) ! Linearized advection speed in X and Y direction
LAMBDA(2) = ONESIXTH *(C(I1,2)+C(I2,2)+C(I3,2))
KELEM(1,IE) = LAMBDA(1) * PDLIB_IEN(1,IE) + LAMBDA(2) * PDLIB_IEN(2,IE) ! K-Values - so called Flux Jacobians
KELEM(2,IE) = LAMBDA(1) * PDLIB_IEN(3,IE) + LAMBDA(2) * PDLIB_IEN(4,IE)
KELEM(3,IE) = LAMBDA(1) * PDLIB_IEN(5,IE) + LAMBDA(2) * PDLIB_IEN(6,IE)
KTMP = KELEM(:,IE) ! Copy
NM(IE) = - 1.D0/MIN(-THR,SUM(MIN(ZERO,KTMP))) ! N-Values
KELEM(:,IE) = MAX(ZERO,KTMP)
FL11 = C(I2,1) * PDLIB_IEN(1,IE) + C(I2,2) * PDLIB_IEN(2,IE) ! Weights for Simpson Integration
FL12 = C(I3,1) * PDLIB_IEN(1,IE) + C(I3,2) * PDLIB_IEN(2,IE)
FL21 = C(I3,1) * PDLIB_IEN(3,IE) + C(I3,2) * PDLIB_IEN(4,IE)
FL22 = C(I1,1) * PDLIB_IEN(3,IE) + C(I1,2) * PDLIB_IEN(4,IE)
FL31 = C(I1,1) * PDLIB_IEN(5,IE) + C(I1,2) * PDLIB_IEN(6,IE)
FL32 = C(I2,1) * PDLIB_IEN(5,IE) + C(I2,2) * PDLIB_IEN(6,IE)
FL111 = 2.d0*FL11+FL12
FL112 = 2.d0*FL12+FL11
FL211 = 2.d0*FL21+FL22
FL212 = 2.d0*FL22+FL21
FL311 = 2.d0*FL31+FL32
FL312 = 2.d0*FL32+FL31
FLALL(1,IE) = (FL311 + FL212) * ONESIXTH + KELEM(1,IE)
FLALL(2,IE) = (FL111 + FL312) * ONESIXTH + KELEM(2,IE)
FLALL(3,IE) = (FL211 + FL112) * ONESIXTH + KELEM(3,IE)
END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 3'
!/DEBUGSOLVER FLUSH(740+IAPROC)
IF (LCALC) THEN
KKSUM = ZERO
DO IE = 1, NE
NI = INE(:,IE)
KKSUM(NI) = KKSUM(NI) + KELEM(:,IE)
END DO
DO IP = 1, npa
IP_glob=iplg(IP)
DTMAXEXP = PDLIB_SI(IP)/MAX(DBLE(10.E-10),KKSUM(IP)*IOBDP(IP))
DTMAX = MIN( DTMAX, DTMAXEXP)
CFLXYMAX(IP) = MAX(CFLXYMAX(IP),DBLE(DT)/DTMAXEXP)
END DO
FIN(1)=DTMAX
CALL MPI_ALLREDUCE(FIN,FOUT,1,rtype,MPI_MIN,MPI_COMM_WCMP,ierr)
DTMAX_GL=FOUT(1)
CFLXY = DBLE(DT)/DTMAX_GL
REST = ABS(MOD(CFLXY,1.0d0))
IF (REST .LT. THR) THEN
ITER(IK,ITH) = ABS(NINT(CFLXY))
ELSE IF (REST .GT. THR .AND. REST .LT. 0.5d0) THEN
ITER(IK,ITH) = ABS(NINT(CFLXY)) + 1
ELSE
ITER(IK,ITH) = ABS(NINT(CFLXY))
END IF
END IF ! LCALC
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 4'
!/DEBUGSOLVER FLUSH(740+IAPROC)
DO IP = 1, npa
DTSI(IP) = DBLE(DT)/DBLE(ITER(IK,ITH))/PDLIB_SI(IP) ! Some precalculations for the time integration.
END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 4.1'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R8(PDLIB_SI, "PDLIB_SI in input")
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 4.2'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R8(DTSI, "DTSI in input")
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 5'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'IK=', IK, ' ITH=', ITH
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ITER=', ITER(IK,ITH)
!/DEBUGSOLVER FLUSH(740+IAPROC)
DO IT = 1, ITER(IK,ITH)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'IK=', IK, ' ITH=', ITH
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'IT=', IT, ' ITER=', ITER(IK,ITH)
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'IT=', IT
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVER END IF
U = DBLE(AC)
ST = ZERO
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 ...
END DO ! IE
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R8(ST, "ST in loop")
!/DEBUGSOLVER END IF
!
! IOBPD=0 : waves coming from land
! IOBPD=1 : waves coming from the coast
!
DO IP = 1, npa
IP_glob=iplg(IP)
U(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*(1-IOBPA(IP_glob)))*DBLE(IOBPD(ITH,IP_glob))*IOBDP(IP_glob)
!/REF1 IF (REFPARS(3).LT.0.5.AND.IOBPD(ITH,IP_glob).EQ.0.AND.IOBPA(IP_glob).EQ.0) U(IP) = AC(IP) ! restores reflected boundary values
END DO
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R8(U, "U in loop")
!/DEBUGSOLVER END IF
AC = REAL(U)
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R4(AC, "AC before synchronization")
!/DEBUGSOLVER END IF
CALL PDLIB_exchange1DREAL(AC)
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R4(AC, "AC after synchronization")
!/DEBUGSOLVER END IF
!
! 5 Update boundaries ... would be better to omit any if clause in this loop ...
! a possibility would be to use NBI = 0 when FLBPI is FALSE and loop on IBI whatever the value of NBI
!
IF ( FLBPI ) THEN
RD1=RD10 - DT * REAL(ITER(IK,ITH)-IT)/REAL(ITER(IK,ITH))
RD2=RD20
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
ELSE
RD1 = 0.
RD2 = 1.
END IF
!/DEBUGSOLVER sumAC=0
!/DEBUGSOLVER sumBPI0=0
!/DEBUGSOLVER sumBPIN=0
!/DEBUGSOLVER sumCG=0
!/DEBUGSOLVER sumCLATS=0
DO IBI = 1, NBI
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) ) &
/ CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
!/DEBUGSOLVER sumAC=sumAC + AC(JX)
!/DEBUGSOLVER sumBPI0=sumBPI0 + BBPI0(ISP,IBI)
!/DEBUGSOLVER sumBPIN=sumBPIN + BBPIN(ISP,IBI)
!/DEBUGSOLVER sumCG=sumCG + CG(IK,ISBPI(IBI))
!/DEBUGSOLVER sumCLATS=sumCLATS + CLATS(ISBPI(IBI))
END IF
END DO
END IF
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NBI=', NBI
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'RD1=', RD1, ' RD2=', RD2
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumAC=', sumAC
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumBPI0=', sumBPI0
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumBPIN=', sumBPIN
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumCG=', sumCG
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumCLATS=', sumCLATS
!/DEBUGSOLVER FLUSH(740+IAPROC)
CALL PDLIB_exchange1DREAL(AC)
!/DEBUGSOLVER IF (testWrite) THEN
!/DEBUGSOLVER CALL SCAL_INTEGRAL_PRINT_R4(AC, "AC after FLBPI")
!/DEBUGSOLVER END IF
END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_W3XYPFSN2, step 6'
!/DEBUGSOLVER FLUSH(740+IAPROC)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPFSPSI2 ( ISP, C, LCALC, RD10, RD20, DT, AC)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Explicit PSI-Scheme
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NK, NTH, NX, &
IEN, CLATS, MAPSF, IOBPD, IOBP, IOBPA, NNZ, IOBDP
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: CG, ITER, CFLXYMAX
USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, ISBPI, BBPI0, BBPIN
USE W3TIMEMD, ONLY: DSEC21
USE W3GDATMD, ONLY: NSEAL, DMIN
!/REF1 USE W3GDATMD, ONLY: REFPARS
USE W3ADATMD, ONLY: MPI_COMM_WCMP
use yowElementpool, ONLY: ne, INE
use YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA, PDLIB_SI, iplg, npa
USE W3ODATMD, ONLY : IAPROC
use yowDatapool, ONLY: rtype
use yowExchangeModule, ONLY : PDLIB_exchange1DREAL
USE MPI, ONLY : MPI_MIN
USE W3PARALL, ONLY : INIT_GET_JSEA_ISPROC
USE W3PARALL, ONLY : ONESIXTH, THR, ZERO
USE yowRankModule, ONLY : IPGL_npa
IMPLICIT NONE
INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber,
! actual Wave Direction
REAL, INTENT(IN) :: DT ! Time interval for which the
! advection should be computed
! for the given velocity field
REAL, INTENT(IN) :: C(npa,2) ! Velocity field in its
! X- and Y- Components,
REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and
! after advection
REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation
! coefficients for boundary
! conditions
LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of
! the max. Global Time step
!/REF1 INTEGER(KIND=1) :: IOBPDR(NX)
!/S INTEGER, SAVE :: IENT = 0
INTEGER :: IP, IE, POS, IT, I1, I2, I3, I, J, ITH, IK
INTEGER :: IBI, NI(3), JX
INTEGER :: ISPROC, IP_glob, JSEA, ierr
REAL :: RD1, RD2
REAL :: UTILDE
REAL :: SUMTHETA
REAL :: FL1, FL2, FL3
REAL :: FT, CFLXY
REAL :: FL11, FL12, FL21, FL22, FL31, FL32
REAL :: FL111, FL112, FL211, FL212, FL311, FL312
REAL :: DTSI(npa), U(npa)
REAL :: DTMAX, DTMAX_GL, DTMAXEXP, REST
REAL :: LAMBDA(2), KTMP(3), TMP(3)
REAL :: THETA_L(3), BET1(3), BETAHAT(3)
REAL :: KELEM(3,NE), FLALL(3,NE)
REAL :: KKSUM(npa), ST(npa)
REAL :: NM(NE)
!/S CALL STRACE (IENT, 'W3XYPFSN')
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
DTMAX = DBLE(10.E10)
!/REF1 IOBPDR(:)=(1-IOBP(:))*(1-IOBPD(ITH,:))
DO IE = 1, NE
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
LAMBDA(1) = ONESIXTH *(C(I1,1)+C(I2,1)+C(I3,1)) ! Linearized advection speed in X and Y direction
LAMBDA(2) = ONESIXTH *(C(I1,2)+C(I2,2)+C(I3,2))
KELEM(1,IE) = LAMBDA(1) * PDLIB_IEN(1,IE) + LAMBDA(2) * PDLIB_IEN(2,IE) ! K-Values - so called Flux Jacobians
KELEM(2,IE) = LAMBDA(1) * PDLIB_IEN(3,IE) + LAMBDA(2) * PDLIB_IEN(4,IE)
KELEM(3,IE) = LAMBDA(1) * PDLIB_IEN(5,IE) + LAMBDA(2) * PDLIB_IEN(6,IE)
KTMP = KELEM(:,IE) ! Copy
NM(IE) = - 1.D0/MIN(-THR,SUM(MIN(ZERO,KTMP))) ! N-Values
KELEM(:,IE) = MAX(ZERO,KTMP)
FL11 = C(I2,1) * PDLIB_IEN(1,IE) + C(I2,2) * PDLIB_IEN(2,IE) ! Weights for Simpson Integration
FL12 = C(I3,1) * PDLIB_IEN(1,IE) + C(I3,2) * PDLIB_IEN(2,IE)
FL21 = C(I3,1) * PDLIB_IEN(3,IE) + C(I3,2) * PDLIB_IEN(4,IE)
FL22 = C(I1,1) * PDLIB_IEN(3,IE) + C(I1,2) * PDLIB_IEN(4,IE)
FL31 = C(I1,1) * PDLIB_IEN(5,IE) + C(I1,2) * PDLIB_IEN(6,IE)
FL32 = C(I2,1) * PDLIB_IEN(5,IE) + C(I2,2) * PDLIB_IEN(6,IE)
FL111 = 2.d0*FL11+FL12
FL112 = 2.d0*FL12+FL11
FL211 = 2.d0*FL21+FL22
FL212 = 2.d0*FL22+FL21
FL311 = 2.d0*FL31+FL32
FL312 = 2.d0*FL32+FL31
FLALL(1,IE) = (FL311 + FL212)! * ONESIXTH + KELEM(1,IE)
FLALL(2,IE) = (FL111 + FL312)! * ONESIXTH + KELEM(2,IE)
FLALL(3,IE) = (FL211 + FL112)! * ONESIXTH + KELEM(3,IE)
END DO
IF (LCALC) THEN
KKSUM = ZERO
DO IE = 1, NE
NI = INE(:,IE)
KKSUM(NI) = KKSUM(NI) + KELEM(:,IE)
END DO
DO IP = 1, npa
DTMAXEXP = PDLIB_SI(IP)/MAX(DBLE(10.E-10),KKSUM(IP)*IOBDP(IP))
DTMAX = MIN( DTMAX, DTMAXEXP)
CFLXYMAX(IP) = MAX(CFLXYMAX(IP),DBLE(DT)/DTMAXEXP)
END DO ! IP
CALL MPI_ALLREDUCE(DTMAX,DTMAX_GL,1,rtype,MPI_MIN,MPI_COMM_WCMP,ierr)
CFLXY = DBLE(DT)/DTMAX_GL
REST = ABS(MOD(CFLXY,1.0d0))
IF (REST .LT. THR) THEN
ITER(IK,ITH) = ABS(NINT(CFLXY))
ELSE IF (REST .GT. THR .AND. REST .LT. 0.5d0) THEN
ITER(IK,ITH) = ABS(NINT(CFLXY)) + 1
ELSE
ITER(IK,ITH) = ABS(NINT(CFLXY))
END IF
END IF
DO IP = 1, npa
DTSI(IP) = DBLE(DT)/DBLE(ITER(IK,ITH))/PDLIB_SI(IP) ! Some precalculations for the time integration.
END DO
DO IT = 1, ITER(IK,ITH)
U = DBLE(AC)
ST = ZERO
DO IE = 1, NE
NI = INE(:,IE)
FT = -ONESIXTH*DOT_PRODUCT(U(NI),FLALL(:,IE))
UTILDE = NM(IE) * ( DOT_PRODUCT(KELEM(:,IE),U(NI)) - FT )
THETA_L(:) = KELEM(:,IE) * (U(NI) - UTILDE)
IF (ABS(FT) .GT. 0.0d0) THEN
BET1(:) = THETA_L(:)/FT
IF (ANY( BET1 .LT. 0.0d0) ) THEN
BETAHAT(1) = BET1(1) + 0.5d0 * BET1(2)
BETAHAT(2) = BET1(2) + 0.5d0 * BET1(3)
BETAHAT(3) = BET1(3) + 0.5d0 * BET1(1)
BET1(1) = MAX(ZERO,MIN(BETAHAT(1),1.d0-BETAHAT(2),1.d0))
BET1(2) = MAX(ZERO,MIN(BETAHAT(2),1.d0-BETAHAT(3),1.d0))
BET1(3) = MAX(ZERO,MIN(BETAHAT(3),1.d0-BETAHAT(1),1.d0))
THETA_L(:) = FT * BET1
END IF
ELSE
THETA_L(:) = ZERO
END IF
ST(NI) = ST(NI) + THETA_L ! the 2nd term are the theta values of each node ...
END DO
DO IP = 1, npa
IP_glob=iplg(IP)
U(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*(1-IOBPA(IP_glob)))*DBLE(IOBPD(ITH,IP_glob))*IOBDP(IP_glob)
!/REF1 IF (REFPARS(3).LT.0.5.AND.IOBPD(ITH,IP_glob).EQ.0.AND.IOBPA(IP_glob).EQ.0) U(IP) = AC(IP) ! restores reflected boundary values
END DO
AC = REAL(U)
!
! 5 Update boundaries ... this should be implemented differently ... it is better to omit any if clause in this loop ...
!
IF ( FLBPI ) THEN
RD1=RD10 - DT * REAL(ITER(IK,ITH)-IT)/REAL(ITER(IK,ITH))
RD2=RD20
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
ELSE
RD1 = 0.
RD2 = 1.
END IF
!
! NB: this treatment of the open boundary (time interpolation) is different from
! the constant boundary in the structured grids ... which restores the boundary
! to the initial value: IF ( MAPSTA(IXY).EQ.2 ) VQ(IXY) = AQ(IXY)
! Why this difference ?
!
DO IBI=1, NBI
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) ) &
/ CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
END IF
ENDDO
END IF
CALL PDLIB_exchange1DREAL(AC)
END DO ! IT
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE HACK_CHECK(string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY : NK, NTH
USE W3WDATMD, ONLY : VA
USE W3GDATMD, ONLY : NSPEC, NX, NY, NSEAL
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
CHARACTER(*), INTENT(in) :: string
INTEGER ITH_F, IK
INTEGER ITH, ISP, JSEA
REAL eVal, eErr
!/S CALL STRACE (IENT, 'HACK_CHECK')
ITH_F=4
WRITE(740+IAPROC,*) 'HACK_CHECK, begin'
DO ITH=1,NTH
IF (ITH .eq. ITH_F) THEN
eVal=0.1
ELSE
eVal=0
END IF
DO IK=1,NK
ISP=ITH + (IK-1)*NTH
DO JSEA=1,NSEAL
eErr=abs(VA(ISP,JSEA) - eVal)
IF (eErr .gt. 0.01) THEN
WRITE(740+IAPROC,*) 'HACK CHECK, str=', string
WRITE(740+IAPROC,*) 'ITH=', ITH
WRITE(740+IAPROC,*) 'IK=', IK
WRITE(740+IAPROC,*) 'ISP=', ISP
WRITE(740+IAPROC,*) 'JSEA=', JSEA
WRITE(740+IAPROC,*) 'eVal=', eVal
WRITE(740+IAPROC,*) 'VA(ISP,JSEA)=', VA(ISP,JSEA)
FLUSH(740+IAPROC)
END IF
END DO
END DO
END DO
WRITE(740+IAPROC,*) 'HACK_CHECK, end'
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE GET_SCAL_INTEGRAL_COH_R4(V, eSum)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY : NSEAL
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
REAL, INTENT(in) :: V(NSEAL)
CHARACTER (len = *), PARAMETER :: string = "Test a view"
REAL, INTENT(out) :: eSum
REAL :: V8(NSEAL)
LOGICAL :: PrintFullValue = .FALSE.
LOGICAL :: PrintBasicData = .FALSE.
!
!/S CALL STRACE (IENT, 'GET_SCAL_INTEGRAL_COH_R4')
V8 = DBLE(V)
CALL SCAL_INTEGRAL_PRINT_GENERAL(eSum, V8, string, PrintFullValue, PrintBasicData)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE SCAL_INTEGRAL_PRINT_R8(V, string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel Debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY : NSEAL
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
REAL, INTENT(in) :: V(NSEAL)
CHARACTER(*), INTENT(in) :: string
REAL :: V8(NSEAL)
LOGICAL :: PrintFullValue = .FALSE.
LOGICAL :: PrintBasicData = .TRUE.
REAL :: eSum
!
!/S CALL STRACE (IENT, 'SCAL_INTEGRAL_PRINT_R8')
V8 = DBLE(V)
CALL SCAL_INTEGRAL_PRINT_GENERAL(eSum, V8, string, PrintFullValue, PrintBasicData)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE SCAL_INTEGRAL_PRINT_R4(V, string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel Debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY : NSEAL
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
REAL, INTENT(in) :: V(NSEAL)
CHARACTER(*), INTENT(in) :: string
LOGICAL :: PrintFullValue = .FALSE.
LOGICAL :: PrintBasicData = .TRUE.
REAL :: eSum
!
!/S CALL STRACE (IENT, 'SCAL_INTEGRAL_PRINT_R4')
CALL SCAL_INTEGRAL_PRINT_GENERAL(eSum, V, string, PrintFullValue, PrintBasicData)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE SCAL_INTEGRAL_PRINT_GENERAL(eSum, V, string, PrintFullValue, PrintBasicData)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel Debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY : NK, NTH, FTE
USE W3GDATMD, ONLY : NSPEC, NX, NY, NSEAL, MAPFS
USE W3ADATMD, ONLY : MPI_COMM_WCMP
USE W3GDATMD, ONLY : GTYPE, UNGTYPE
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
use yowDatapool, ONLY: rtype, istatus
USE YOWNODEPOOL, ONLY: npa, iplg
USE W3PARALL, ONLY: INIT_GET_ISEA
!/
IMPLICIT NONE
INCLUDE "mpif.h"
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!! INCLUDE "mpif.h"
!
REAL, INTENT(out) :: eSum
REAL, INTENT(in) :: V(NSEAL)
CHARACTER(*), INTENT(in) :: string
LOGICAL, INTENT(in) :: PrintFullValue, PrintBasicData
!
REAL Vcoll(NX), rVect(NX)
REAL CoherencyError, eVal1, eVal2, eErr
INTEGER rStatus(NX), Status(NX)
INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
INTEGER nbIncorr
INTEGER ITH, IK
!/S CALL STRACE (IENT, 'SCAL_INTEGRAL_PRINT_GENERAL')
WRITE(740+IAPROC,*) 'IAPROC=', IAPROC
WRITE(740+IAPROC,*) 'NAPROC=', NAPROC
WRITE(740+IAPROC,*) 'NTPROC=', NTPROC
WRITE(740+IAPROC,*) 'Beginning of routine'
WRITE(740+IAPROC,*) 'string=', string
WRITE(740+IAPROC,*) 'NX=', NX
FLUSH(740+IAPROC)
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
IF (GTYPE .ne. UNGTYPE) THEN
RETURN
END IF
Vcoll=0
Status=0
! WRITE(740+IAPROC,*) 'PrintFullValue=', PrintFullValue
! FLUSH(740+IAPROC)
DO JSEA=1,NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA=MAPFS(1,IP_glob)
Vcoll(IP_glob)=V(JSEA)
Status(IP_glob)=1
END DO
!
! Now find global arrays
!
CoherencyError=0
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
WRITE(740+IAPROC,*) 'Before MPI_RECV, 1'
FLUSH(740+IAPROC)
CALL MPI_RECV(rVect,NX,MPI_REAL8, iProc-1, 370, MPI_COMM_WCMP, istatus, ierr)
WRITE(740+IAPROC,*) 'Before MPI_RECV, 2'
FLUSH(740+IAPROC)
CALL MPI_RECV(rStatus,NX,MPI_INTEGER, iProc-1, 430, MPI_COMM_WCMP, istatus, ierr)
WRITE(740+IAPROC,*) 'After MPI_RECV, 2'
FLUSH(740+IAPROC)
DO I=1,NX
IF (rStatus(I) .eq. 1) THEN
eVal1=Vcoll(I)
eVal2=rVect(I)
Vcoll(I)=rVect(I)
IF (Status(I) .eq. 1) THEN
eErr=abs(eVal1 - eVal2)
CoherencyError = CoherencyError + eErr
ELSE
Vcoll(I)=eVal2
END IF
Status(I)=1
END IF
END DO
END DO
ELSE
WRITE(740+IAPROC,*) 'Before MPI_SEND, 1'
FLUSH(740+IAPROC)
CALL MPI_SEND(Vcoll,NX,MPI_REAL8, 0, 370, MPI_COMM_WCMP, ierr)
WRITE(740+IAPROC,*) 'Before MPI_SEND, 2'
FLUSH(740+IAPROC)
CALL MPI_SEND(Status,NX,MPI_INTEGER, 0, 430, MPI_COMM_WCMP, ierr)
WRITE(740+IAPROC,*) 'After MPI_SEND, 2'
FLUSH(740+IAPROC)
END IF
eSum=sum(Vcoll)
IF ((IAPROC .eq. 1).and. PrintBasicData) THEN
nbIncorr=0
DO IX=1,NX
ISEA=MAPFS(1,IX)
IF (ISEA .gt. 0) THEN
IF (Status(IX) .eq. 0) THEN
nbIncorr=nbIncorr+1
END IF
END IF
END DO
IF (nbIncorr .gt. 0) THEN
WRITE(*,*) ' nbIncorr=', nbIncorr
WRITE(*,*) ' NX=', NX
WRITE(*,*) ' npa=', npa
STOP
END IF
WRITE(740+IAPROC,*) 'VA_INTEGRAL sum,coh=', sum(Vcoll), CoherencyError, TRIM(string)
FLUSH(740+IAPROC)
IF (PrintFullValue) THEN
WRITE(740+IAPROC,*) 'Value of V at nodes'
DO IX=1,NX
WRITE(740+IAPROC,*) 'IX=', IX, ' V=', Vcoll(IX)
END DO
FLUSH(740+IAPROC)
END IF
END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE ALL_VAOLD_INTEGRAL_PRINT(string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel Debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NSEAL
USE W3WDATMD, ONLY : VAOLD
USE W3ODATMD, ONLY : IAPROC
USE W3GDATMD, ONLY : NSPEC
IMPLICIT NONE
REAL :: FIELD(NSPEC,NSEAL)
CHARACTER(*), INTENT(in) :: string
! LOGICAL :: PrintHsNode = .FALSE.
LOGICAL :: PrintHs = .TRUE.
LOGICAL :: PrintHsNode = .FALSE.
INTEGER :: ScalMeth = 4
INTEGER ISPEC, JSEA
DO JSEA=1,NSEAL
DO ISPEC=1,NSPEC
FIELD(ISPEC,JSEA) = VAOLD(ISPEC,JSEA)
END DO
END DO
CALL ALL_FIELD_INTEGRAL_PRINT_GENERAL(FIELD, string, PrintHs, PrintHsNode, ScalMeth)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE ALL_VA_INTEGRAL_PRINT(IMOD, string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NSEAL, NSEA, NX, NY
USE W3WDATMD, ONLY : VA
USE W3ODATMD, ONLY : IAPROC, NAPROC
USE W3GDATMD, ONLY : NSPEC, GRIDS, GTYPE, UNGTYPE
IMPLICIT NONE
REAL :: FIELD(NSPEC,NSEAL)
INTEGER, INTENT(in) :: IMOD
CHARACTER(*), INTENT(in) :: string
! LOGICAL :: PrintHsNode = .FALSE.
LOGICAL :: PrintHs = .TRUE.
LOGICAL :: PrintHsNode = .FALSE.
INTEGER :: ScalMeth = 4
INTEGER ISTAT
INTEGER ISPEC, JSEA
IF (GRIDS(IMOD)%GTYPE .ne. UNGTYPE) THEN
RETURN
END IF
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
! WRITE(*,*) 'IMOD=', IMOD, ' NSEA=', NSEA, ' NSEAL=', NSEAL
! WRITE(*,*) 'GRIDS(IMOD)%NSEA = ', GRIDS(IMOD)%NSEA
! WRITE(*,*) 'GRIDS(IMOD)%NSEAL = ', GRIDS(IMOD)%NSEAL
! WRITE(*,*) 'NX=', NX, ' NY=', NY
! WRITE(*,*) 'Bonjour'
! WRITE(740+IAPROC,*) 'Bonjour, ALL_VA_INTEGRAL'
! FLUSH(740+IAPROC)
! RETURN
WRITE(740+IAPROC,*) 'Entering ALL_INTEGRAL_PRINT, NSEAL=', NSEAL
FLUSH(740+IAPROC)
DO JSEA=1,NSEAL
! WRITE(740+IAPROC,*) 'JSEA=', JSEA
! FLUSH(740+IAPROC)
DO ISPEC=1,NSPEC
FIELD(ISPEC,JSEA) = VA(ISPEC,JSEA)
END DO
END DO
WRITE(740+IAPROC,*) 'Before call to ALL_FIELD_INTEGRAL_PRINT_GENERAL'
WRITE(740+IAPROC,*) 'NSPEC=', NSPEC, ' NX=', NX
FLUSH(740+IAPROC)
CALL ALL_FIELD_INTEGRAL_PRINT_GENERAL(FIELD, string, PrintHs, PrintHsNode, ScalMeth)
WRITE(740+IAPROC,*) 'After call to ALL_FIELD_INTEGRAL_PRINT_GENERAL'
FLUSH(740+IAPROC)
! IF (NSEAL >= 40) THEN
! WRITE(740+IAPROC,*) 'min/max/sum(VA(:,TESTNODE))=', minval(VA(:,TESTNODE)), maxval(VA(:,TESTNODE)), sum(VA(:,TESTNODE))
! FLUSH(740+IAPROC)
! END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE ALL_FIELD_INTEGRAL_PRINT(FIELD, string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NSEAL
USE W3WDATMD, ONLY : VA
USE W3ODATMD, ONLY : IAPROC
USE W3GDATMD, ONLY : NSPEC
IMPLICIT NONE
REAL, INTENT(in) :: FIELD(NSPEC,NSEAL)
CHARACTER(*), INTENT(in) :: string
! LOGICAL :: PrintHsNode = .FALSE.
LOGICAL :: PrintHs = .FALSE.
LOGICAL :: PrintHsNode = .FALSE.
INTEGER :: ScalMeth = 1
CALL ALL_FIELD_INTEGRAL_PRINT_GENERAL(FIELD, string, PrintHs, PrintHsNode, ScalMeth)
!/
!/ End of JACOBI_FINALIZE -------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE ALL_FIELD_INTEGRAL_PRINT_GENERAL(FIELD, string, PrintHs, PrintHsNode, ScalMeth)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : CLATS, SIG
USE W3ADATMD, ONLY : CG
USE W3GDATMD, ONLY : NK, NTH, FTE
USE W3WDATMD, ONLY : VA
USE W3GDATMD, ONLY : NSPEC, NX, NY, NSEAL, MAPFS
USE W3ADATMD, ONLY : MPI_COMM_WCMP
USE W3GDATMD, ONLY : GTYPE, UNGTYPE
USE W3GDATMD, ONLY : DTH, DSII
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
USE CONSTANTS, ONLY: TPIINV
use yowDatapool, ONLY: rtype, istatus
USE YOWNODEPOOL, ONLY: npa, iplg
USE W3PARALL, ONLY: INIT_GET_ISEA
IMPLICIT NONE
INCLUDE "mpif.h"
REAL sumHS, avgHS, maxHS
REAL ET, HSIG, EBND
REAL, INTENT(in) :: FIELD(NSPEC,NSEAL)
CHARACTER(*), INTENT(in) :: string
LOGICAL, INTENT(in) :: PrintHs, PrintHsNode
INTEGER, INTENT(in) :: ScalMeth
REAL Vcoll(NSPEC,NX)
REAL VcollExp(NSPEC*NX)
REAL rVect(NSPEC*NX)
REAL Vhs(NX), rHs(NX)
REAL VminCG(NX), rMinCG(NX)
REAL VmaxCG(NX), rMaxCG(NX)
REAL eFact
REAL CoherencyError, eVal1, eVal2, eErr
INTEGER rStatus(NX), Status(NX)
INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
INTEGER nbIncorr
INTEGER ITH, IK
! WRITE(740+IAPROC,*) 'IAPROC=', IAPROC
! WRITE(740+IAPROC,*) 'NAPROC=', NAPROC
! WRITE(740+IAPROC,*) 'NTPROC=', NTPROC
! WRITE(740+IAPROC,*) 'Beginning of routine'
! WRITE(740+IAPROC,*) 'string=', string
! FLUSH(740+IAPROC)
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
IF (GTYPE .ne. UNGTYPE) THEN
RETURN
END IF
VcollExp=0
Vhs=0
Status=0
WRITE(740+IAPROC,*) 'ScalMeth=', ScalMeth
FLUSH(740+IAPROC)
WRITE(740+IAPROC,*) 'min(CLATS)=', minval(CLATS)
FLUSH(740+IAPROC)
WRITE(740+IAPROC,*) 'max(CLATS)=', maxval(CLATS)
FLUSH(740+IAPROC)
WRITE(740+IAPROC,*) 'min(CG)=', minval(CG)
FLUSH(740+IAPROC)
WRITE(740+IAPROC,*) 'max(CG)=', maxval(CG)
FLUSH(740+IAPROC)
WRITE(740+IAPROC,*) 'FTE=', FTE
FLUSH(740+IAPROC)
DO JSEA=1,NSEAL
! WRITE(740+IAPROC,*) 'JSEA=', JSEA
! FLUSH(740+IAPROC)
IP = JSEA
! WRITE(740+IAPROC,*) 'IP=', IP
! FLUSH(740+IAPROC)
IP_glob = iplg(IP)
! WRITE(740+IAPROC,*) 'IP_glob=', IP_glob
! FLUSH(740+IAPROC)
ISEA=MAPFS(1,IP_glob)
! WRITE(740+IAPROC,*) 'ISEA=', ISEA
! FLUSH(740+IAPROC)
DO ISP=1,NSPEC
VcollExp(ISP+NSPEC*(IP_glob-1))=FIELD(ISP,JSEA)
END DO
! WRITE(740+IAPROC,*) 'After VcollExp assignation'
! FLUSH(740+IAPROC)
!
ET = 0
DO IK=1, NK
EBND = 0.
DO ITH=1, NTH
ISP = ITH + (IK-1)*NTH
IF (ScalMeth .eq. 1) THEN
eFact=1
END IF
IF (ScalMeth .eq. 2) THEN
eFact=SIG(IK)/CG(IK,ISEA)
END IF
IF (ScalMeth .eq. 3) THEN
eFact=CLATS(ISEA)/CG(IK,ISEA)
END IF
IF (ScalMeth .eq. 4) THEN
eFact=DTH * DSII(IK) * TPIINV
END IF
EBND = EBND + FIELD(ISP,JSEA)*eFact
END DO
ET = ET + EBND
END DO
! WRITE(740+IAPROC,*) 'After IK loop'
! FLUSH(740+IAPROC)
ET = ET + FTE *EBND
IF (PrintHs) THEN
HSIG = 4. * SQRT ( ET )
Vhs(IP_glob)=HSIG
END IF
VminCG(IP_glob)=minval(CG(:,ISEA))
VmaxCG(IP_glob)=maxval(CG(:,ISEA))
Status(IP_glob)=1
! WRITE(740+IAPROC,*) 'After Status assignation'
! FLUSH(740+IAPROC)
END DO
WRITE(740+IAPROC,*) 'After status settings'
FLUSH(740+IAPROC)
!
! Now find global arrays
!
CoherencyError=0
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
WRITE(740+IAPROC,*) '1: iProc=', iProc
FLUSH(740+IAPROC)
CALL MPI_RECV(rVect,NSPEC*NX,MPI_REAL, iProc-1, 37, MPI_COMM_WCMP, istatus, ierr)
CALL MPI_RECV(rHs,NX,MPI_REAL, iProc-1, 38, MPI_COMM_WCMP, istatus, ierr)
CALL MPI_RECV(rMinCG,NX,MPI_REAL, iProc-1, 39, MPI_COMM_WCMP, istatus, ierr)
CALL MPI_RECV(rMaxCG,NX,MPI_REAL, iProc-1, 40, MPI_COMM_WCMP, istatus, ierr)
WRITE(740+IAPROC,*) '2: iProc=', iProc
FLUSH(740+IAPROC)
CALL MPI_RECV(rStatus,NX,MPI_INTEGER, iProc-1, 43, MPI_COMM_WCMP, istatus, ierr)
WRITE(740+IAPROC,*) '3: iProc=', iProc
FLUSH(740+IAPROC)
DO I=1,NX
IF (rStatus(I) .eq. 1) THEN
DO ISP=1,NSPEC
eVal1=VcollExp(ISP+NSPEC*(I-1))
eVal2=rVect(ISP+NSPEC*(I-1))
VcollExp(ISP+NSPEC*(I-1))=rVect(ISP+NSPEC*(I-1))
Vhs(I)=rHs(I)
VminCG(I)=rMinCG(I)
VmaxCG(I)=rMaxCG(I)
IF (Status(I) .eq. 1) THEN
eErr=abs(eVal1 - eVal2)
CoherencyError = CoherencyError + eErr
ELSE
VcollExp(ISP+NSPEC*(I-1))=eVal2
END IF
END DO
Status(I)=1
END IF
END DO
END DO
ELSE
WRITE(740+IAPROC,*) 'Second case 1'
FLUSH(740+IAPROC)
CALL MPI_SEND(VcollExp,NSPEC*NX,MPI_REAL, 0, 37, MPI_COMM_WCMP, ierr)
CALL MPI_SEND(Vhs,NX,MPI_REAL, 0, 38, MPI_COMM_WCMP, ierr)
CALL MPI_SEND(VminCG,NX,MPI_REAL, 0, 39, MPI_COMM_WCMP, ierr)
CALL MPI_SEND(VmaxCG,NX,MPI_REAL, 0, 40, MPI_COMM_WCMP, ierr)
WRITE(740+IAPROC,*) 'Second case 2'
FLUSH(740+IAPROC)
CALL MPI_SEND(Status,NX,MPI_INTEGER, 0, 43, MPI_COMM_WCMP, ierr)
WRITE(740+IAPROC,*) 'Second case 3'
FLUSH(740+IAPROC)
END IF
IF (IAPROC .eq. 1) THEN
DO I=1,NX
DO ISP=1,NSPEC
Vcoll(ISP,I)=VcollExp(ISP + NSPEC*(I-1))
END DO
END DO
nbIncorr=0
DO IX=1,NX
ISEA=MAPFS(1,IX)
IF (ISEA .gt. 0) THEN
IF (Status(IX) .eq. 0) THEN
nbIncorr=nbIncorr+1
END IF
END IF
END DO
IF (nbIncorr .gt. 0) THEN
WRITE(*,*) ' nbIncorr=', nbIncorr
WRITE(*,*) ' NX=', NX
WRITE(*,*) ' npa=', npa
STOP
END IF
WRITE(740+IAPROC,*) 'FIELD_INTEGRAL sum,coh=', sum(Vcoll), CoherencyError, TRIM(string)
IF (PrintHs) THEN
sumHS=sum(Vhs)
avgHS=sumHS/REAL(NX)
maxHS=maxval(Vhs)
WRITE(740+IAPROC,*) 'HS max,avg=', maxHS, avgHS
END IF
IF (PrintHsNode) THEN
WRITE(740+IAPROC,*) 'Value of HS at nodes'
DO IX=1,NX
! WRITE(740+IAPROC,*) 'IX=', IX, ' hs/cg(mm)=', Vhs(IX), VminCG(IX), VmaxCG(IX)
WRITE(740+IAPROC,*) 'IX=', IX, ' hs/sumV=', Vhs(IX), sum(Vcoll(:,IX))
END DO
FLUSH(740+IAPROC)
END IF
FLUSH(740+IAPROC)
END IF
!/
!/ End of JACOBI_FINALIZE -------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE TEST_MPI_STATUS(string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3ADATMD, ONLY : MPI_COMM_WCMP
USE W3GDATMD, ONLY : GTYPE, UNGTYPE
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
use yowDatapool, ONLY: rtype, istatus
IMPLICIT NONE
INCLUDE "mpif.h"
CHARACTER(*), INTENT(in) :: string
REAL VcollExp(1)
REAL rVect(1)
INTEGER iProc, ierr
WRITE(740+IAPROC,*) 'TEST_MPI_STATUS, at string=', string
FLUSH(740+IAPROC)
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
WRITE(740+IAPROC,*) 'After status settings'
FLUSH(740+IAPROC)
!
! Now find global arrays
!
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
WRITE(740+IAPROC,*) '1: iProc=', iProc
FLUSH(740+IAPROC)
CALL MPI_RECV(rVect,1,MPI_REAL, iProc-1, 37, MPI_COMM_WCMP, istatus, ierr)
END DO
ELSE
WRITE(740+IAPROC,*) 'Second case 1'
FLUSH(740+IAPROC)
CALL MPI_SEND(VcollExp,1,MPI_REAL, 0, 37, MPI_COMM_WCMP, ierr)
END IF
WRITE(740+IAPROC,*) 'Leaving the TEST_MPI_STATUS'
FLUSH(740+IAPROC)
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
!* maxidx should be "np" or "npa" *
SUBROUTINE CHECK_ARRAY_INTEGRAL_NX_R8(TheARR, string, maxidx)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NSPEC
USE YOWNODEPOOL, ONLY: npa
CHARACTER(*), INTENT(in) :: string
INTEGER, INTENT(in) :: maxidx
REAL, INTENT(in) :: TheARR(NSPEC, npa)
LOGICAL :: PrintMinISP, LocalizeMaximum
PrintMinISP = .false.
LocalizeMaximum = .false.
CALL CHECK_ARRAY_INTEGRAL_NX_R8_MaxFunct(TheARR, string, maxidx, PrintMinISP, LocalizeMaximum)
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
!* maxidx should be "np" or "npa" *
SUBROUTINE CHECK_ARRAY_INTEGRAL_NX_R8_MaxFunct(TheARR, string, maxidx, PrintMinISP, LocalizeMaximum)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NK, NTH
USE W3WDATMD, ONLY : VA
USE W3GDATMD, ONLY : NSPEC, NX, NY, NSEAL, MAPFS
USE W3ADATMD, ONLY : MPI_COMM_WCMP
USE W3GDATMD, ONLY : GTYPE, UNGTYPE
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
use yowDatapool, ONLY: rtype, istatus
USE YOWNODEPOOL, ONLY: npa, iplg
USE W3PARALL, ONLY: INIT_GET_ISEA
IMPLICIT NONE
INCLUDE "mpif.h"
CHARACTER(*), INTENT(in) :: string
INTEGER, INTENT(in) :: maxidx
REAL, INTENT(in) :: TheARR(NSPEC, npa)
LOGICAL, INTENT(in) :: PrintMinISP, LocalizeMaximum
!
REAL Vcoll(NSPEC,NX), VcollExp(NSPEC*NX), rVect(NSPEC*NX)
REAL CoherencyError, eVal1, eVal2, eErr
INTEGER rStatus(NX), Status(NX)
INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
REAL :: mval, eVal, eSum
REAL :: TheMax, TheSum, TheNb, TheAvg
REAL :: eFact, Threshold
LOGICAL :: IsFirst
INTEGER nbIncorr
INTEGER ITH, IK
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
IF (GTYPE .ne. UNGTYPE) THEN
RETURN
END IF
VcollExp=0
Status=0
DO IP=1,maxidx
IP_glob=iplg(IP)
DO ISP=1,NSPEC
VcollExp(ISP+NSPEC*(IP_glob-1))=TheARR(ISP,IP)
END DO
Status(IP_glob)=1
END DO
!
! Now find global arrays
!
CoherencyError=0
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
CALL MPI_RECV(rVect ,NSPEC*NX,MPI_REAL , iProc-1, 37, MPI_COMM_WCMP, istatus, ierr)
CALL MPI_RECV(rStatus,NX ,MPI_INTEGER, iProc-1, 43, MPI_COMM_WCMP, istatus, ierr)
DO I=1,NX
IF (rStatus(I) .eq. 1) THEN
DO ISP=1,NSPEC
eVal1=VcollExp(ISP+NSPEC*(I-1))
eVal2=rVect(ISP+NSPEC*(I-1))
VcollExp(ISP+NSPEC*(I-1))=rVect(ISP+NSPEC*(I-1))
IF (Status(I) .eq. 1) THEN
eErr=abs(eVal1 - eVal2)
CoherencyError = CoherencyError + eErr
ELSE
VcollExp(ISP+NSPEC*(I-1))=eVal2
END IF
END DO
Status(I)=1
END IF
END DO
END DO
ELSE
CALL MPI_SEND(VcollExp,NSPEC*NX,MPI_REAL , 0, 37, MPI_COMM_WCMP, ierr)
CALL MPI_SEND(Status ,NX ,MPI_INTEGER, 0, 43, MPI_COMM_WCMP, ierr)
END IF
IF (IAPROC .eq. 1) THEN
DO I=1,NX
DO ISP=1,NSPEC
Vcoll(ISP,I)=VcollExp(ISP + NSPEC*(I-1))
END DO
END DO
nbIncorr=0
DO IX=1,NX
ISEA=MAPFS(1,IX)
IF (ISEA .gt. 0) THEN
IF (Status(IX) .eq. 0) THEN
nbIncorr=nbIncorr+1
END IF
END IF
END DO
IF (nbIncorr .gt. 0) THEN
WRITE(*,*) ' nbIncorr=', nbIncorr
WRITE(*,*) ' NX=', NX
WRITE(*,*) ' npa=', npa
STOP
END IF
WRITE(740+IAPROC,*) 'ARRAY_NX sum,coh=', sum(Vcoll), CoherencyError, TRIM(string)
IF (PrintMinISP) THEN
DO ISP=1,NSPEC
IsFirst=.true.
eSum=0
DO IP=1,maxidx
eVal=abs(Vcoll(ISP, IP))
eSum=eSum + eVal
IF (IsFirst.eqv. .true.) then
mval=eVal
ELSE
IF (eVal .lt. mval) THEN
mval=eVal
ENDIF
ENDIF
IsFirst=.false.
END DO
WRITE(740+IAPROC,*) 'ISP=', ISP, ' mval/sum=', mval, eSum
END DO
FLUSH(740+IAPROC)
END IF
IF (LocalizeMaximum) THEN
TheMax=0
TheNb=0
TheSum=0
DO IP=1,maxidx
DO ISP=1,NSPEC
eVal = abs(Vcoll(ISP, IP))
TheSum = TheSum + eVal
TheNb = TheNb + 1
IF (eVal .gt. TheMax) THEN
TheMax=eVal
END IF
END DO
END DO
TheAvg = TheSum / TheNb
WRITE(740+IAPROC,*) 'TheAvg/TheMax=', TheAvg, TheMax
eFact=0.5
Threshold=eFact * TheMax
DO IP=1,maxidx
DO ISP=1,NSPEC
eVal = abs(Vcoll(ISP, IP))
IF (eVal .gt. Threshold) THEN
WRITE(740+IAPROC,*) 'ISP/IP/val=', ISP, IP, eVal
END IF
END DO
END DO
FLUSH(740+IAPROC)
END IF
END IF
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE CHECK_ARRAY_INTEGRAL_NX_R8_MaxFunct
!/ ------------------------------------------------------------------- /
SUBROUTINE COLLECT_AND_PRINT(V, ifile)
!/ ------------------------------------------------------------------- /
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NX, NSEAL, NSEA
USE W3ADATMD, ONLY: MPI_COMM_WCMP
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
USE YOWNODEPOOL, ONLY: iplg, npa
use yowDatapool, ONLY: rtype, istatus
IMPLICIT NONE
INCLUDE "mpif.h"
INTEGER, INTENT(in) :: ifile
REAL, INTENT(in) :: V(npa)
REAL Vcoll(NX), rVect(NX)
INTEGER rStatus(NX), Status(NX)
INTEGER JSEA, ISEA, iProc, I, ierr
INTEGER nbIncorr
Vcoll=0
Status=0
DO JSEA=1,npa
ISEA=iplg(JSEA)
Vcoll(ISEA)=V(JSEA)
Status(ISEA)=1
END DO
!
! Now find global arrays
!
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
CALL MPI_RECV(rVect,NX,rtype, iProc-1, 19, MPI_COMM_WCMP, istatus, ierr)
CALL MPI_RECV(rStatus,NX,MPI_INTEGER, iProc-1, 23, MPI_COMM_WCMP, istatus, ierr)
DO I=1,NX
IF (rStatus(I) .eq. 1) THEN
Vcoll(I)=rVect(I)
Status(I)=1
END IF
END DO
END DO
ELSE
CALL MPI_SEND(Vcoll,NX,rtype, 0, 19, MPI_COMM_WCMP, ierr)
CALL MPI_SEND(Status,NX,MPI_INTEGER, 0, 23, MPI_COMM_WCMP, ierr)
END IF
IF (IAPROC .eq. 1) THEN
nbIncorr=0
DO I=1,NX
IF (Status(I) .eq. 0) THEN
nbIncorr=nbIncorr+1
END IF
END DO
IF (nbIncorr .gt. 0) THEN
WRITE(740+IAPROC,*) ' nbIncorr=', nbIncorr
WRITE(740+IAPROC,*) 'NX - NSEAL=', NX - NSEAL
FLUSH(740+IAPROC)
STOP
END IF
!/DEBUGSOLVER DO I=1,NX
!/DEBUGSOLVER WRITE(ifile,*) 'I=', I, ' var=', Vcoll(I)
!/DEBUGSOLVER END DO
END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE DOMAIN_INTEGRAL(V, eScal)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NX
USE W3ADATMD, ONLY: MPI_COMM_WCMP
use yowDatapool, ONLY: rtype, istatus
use YOWNODEPOOL, ONLY: np, npa
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
IMPLICIT NONE
REAL, INTENT(in) :: V(npa)
REAL, INTENT(inout) :: eScal
INTEGER IP
REAL :: rScal(1), lScal(1)
INTEGER iProc
INTEGER ierr
lScal=0
DO IP=1,np
lScal(1)=lScal(1) + V(IP)
END DO
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
CALL MPI_RECV(rScal,1,rtype, iProc-1, 19, MPI_COMM_WCMP, istatus, ierr)
lScal = lScal + rScal
END DO
DO iProc=2,NAPROC
CALL MPI_SEND(lScal,1,rtype, iProc-1, 23, MPI_COMM_WCMP, ierr)
END DO
ELSE
CALL MPI_SEND(lScal,1,rtype, 0, 19, MPI_COMM_WCMP, ierr)
CALL MPI_RECV(lScal,1,rtype, 0, 23, MPI_COMM_WCMP, istatus, ierr)
END IF
eScal=lScal(1)
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPFSFCT2 ( ISP, C, LCALC, RD10, RD20, DT, AC)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Explicit LF-FCT scheme
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NK, NTH, NX, &
IEN, CLATS, MAPSF, IOBPD, IOBP, TRIA, IOBDP
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: CG, ITER, CFLXYMAX
USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, ISBPI, BBPI0, BBPIN
USE W3TIMEMD, ONLY: DSEC21
USE W3GDATMD, ONLY: NSEAL, IOBPA
!/REF1 USE W3GDATMD, ONLY: REFPARS
USE W3ADATMD, ONLY: MPI_COMM_WCMP
use yowElementpool, ONLY: ne, INE
use YOWNODEPOOL, ONLY: PDLIB_SI, PDLIB_IEN, PDLIB_TRIA
use YOWNODEPOOL, ONLY: iplg, npa
use yowDatapool, ONLY: rtype
USE W3ODATMD, ONLY : IAPROC
USE MPI, ONLY : MPI_MIN
USE W3PARALL, ONLY : INIT_GET_JSEA_ISPROC, ONESIXTH, ZERO
USE W3PARALL, ONLY : THR
use yowExchangeModule, ONLY : PDLIB_exchange1DREAL
USE yowRankModule, ONLY : IPGL_npa
IMPLICIT NONE
INTEGER, INTENT(IN) :: ISP ! Actual Frequency/Wavenumber,
! actual Wave Direction
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 its
! X- and Y- Components,
REAL, INTENT(INOUT) :: AC(npa) ! Wave Action before and after
! advection
REAL, INTENT(IN) :: RD10, RD20 ! Time interpolation
! coefficients for boundary
! condition
LOGICAL, INTENT(IN) :: LCALC ! Switch for the calculation of
! the max. Global Time step
INTEGER :: IP, IE, POS, IT, I1, I2, I3, I, J, ITH, IK
INTEGER :: IBI, NI(3), JX
REAL :: RD1, RD2
REAL :: UTILDE
REAL :: SUMTHETA
REAL :: FL1, FL2, FL3
REAL :: FT, CFLXY
REAL :: FL11, FL12, FL21, FL22, FL31, FL32
REAL :: FL111, FL112, FL211, FL212, FL311, FL312
REAL :: DTSI(npa), U(npa), DT4AI, TMP1
REAL :: DTMAX_GL, DTMAX, DTMAXEXP, REST
REAL :: LAMBDA(2), KTMP(3), TMP(3)
REAL :: BET1(3), BETAHAT(3)
REAL :: THETA_L(3,NE), THETA_H(3,NE), THETA_ACE(3,NE), UTMP(3)
REAL :: WII(2,npa), UL(npa), USTARI(2,npa)
REAL :: PM(npa), PP(npa), UIM(npa), UIP(npa)
REAL :: KELEM(3,NE), FLALL(3,NE)
REAL :: KKSUM(npa), ST(npa), BETA
REAL :: NM(NE)
INTEGER :: ISproc, IP_glob, JSEA, ierr
REAL :: eScal
!/REF1 INTEGER(KIND=1) :: IOBPDR(NX)
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
DTMAX = DBLE(10.E10)
!/REF1 IOBPDR(:)=(1-IOBP(:))*(1-IOBPD(ITH,:))
DO IE = 1, NE
I1 = INE(1,IE) ! Index of the Element Nodes
I2 = INE(2,IE)
I3 = INE(3,IE)
LAMBDA(1) = ONESIXTH *(C(I1,1)+C(I2,1)+C(I3,1)) ! Linearized advection speed in X and Y direction
LAMBDA(2) = ONESIXTH *(C(I1,2)+C(I2,2)+C(I3,2))
KELEM(1,IE) = LAMBDA(1) * PDLIB_IEN(1,IE) + LAMBDA(2) * PDLIB_IEN(2,IE) ! K-Values - so called Flux Jacobians
KELEM(2,IE) = LAMBDA(1) * PDLIB_IEN(3,IE) + LAMBDA(2) * PDLIB_IEN(4,IE)
KELEM(3,IE) = LAMBDA(1) * PDLIB_IEN(5,IE) + LAMBDA(2) * PDLIB_IEN(6,IE)
KTMP = KELEM(:,IE) ! Copy
NM(IE) = - 1.D0/MIN(-THR,SUM(MIN(ZERO,KTMP))) ! N-Values
FL11 = C(I2,1) * PDLIB_IEN(1,IE) + C(I2,2) * PDLIB_IEN(2,IE) ! Weights for Simpson Integration
FL12 = C(I3,1) * PDLIB_IEN(1,IE) + C(I3,2) * PDLIB_IEN(2,IE)
FL21 = C(I3,1) * PDLIB_IEN(3,IE) + C(I3,2) * PDLIB_IEN(4,IE)
FL22 = C(I1,1) * PDLIB_IEN(3,IE) + C(I1,2) * PDLIB_IEN(4,IE)
FL31 = C(I1,1) * PDLIB_IEN(5,IE) + C(I1,2) * PDLIB_IEN(6,IE)
FL32 = C(I2,1) * PDLIB_IEN(5,IE) + C(I2,2) * PDLIB_IEN(6,IE)
FL111 = 2.d0*FL11+FL12
FL112 = 2.d0*FL12+FL11
FL211 = 2.d0*FL21+FL22
FL212 = 2.d0*FL22+FL21
FL311 = 2.d0*FL31+FL32
FL312 = 2.d0*FL32+FL31
FLALL(1,IE) = (FL311 + FL212)! * ONESIXTH + KELEM(1,IE)
FLALL(2,IE) = (FL111 + FL312)! * ONESIXTH + KELEM(2,IE)
FLALL(3,IE) = (FL211 + FL112)! * ONESIXTH + KELEM(3,IE)
END DO
! If the current field or water level changes estimate the iteration
! number based on the new flow field and the CFL number of the scheme
IF (LCALC) THEN
KKSUM = ZERO
DO IE = 1, NE
NI = INE(:,IE)
KKSUM(NI) = KKSUM(NI) + KELEM(:,IE)
END DO ! IE
DO IP = 1, npa
DTMAXEXP = PDLIB_SI(IP)/MAX(DBLE(10.E-10),KKSUM(IP)*IOBDP(IP))
DTMAX = MIN( DTMAX, DTMAXEXP)
CFLXYMAX(IP) = MAX(CFLXYMAX(IP),DBLE(DT)/DTMAXEXP)
END DO
CALL MPI_ALLREDUCE(DTMAX,DTMAX_GL,1,rtype,MPI_MIN,MPI_COMM_WCMP,ierr)
CFLXY = DBLE(DT)/DTMAX_GL
REST = ABS(MOD(CFLXY,1.0d0))
IF (REST .LT. THR) THEN
ITER(IK,ITH) = ABS(NINT(CFLXY))
ELSE IF (REST .GT. THR .AND. REST .LT. 0.5d0) THEN
ITER(IK,ITH) = ABS(NINT(CFLXY)) + 1
ELSE
ITER(IK,ITH) = ABS(NINT(CFLXY))
END IF
END IF ! LCALC
DT4AI = DBLE(DT)/DBLE(ITER(IK,ITH))
DTSI(:) = DT4AI/PDLIB_SI(:) ! Some precalculations for the time integration.
U = DBLE(AC) ! correct
UL = U
DO IT = 1, ITER(IK,ITH)
ST = ZERO
DO IE = 1, NE
NI = INE(:,IE)
UTMP = U(NI)
FT = - ONESIXTH*DOT_PRODUCT(UTMP,FLALL(:,IE))
TMP = MAX(ZERO,KELEM(:,IE))
UTILDE = NM(IE) * ( DOT_PRODUCT(TMP,UTMP) - FT )
THETA_L(:,IE) = TMP * ( UTMP - UTILDE )
IF (ABS(FT) .GT. THR) THEN
BET1(:) = THETA_L(:,IE)/FT
IF (ANY( BET1 .LT. 0.0d0) ) THEN
BETAHAT(1) = BET1(1) + 0.5d0 * BET1(2)
BETAHAT(2) = BET1(2) + 0.5d0 * BET1(3)
BETAHAT(3) = BET1(3) + 0.5d0 * BET1(1)
BET1(1) = MAX(ZERO,MIN(BETAHAT(1),1.d0-BETAHAT(2),1.d0))
BET1(2) = MAX(ZERO,MIN(BETAHAT(2),1.d0-BETAHAT(3),1.d0))
BET1(3) = MAX(ZERO,MIN(BETAHAT(3),1.d0-BETAHAT(1),1.d0))
THETA_L(:,IE) = FT * BET1
END IF
ELSE
THETA_L(:,IE) = ZERO
END IF
! THETA_H(:,IE) = (ONETHIRD+DT4AI/(2.d0*PDLIB_TRIA(IE)) * KELEM(:,IE))*FT ! LAX-WENDROFF
THETA_H(:,IE) = (1./3.+2./3.* KELEM(:,IE)/SUM(ABS(KELEM(:,IE))) )*FT ! CENTRAL SCHEME
! Antidiffusive residual according to the higher order nonmonotone scheme
THETA_ACE(:,IE) = ((THETA_H(:,IE) - THETA_L(:,IE))) * DT4AI/PDLIB_SI(NI)
ST(NI) = ST(NI) + THETA_L(:,IE)*DT4AI/PDLIB_SI(NI)
END DO
DO IP = 1,npa
IP_glob=iplg(IP)
!UL(IP) = MAX(ZERO,U(IP)-ST(IP))*DBLE(IOBPD(ITH,IP_glob))
UL(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*(1-IOBPA(IP_glob)))*DBLE(IOBPD(ITH,IP_glob))*IOBDP(IP_glob)
END DO
USTARI(1,:) = MAX(UL,U)
USTARI(2,:) = MIN(UL,U)
UIP = -THR
UIM = THR
PP = ZERO
PM = ZERO
DO IE = 1, NE
NI = INE(:,IE)
PP(NI) = PP(NI) + MAX( THR, -THETA_ACE(:,IE))
PM(NI) = PM(NI) + MIN( -THR, -THETA_ACE(:,IE))
UIP(NI) = MAX (UIP(NI), MAXVAL( USTARI(1,NI) ))
UIM(NI) = MIN (UIM(NI), MINVAL( USTARI(2,NI) ))
END DO
WII(1,:) = MIN(1.0d0,(UIP-UL) / PP)
WII(2,:) = MIN(1.0d0,(UIM-UL) / PM)
ST = ZERO
DO IE = 1, NE
DO I = 1, 3
IP = INE(I,IE)
IF (-THETA_ACE(I,IE) .GE. 0.) THEN
TMP(I) = WII(1,IP)
ELSE
TMP(I) = WII(2,IP)
END IF
END DO
BETA = MINVAL(TMP)
NI = INE(:,IE)
ST(NI) = ST(NI) + BETA * THETA_ACE(:,IE)
END DO
!
! IOBPD is the switch for removing energy coming from the shoreline
!
DO IP = 1,npa
IP_glob=iplg(IP)
U(IP) = MAX(ZERO,U(IP)-DTSI(IP)*ST(IP)*(1-IOBPA(IP_glob)))*DBLE(IOBPD(ITH,IP_glob))*IOBDP(IP_glob)
!/REF1 IF (REFPARS(3).LT.0.5.AND.IOBPD(ITH,IP_glob).EQ.0.AND.IOBPA(IP_glob).EQ.0) U(IP) = AC(IP) ! restores reflected boundary values
END DO
AC = REAL(U)
!
! 5 Update open boundaries ... this should be implemented differently ... it is better to omit any if clause in this loop ...
!
IF ( FLBPI ) THEN
RD1=RD10 - DT * REAL(ITER(IK,ITH)-IT)/REAL(ITER(IK,ITH))
RD2=RD20
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
ELSE
RD1 = 0.
RD2 = 1.
END IF
!
! NB: this treatment of the open boundary (time interpolation) is different from
! the constant boundary in the structured grids ... which restores the boundary
! to the initial value: IF ( MAPSTA(IXY).EQ.2 ) VQ(IXY) = AQ(IXY)
! Why this difference ?
!
DO IBI=1, NBI
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) ) &
/ CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
END IF
ENDDO
END IF
CALL PDLIB_exchange1DREAL(AC)
U = DBLE(AC)
END DO ! IT
! CALL EXTCDE ( 99 )
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE PDLIB_W3XYPFSFCT2
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPUG_BLOCK_IMPLICIT(FACX, FACY, DTG, VGX, VGY)
!/ ------------------------------------------------------------------- /
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Block Explicit N-Scheme
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3ODATMD, ONLY: IAPROC
USE W3GDATMD, ONLY: B_JGS_USE_JACOBI
IMPLICIT NONE
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
INTEGER DoSomething
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_USE_JACOBI=', B_JGS_USE_JACOBI
!/DEBUGSOLVER FLUSH(740+IAPROC)
DoSomething=0
IF (B_JGS_USE_JACOBI) THEN
CALL PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(FACX, FACY, DTG, VGX, VGY)
DoSomething=1
END IF
IF (DoSomething .eq. 0) THEN
WRITE(*,*) 'Error: You need to use with JGS_USE_JACOBI'
STOP 'Correct your implicit solver options'
END IF
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPUG_BLOCK_EXPLICIT(FACX, FACY, DTG, VGX, VGY)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Driver for block explicit routine
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3ODATMD, ONLY: IAPROC
USE W3GDATMD, ONLY: B_JGS_USE_JACOBI
IMPLICIT NONE
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
CALL PDLIB_EXPLICIT_BLOCK(FACX, FACY, DTG, VGX, VGY)
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
END SUBROUTINE
!/ --------------------------------------------------------------------- /
SUBROUTINE PRINT_WN_STATISTIC(string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3ODATMD, ONLY : IAPROC
USE W3GDATMD, ONLY: NK
USE W3ADATMD, ONLY: WN
USE W3GDATMD, ONLY: NSEAL
USE YOWNODEPOOL, ONLY: NP
IMPLICIT NONE
CHARACTER(*), INTENT(in) :: string
REAL TotalSumDMM, eDMM, sumDMM
INTEGER IP, IK, ISEA
WRITE(740+IAPROC,*) 'PRINT_WN_STATISTIC'
TotalSumDMM=0
DO ISEA=1,NSEAL
sumDMM=0
DO IK=0, NK
eDMM = WN(IK+1,ISEA) - WN(IK,ISEA)
sumDMM=sumDMM + abs(eDMM)
END DO
IF (ISEA .eq. 1) THEN
WRITE(740+IAPROC,*) 'ISEA=', ISEA
WRITE(740+IAPROC,*) 'sumDMM=', sumDMM
END IF
TotalSumDMM = TotalSumDMM + sumDMM
END DO
WRITE(740+IAPROC,*) 'string=', string
WRITE(740+IAPROC,*) 'TotalSumDMM=', TotalSumDMM
FLUSH(740+IAPROC)
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE WRITE_VAR_TO_TEXT_FILE(TheArr, eFile)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY : NK, NTH
USE W3WDATMD, ONLY : VA
USE W3GDATMD, ONLY : NSPEC, NX, NY, NSEAL, MAPFS
USE W3ADATMD, ONLY : MPI_COMM_WCMP
USE W3GDATMD, ONLY : GTYPE, UNGTYPE
USE W3ODATMD, ONLY : IAPROC, NAPROC, NTPROC
use yowDatapool, ONLY: rtype, istatus
USE YOWNODEPOOL, ONLY: npa, iplg, np
USE W3PARALL, ONLY: INIT_GET_ISEA
IMPLICIT NONE
INCLUDE "mpif.h"
CHARACTER(*), INTENT(in) :: eFile
REAL, INTENT(in) :: TheARR(NSPEC, npa)
!
REAL Vcoll(NSPEC,NX), VcollExp(NSPEC*NX), rVect(NSPEC*NX)
REAL CoherencyError, eVal1, eVal2, eErr
INTEGER rStatus(NX), Status(NX)
INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
INTEGER nbIncorr
INTEGER ITH, IK
INTEGER fhndl
REAL eSum
IF (IAPROC .gt. NAPROC) THEN
RETURN
END IF
IF (GTYPE .ne. UNGTYPE) THEN
RETURN
END IF
VcollExp=0
Status=0
DO IP=1,np
IP_glob=iplg(IP)
DO ISP=1,NSPEC
VcollExp(ISP+NSPEC*(IP_glob-1))=TheARR(ISP,IP)
END DO
Status(IP_glob)=1
END DO
!
! Now find global arrays
!
CoherencyError=0
IF (IAPROC .eq. 1) THEN
DO iProc=2,NAPROC
CALL MPI_RECV(rVect ,NSPEC*NX,MPI_DOUBLE , iProc-1, 37, MPI_COMM_WCMP, istatus, ierr)
CALL MPI_RECV(rStatus,NX ,MPI_INTEGER, iProc-1, 43, MPI_COMM_WCMP, istatus, ierr)
DO I=1,NX
IF (rStatus(I) .eq. 1) THEN
DO ISP=1,NSPEC
eVal1=VcollExp(ISP+NSPEC*(I-1))
eVal2=rVect(ISP+NSPEC*(I-1))
VcollExp(ISP+NSPEC*(I-1))=rVect(ISP+NSPEC*(I-1))
IF (Status(I) .eq. 1) THEN
eErr=abs(eVal1 - eVal2)
CoherencyError = CoherencyError + eErr
ELSE
VcollExp(ISP+NSPEC*(I-1))=eVal2
END IF
END DO
Status(I)=1
END IF
END DO
END DO
ELSE
CALL MPI_SEND(VcollExp,NSPEC*NX,MPI_DOUBLE , 0, 37, MPI_COMM_WCMP, ierr)
CALL MPI_SEND(Status ,NX ,MPI_INTEGER, 0, 43, MPI_COMM_WCMP, ierr)
END IF
IF (IAPROC .eq. 1) THEN
DO I=1,NX
DO ISP=1,NSPEC
Vcoll(ISP,I)=VcollExp(ISP + NSPEC*(I-1))
END DO
END DO
OPEN(fhndl, FILE=eFile)
DO IX=1,NX
eSum=sum(VColl(:,IX))
WRITE(fhndl,*) 'IX=', IX, 'eSum=', eSum
END DO
CLOSE(fhndl)
END IF
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PrintTotalOffContrib(string)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Source code for parallel debugging
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE YOWNODEPOOL, ONLY: PDLIB_CCON, NPA, PDLIB_I_DIAG, PDLIB_JA, PDLIB_IA_P
USE W3GDATMD, ONLY: NSPEC
USE W3ODATMD, ONLY : IAPROC
IMPLICIT NONE
CHARACTER(*), INTENT(in) :: string
INTEGER J, IP, JP, I, ISP
REAL TheSum1, TheSum2
J = 0
TheSum1=0
DO IP = 1, npa
DO I = 1, PDLIB_CCON(IP)
J = J + 1
IF (J .ne. PDLIB_I_DIAG(IP)) THEN
DO ISP=1,NSPEC
TheSum1=TheSum1 + abs(ASPAR_JAC(ISP,J))
END DO
END IF
END DO
END DO
!
TheSum2=0
DO IP = 1, npa
DO i = PDLIB_IA_P(IP)+1, PDLIB_IA_P(IP+1)
JP=PDLIB_JA(I)
IF (JP .ne. IP) THEN
DO ISP=1,NSPEC
TheSum2=TheSum2 + abs(ASPAR_JAC(ISP,I))
END DO
END IF
END DO
END DO
WRITE(740+IAPROC,*) 'string=', string
WRITE(740+IAPROC,*) 'TheSum12=', TheSum1, TheSum2
FLUSH(740+IAPROC)
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE COMPUTE_MEAN_PARAM (A, CG, WN, EMEAN, FMEAN, WNMEAN, AMAX)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute mean prarameter
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE CONSTANTS
USE W3GDATMD, ONLY: NK, NTH, SIG, DDEN, FTE, FTF, FTWN
!/T USE W3ODATMD, ONLY: NDST
!/S USE W3SERVMD, ONLY: STRACE
!
IMPLICIT NONE
REAL, INTENT(IN) :: A(NTH,NK), CG(NK), WN(NK)
REAL, INTENT(OUT) :: EMEAN, FMEAN, WNMEAN, AMAX
INTEGER :: IK, ITH
!/S INTEGER, SAVE :: IENT = 0
REAL :: EB(NK), EBAND
!/S CALL STRACE (IENT, 'W3SPR0')
!
EMEAN = 0.
FMEAN = 0.
WNMEAN = 0.
AMAX = 0.
!
! 1. Integral over directions
!
DO IK=1, NK
EB(IK) = 0.
DO ITH=1, NTH
EB(IK) = EB(IK) + A(ITH,IK)
AMAX = MAX ( AMAX , A(ITH,IK) )
END DO
END DO
!
! 2. Integrate over directions
!
DO IK=1, NK
EB(IK) = EB(IK) * DDEN(IK) / CG(IK)
EMEAN = EMEAN + EB(IK)
FMEAN = FMEAN + EB(IK) / SIG(IK)
WNMEAN = WNMEAN + EB(IK) / SQRT(WN(IK))
END DO
!
! 3. Add tail beyond discrete spectrum
! ( DTH * SIG absorbed in FTxx )
!
EBAND = EB(NK) / DDEN(NK)
EMEAN = EMEAN + EBAND * FTE
FMEAN = FMEAN + EBAND * FTF
WNMEAN = WNMEAN + EBAND * FTWN
!
! 4. Final processing
!
FMEAN = TPIINV * EMEAN / MAX ( 1.E-7 , FMEAN )
WNMEAN = ( EMEAN / MAX ( 1.E-7 , WNMEAN ) )**2
!
!/T WRITE (NDST,9000) EMEAN, FMEAN, WNMEAN
!
RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3SPR0 : E,F,WN MEAN ',3E10.3)
!/
!/ End of W3SPR0 ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI(DTG,FACX,FACY,VGX,VGY)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for advection part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, FACHFA, IOBPD, DMIN
USE W3GDATMD, ONLY: NSEAL, IOBDP, IOBPA, CLATS
USE W3GDATMD, ONLY: MAPSTA
USE W3WDATMD, ONLY: VA
USE W3ADATMD, ONLY: CG, DW, WN, CX, CY
!/MEMCHECK USE W3ADATMD, ONLY: MALLINFOS
USE W3IDATMD, ONLY: FLCUR
USE W3GDATMD, ONLY: ECOS, ESIN, MAPFS
USE W3PARALL, ONLY : ONESIXTH, ZERO, THR
use yowElementpool, ONLY: ne, INE
USE YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA, &
PDLIB_IE_CELL, PDLIB_POS_CELL, PDLIB_CCON, NP, NPA, &
PDLIB_IA_P, PDLIB_POSI, PDLIB_IA, PDLIB_NNZ, iplg, &
PDLIB_I_DIAG, PDLIB_JA
USE W3GDATMD, ONLY: IOBP
USE W3ODATMD, ONLY : IAPROC
USE W3PARALL, ONLY : ZERO
!/MEMCHECK USE MallocInfo_m
!/DB1 USE W3SDB1MD
!/BT1 USE W3SBT1MD
!/BT4 USE W3SBT4MD
!/BT8 USE W3SBT8MD
!/BT9 USE W3SBT9MD
!/IC1 USE W3SIC1MD
!/IC2 USE W3SIC2MD
!/IC3 USE W3SIC3MD
!/TR1 USE W3STR1MD
implicit none
REAL, INTENT(in) :: DTG, FACX, FACY, VGX, VGY
INTEGER :: IP, ISP, ISEA, IP_glob
INTEGER :: idx, IS
INTEGER :: I, J, ITH, IK, J2
INTEGER :: IE, POS, JSEA
INTEGER :: I1, I2, I3, NI(3)
INTEGER :: counter
!/REF1 INTEGER :: eIOBPDR
INTEGER :: POS_TRICK(3,2)
REAL :: DTK, TMP3
REAL :: LAMBDA(2)
REAL :: FL11, FL12
REAL :: FL21, FL22
REAL :: FL31, FL32
REAL :: CRFS(3), K(3)
REAL :: KP(3,NSPEC,NE)
REAL :: KM(3), CXY(3,2)
REAL :: K1, eSI, eVS, eVD
REAL :: eVal1, eVal2, eVal3
REAL :: DELTAL(3,NSPEC,NE)
REAL :: NM(NSPEC,NE)
REAL :: TRIA03, SIDT, CCOS, CSIN
REAL :: SPEC(NSPEC), DEPTH
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'calcARRAY_JACOBI, begin'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!!/DEBUGSRC DO JSEA=1,NSEAL
!!/DEBUGSRC WRITE(740+IAPROC,*) 'JSEA=', JSEA
!!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VS)=', minval(VSTOT(:,JSEA)), maxval(VSTOT(:,JSEA)), sum(VSTOT(:,JSEA))
!!/DEBUGSRC END DO
I = 0
IE = 0
POS = 0
I1 = 0
I2 = 0
I3 = 0
DTK = 0
TMP3 = 0
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_JACOBI SECTION 0'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
DO IE = 1, NE
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
NI = INE(:,IE)
DO IS = 1, NSPEC
ITH = 1 + MOD(IS-1,NTH)
IK = 1 + (IS-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
CXY(:,1) = CCOS * CG(IK,NI) / CLATS(NI)
CXY(:,2) = CSIN * CG(IK,NI)
IF (FLCUR) THEN
CXY(:,1) = CXY(:,1) + FACX * CX(NI)/CLATS(NI)
CXY(:,2) = CXY(:,2) + FACY * CY(NI)
ENDIF
!/MGP CXY(:,1) = CXY(:,1) - CCURX*VGX/CLATS(ISEA)
!/MGP CXY(:,2) = CXY(:,2) - CCURY*VGY
FL11 = CXY(2,1)*PDLIB_IEN(1,IE)+CXY(2,2)*PDLIB_IEN(2,IE)
FL12 = CXY(3,1)*PDLIB_IEN(1,IE)+CXY(3,2)*PDLIB_IEN(2,IE)
FL21 = CXY(3,1)*PDLIB_IEN(3,IE)+CXY(3,2)*PDLIB_IEN(4,IE)
FL22 = CXY(1,1)*PDLIB_IEN(3,IE)+CXY(1,2)*PDLIB_IEN(4,IE)
FL31 = CXY(1,1)*PDLIB_IEN(5,IE)+CXY(1,2)*PDLIB_IEN(6,IE)
FL32 = CXY(2,1)*PDLIB_IEN(5,IE)+CXY(2,2)*PDLIB_IEN(6,IE)
CRFS(1) = - ONESIXTH * (2.0d0 *FL31 + FL32 + FL21 + 2.0d0 * FL22 )
CRFS(2) = - ONESIXTH * (2.0d0 *FL32 + 2.0d0 * FL11 + FL12 + FL31 )
CRFS(3) = - ONESIXTH * (2.0d0 *FL12 + 2.0d0 * FL21 + FL22 + FL11 )
LAMBDA(1) = ONESIXTH * SUM(CXY(:,1))
LAMBDA(2) = ONESIXTH * SUM(CXY(:,2))
K(1) = LAMBDA(1) * PDLIB_IEN(1,IE) + LAMBDA(2) * PDLIB_IEN(2,IE)
K(2) = LAMBDA(1) * PDLIB_IEN(3,IE) + LAMBDA(2) * PDLIB_IEN(4,IE)
K(3) = LAMBDA(1) * PDLIB_IEN(5,IE) + LAMBDA(2) * PDLIB_IEN(6,IE)
KP(:,IS,IE) = MAX(ZERO,K(:))
DELTAL(:,IS,IE) = CRFS(:)- KP(:,IS,IE)
KM(:) = MIN(ZERO,K(:))
NM(IS,IE) = 1.d0/MIN(-THR,SUM(KM))
ENDDO
END DO
J = 0
DO IP = 1, npa
IP_glob=iplg(IP)
DO I = 1, PDLIB_CCON(IP)
J = J + 1
IE = PDLIB_IE_CELL(J)
POS = PDLIB_POS_CELL(J)
I1 = PDLIB_POSI(1,J)
I2 = PDLIB_POSI(2,J)
I3 = PDLIB_POSI(3,J)
!/DEBUGSRC WRITE(740+IAPROC,*) 'I1=', I1, ' PDLIB_I_DIAG=', PDLIB_I_DIAG(IP)
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
K1 = KP(POS,ISP,IE)
!/REF1 eIOBPDR=(1-IOBP(IP_glob))*(1-IOBPD(ITH,IP_glob))
!/REF1 IF (eIOBPDR .eq. 1) THEN
!/REF1 K1=ZERO
!/REF1 END IF
TRIA03= 1./3. * PDLIB_TRIA(IE)
DTK = K1 * DTG * IOBPD(ITH,IP_glob) * (1-IOBPA(IP_glob))
TMP3 = DTK * NM(ISP,IE)
IF (FSGEOADVECT) THEN
ASPAR_JAC(ISP,I1) = ASPAR_JAC(ISP,I1) + TRIA03 + DTK - TMP3*DELTAL(POS,ISP,IE)
ASPAR_JAC(ISP,I2) = ASPAR_JAC(ISP,I2) - TMP3*DELTAL(POS_TRICK(POS,1),ISP,IE)
ASPAR_JAC(ISP,I3) = ASPAR_JAC(ISP,I3) - TMP3*DELTAL(POS_TRICK(POS,2),ISP,IE)
ELSE
ASPAR_JAC(ISP,I1) = ASPAR_JAC(ISP,I1) + TRIA03
END IF
B_JAC(ISP,IP)=B_JAC(ISP,IP) + TRIA03 * VA(ISP,IP) * IOBDP(IP_glob) * IOBPD(ITH,IP_glob)
END DO
END DO
END DO
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_JACOBI SECTION 1'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'sum(VA)=', sum(VA)
!/DEBUGSOLVER CALL PrintTotalOffContrib("Offdiag after the geo advection")
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI2(DTG,FACX,FACY,VGX,VGY)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for advection part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, FACHFA, IOBPD, DMIN
USE W3GDATMD, ONLY: NSEAL, IOBDP, IOBPA, CLATS
USE W3GDATMD, ONLY: MAPSTA
USE W3WDATMD, ONLY: VA, VAOLD
USE W3ADATMD, ONLY: CG, DW, WN, CX, CY
!/MEMCHECK USE MallocInfo_m
!/MEMCHECK USE W3ADATMD, ONLY: MALLINFOS
USE W3IDATMD, ONLY: FLCUR
USE W3GDATMD, ONLY: ECOS, ESIN, MAPFS
USE W3PARALL, ONLY : ONESIXTH, ZERO, THR, IMEM
use yowElementpool, ONLY: ne, INE
USE YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA, &
PDLIB_IE_CELL, PDLIB_POS_CELL, PDLIB_CCON, NP, NPA, &
PDLIB_IA_P, PDLIB_POSI, PDLIB_IA, PDLIB_NNZ, iplg, &
PDLIB_I_DIAG, PDLIB_JA
USE W3GDATMD, ONLY: IOBP
USE W3ODATMD, ONLY : IAPROC
!/DB1 USE W3SDB1MD
!/BT1 USE W3SBT1MD
!/BT4 USE W3SBT4MD
!/BT8 USE W3SBT8MD
!/BT9 USE W3SBT9MD
!/IC1 USE W3SIC1MD
!/IC2 USE W3SIC2MD
!/IC3 USE W3SIC3MD
!/TR1 USE W3STR1MD
implicit none
REAL, INTENT(in) :: DTG, FACX, FACY, VGX, VGY
INTEGER :: IP, ISP, ISEA, IP_glob
INTEGER :: idx, IS
INTEGER :: I, J, ITH, IK, J2
INTEGER :: IE, POS, JSEA
INTEGER :: I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
INTEGER :: counter
!/REF1 INTEGER :: eIOBPDR
INTEGER :: POS_TRICK(3,2), IP1, IP2, IPP1, IPP2
REAL :: DTK, TMP3
REAL :: LAMBDA(2)
REAL :: FL11, FL12
REAL :: FL21, FL22
REAL :: FL31, FL32
REAL :: CRFS(3), K(3)
REAL :: KP(3)
REAL :: KM(3), CXY(3,2)
REAL :: K1, eSI, eVS, eVD
REAL :: eVal1, eVal2, eVal3
REAL :: DELTAL(3)
REAL :: NM
REAL :: IEN_LOCAL(6)
REAL :: TRIA03, SIDT, CCOS, CSIN
REAL :: SPEC(NSPEC), DEPTH
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_JACOBI SECTION 0'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
J = 0
DO IP = 1, npa
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
DO I = 1, PDLIB_CCON(IP)
J = J + 1
IE = PDLIB_IE_CELL(J)
IEN_LOCAL = PDLIB_IEN(:,IE)
POS = PDLIB_POS_CELL(J)
I1 = PDLIB_POSI(1,J)
I2 = PDLIB_POSI(2,J)
I3 = PDLIB_POSI(3,J)
IP1 = INE(POS_TRICK(POS,1),IE)
IP2 = INE(POS_TRICK(POS,2),IE)
IPP1 = POS_TRICK(POS,1)
IPP2 = POS_TRICK(POS,2)
NI = INE(:,IE)
NI_GLOB = iplg(NI)
NI_ISEA = MAPFS(1,NI_GLOB)
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
CXY(:,1) = CCOS * CG(IK,NI_ISEA) / CLATS(NI_ISEA)
CXY(:,2) = CSIN * CG(IK,NI_ISEA)
IF (FLCUR) THEN
CXY(:,1) = CXY(:,1) + FACX * CX(NI_ISEA)/CLATS(NI_ISEA)
CXY(:,2) = CXY(:,2) + FACY * CY(NI_ISEA)
ENDIF
!/MGP CXY(:,1) = CXY(:,1) - CCURX*VGX/CLATS(ISEA)
!/MGP CXY(:,2) = CXY(:,2) - CCURY*VGY
FL11 = CXY(2,1)*IEN_LOCAL(1)+CXY(2,2)*IEN_LOCAL(2)
FL12 = CXY(3,1)*IEN_LOCAL(1)+CXY(3,2)*IEN_LOCAL(2)
FL21 = CXY(3,1)*IEN_LOCAL(3)+CXY(3,2)*IEN_LOCAL(4)
FL22 = CXY(1,1)*IEN_LOCAL(3)+CXY(1,2)*IEN_LOCAL(4)
FL31 = CXY(1,1)*IEN_LOCAL(5)+CXY(1,2)*IEN_LOCAL(6)
FL32 = CXY(2,1)*IEN_LOCAL(5)+CXY(2,2)*IEN_LOCAL(6)
CRFS(1) = - ONESIXTH * (2.0d0 *FL31 + FL32 + FL21 + 2.0d0 * FL22 )
CRFS(2) = - ONESIXTH * (2.0d0 *FL32 + 2.0d0 * FL11 + FL12 + FL31 )
CRFS(3) = - ONESIXTH * (2.0d0 *FL12 + 2.0d0 * FL21 + FL22 + FL11 )
LAMBDA(1) = ONESIXTH * SUM(CXY(:,1))
LAMBDA(2) = ONESIXTH * SUM(CXY(:,2))
K(1) = LAMBDA(1) * IEN_LOCAL(1) + LAMBDA(2) * IEN_LOCAL(2)
K(2) = LAMBDA(1) * IEN_LOCAL(3) + LAMBDA(2) * IEN_LOCAL(4)
K(3) = LAMBDA(1) * IEN_LOCAL(5) + LAMBDA(2) * IEN_LOCAL(6)
KP(:) = MAX(ZERO,K(:))
DELTAL(:) = CRFS(:) - KP(:)
KM(:) = MIN(ZERO,K(:))
NM = 1.d0/MIN(-THR,SUM(KM))
K1 = KP(POS)
!/REF1 eIOBPDR=(1-IOBP(IP_glob))*(1-IOBPD(ITH,IP_glob))
!/REF1 IF (eIOBPDR .eq. 1) THEN
!/REF1 K1=ZERO
!/REF1 END IF
TRIA03 = 1./3. * PDLIB_TRIA(IE)
DTK = K1 * DTG * IOBDP(IP_glob) * IOBPD(ITH,IP_glob) * (1-IOBPA(IP_glob))
TMP3 = DTK * NM
IF (FSGEOADVECT) THEN
ASPAR_JAC(ISP,I1) = ASPAR_JAC(ISP,I1) + TRIA03 + DTK - TMP3*DELTAL(POS)
ASPAR_JAC(ISP,I2) = ASPAR_JAC(ISP,I2) - TMP3*DELTAL(IPP1)
ASPAR_JAC(ISP,I3) = ASPAR_JAC(ISP,I3) - TMP3*DELTAL(IPP2)
ELSE
ASPAR_JAC(ISP,I1) = ASPAR_JAC(ISP,I1) + TRIA03
END IF
B_JAC(ISP,IP)=B_JAC(ISP,IP) + TRIA03 * VA(ISP,IP) * IOBDP(IP_glob) * IOBPD(ITH,IP_glob)
END DO
END DO
END DO
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_JACOBI SECTION 1'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI3(IP,J,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for advection part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, FACHFA, IOBPD, DMIN
USE W3GDATMD, ONLY: NSEAL, IOBDP, IOBPA, CLATS
USE W3GDATMD, ONLY: MAPSTA
USE W3WDATMD, ONLY: VA, VAOLD
USE W3ADATMD, ONLY: CG, DW, WN, CX, CY
!/MEMCHECK USE W3ADATMD, ONLY: MALLINFOS
USE W3IDATMD, ONLY: FLCUR
USE W3GDATMD, ONLY: ECOS, ESIN, MAPFS
USE W3PARALL, ONLY : ONESIXTH, ZERO, THR, ONETHIRD
use yowElementpool, ONLY: ne, INE
USE YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA, &
PDLIB_IE_CELL, PDLIB_POS_CELL, PDLIB_CCON, NP, NPA, &
PDLIB_IA_P, PDLIB_POSI, PDLIB_IA, PDLIB_NNZ, iplg, &
PDLIB_I_DIAG, PDLIB_JA
USE W3GDATMD, ONLY: IOBP
USE W3ODATMD, ONLY : IAPROC
!/MEMCHECK USE MallocInfo_m
!/DB1 USE W3SDB1MD
!/BT1 USE W3SBT1MD
!/BT4 USE W3SBT4MD
!/BT8 USE W3SBT8MD
!/BT9 USE W3SBT9MD
!/IC1 USE W3SIC1MD
!/IC2 USE W3SIC2MD
!/IC3 USE W3SIC3MD
!/TR1 USE W3STR1MD
implicit none
INTEGER, INTENT(IN) :: IP
INTEGER, INTENT(INOUT) :: J
REAL, INTENT(in) :: DTG, FACX, FACY, VGX, VGY
REAL, INTENT(out) :: ASPAR_DIAG_LOCAL(NSPEC), B_JAC_LOCAL(NSPEC), ASPAR_OFF_DIAG_LOCAL(NSPEC)
INTEGER :: ISP, ISEA, IP_glob, IPP1, IPP2
INTEGER :: idx, IS, IP1, IP2
INTEGER :: I, ITH, IK, J2
INTEGER :: IE, POS, JSEA
INTEGER :: I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
INTEGER :: counter
!/REF1 INTEGER :: eIOBPDR
INTEGER :: POS_TRICK(3,2)
REAL*8 :: DTK, TMP3
REAL*8 :: LAMBDA(2)
REAL*8 :: FL11, FL12
REAL*8 :: FL21, FL22
REAL*8 :: FL31, FL32
REAL*8 :: CRFS(3), K(3)
REAL*8 :: KP(3)
REAL*8 :: KM(3), CXY(3,2)
REAL*8 :: K1, eSI, eVS, eVD
REAL*8 :: eVal1, eVal2, eVal3
REAL*8 :: ien_local(6)
REAL*8 :: DELTAL(3)
REAL*8 :: NM
REAL*8 :: TRIA03, SIDT, CCOS, CSIN
REAL*8 :: DEPTH
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
ASPAR_DIAG_LOCAL = 0.d0
B_JAC_LOCAL = 0.d0
ASPAR_OFF_DIAG_LOCAL = 0.d0
IP_glob=iplg(IP)
DO I = 1, PDLIB_CCON(IP)
J = J + 1
IE = PDLIB_IE_CELL(J)
IEN_LOCAL = PDLIB_IEN(:,IE)
POS = PDLIB_POS_CELL(J)
I1 = PDLIB_POSI(1,J)
I2 = PDLIB_POSI(2,J)
I3 = PDLIB_POSI(3,J)
IP1 = INE(POS_TRICK(POS,1),IE)
IP2 = INE(POS_TRICK(POS,2),IE)
IPP1 = POS_TRICK(POS,1)
IPP2 = POS_TRICK(POS,2)
NI = INE(:,IE)
NI_GLOB = iplg(NI)
NI_ISEA = MAPFS(1,NI_GLOB)
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
CXY(:,1) = CCOS * CG(IK,NI_ISEA) / CLATS(NI_ISEA)
CXY(:,2) = CSIN * CG(IK,NI_ISEA)
IF (FLCUR) THEN
CXY(:,1) = CXY(:,1) + FACX * CX(NI_ISEA)/CLATS(NI_ISEA)
CXY(:,2) = CXY(:,2) + FACY * CY(NI_ISEA)
ENDIF
!/MGP CXY(:,1) = CXY(:,1) - CCURX*VGX/CLATS(ISEA)
!/MGP CXY(:,2) = CXY(:,2) - CCURY*VGY
FL11 = CXY(2,1)*IEN_LOCAL(1)+CXY(2,2)*IEN_LOCAL(2)
FL12 = CXY(3,1)*IEN_LOCAL(1)+CXY(3,2)*IEN_LOCAL(2)
FL21 = CXY(3,1)*IEN_LOCAL(3)+CXY(3,2)*IEN_LOCAL(4)
FL22 = CXY(1,1)*IEN_LOCAL(3)+CXY(1,2)*IEN_LOCAL(4)
FL31 = CXY(1,1)*IEN_LOCAL(5)+CXY(1,2)*IEN_LOCAL(6)
FL32 = CXY(2,1)*IEN_LOCAL(5)+CXY(2,2)*IEN_LOCAL(6)
CRFS(1) = - ONESIXTH * (2.0d0 *FL31 + FL32 + FL21 + 2.0d0 * FL22 )
CRFS(2) = - ONESIXTH * (2.0d0 *FL32 + 2.0d0 * FL11 + FL12 + FL31 )
CRFS(3) = - ONESIXTH * (2.0d0 *FL12 + 2.0d0 * FL21 + FL22 + FL11 )
LAMBDA(1) = ONESIXTH * SUM(CXY(:,1))
LAMBDA(2) = ONESIXTH * SUM(CXY(:,2))
K(1) = LAMBDA(1) * IEN_LOCAL(1) + LAMBDA(2) * IEN_LOCAL(2)
K(2) = LAMBDA(1) * IEN_LOCAL(3) + LAMBDA(2) * IEN_LOCAL(4)
K(3) = LAMBDA(1) * IEN_LOCAL(5) + LAMBDA(2) * IEN_LOCAL(6)
KP(:) = MAX(ZERO,K(:))
DELTAL(:) = CRFS(:) - KP(:)
KM(:) = MIN(ZERO,K(:))
NM = 1.d0/MIN(-THR,SUM(KM))
!/REF1 eIOBPDR=(1-IOBP(IP_glob))*(1-IOBPD(ITH,IP_glob))
!/REF1 IF (eIOBPDR .eq. 1) THEN
!/REF1 K1=ZERO
!/REF1 END IF
TRIA03 = ONETHIRD * PDLIB_TRIA(IE)
DTK = KP(POS) * DBLE(DTG) * IOBDP(IP_glob) * IOBPD(ITH,IP_glob) * (1-IOBPA(IP_glob))
TMP3 = DTK * NM
! IF (IP == 224 .AND. ISP == 121) WRITE(10006,'(I10,20F20.15)') ISP, KP(POS), DTK, TMP3, DELTAL(POS)
IF (FSGEOADVECT) THEN
ASPAR_DIAG_LOCAL(ISP) = ASPAR_DIAG_LOCAL(ISP) + TRIA03 + DTK - TMP3*DELTAL(POS)
ASPAR_OFF_DIAG_LOCAL(ISP) = ASPAR_OFF_DIAG_LOCAL(ISP) - TMP3*DELTAL(IPP1)*VA(ISP,IP1)
ASPAR_OFF_DIAG_LOCAL(ISP) = ASPAR_OFF_DIAG_LOCAL(ISP) - TMP3*DELTAL(IPP2)*VA(ISP,IP2)
ELSE
ASPAR_DIAG_LOCAL(ISP) = ASPAR_DIAG_LOCAL(ISP) + TRIA03
END IF
!IF (IP == 2) WRITE(10005,'(2I10,10G20.10)') ISP, IP, VAOLD(ISP,IP)
B_JAC_LOCAL(ISP) = B_JAC_LOCAL(ISP) + TRIA03 * VAOLD(ISP,IP) * IOBDP(IP_glob) * IOBPD(ITH,IP_glob)
END DO
END DO
!/
!/ End of W3XYPFSN --------------------------------------------------- /
!/
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI4(IP,J,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for advection part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, FACHFA, IOBPD, DMIN
USE W3GDATMD, ONLY: NSEAL, IOBDP, IOBPA, CLATS
USE W3GDATMD, ONLY: MAPSTA, NK
USE W3WDATMD, ONLY: VA, VAOLD
USE W3ADATMD, ONLY: CG, DW, WN, CX, CY
!/MEMCHECK USE W3ADATMD, ONLY: MALLINFOS
USE W3IDATMD, ONLY: FLCUR
USE W3GDATMD, ONLY: ECOS, ESIN, MAPFS
USE W3PARALL, ONLY : ONESIXTH, ZERO, THR, ONETHIRD
use yowElementpool, ONLY: ne, INE
USE YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA, &
PDLIB_IE_CELL, PDLIB_POS_CELL, PDLIB_CCON, NP, NPA, &
PDLIB_IA_P, PDLIB_POSI, PDLIB_IA, PDLIB_NNZ, iplg, &
PDLIB_I_DIAG, PDLIB_JA
USE W3GDATMD, ONLY: IOBP
USE W3ODATMD, ONLY : IAPROC
!/MEMCHECK USE MallocInfo_m
!/DB1 USE W3SDB1MD
!/BT1 USE W3SBT1MD
!/BT4 USE W3SBT4MD
!/BT8 USE W3SBT8MD
!/BT9 USE W3SBT9MD
!/IC1 USE W3SIC1MD
!/IC2 USE W3SIC2MD
!/IC3 USE W3SIC3MD
!/TR1 USE W3STR1MD
implicit none
INTEGER, INTENT(IN) :: IP
INTEGER, INTENT(INOUT) :: J
REAL, INTENT(in) :: DTG, FACX, FACY, VGX, VGY
REAL, INTENT(out) :: ASPAR_DIAG_LOCAL(NSPEC), B_JAC_LOCAL(NSPEC), ASPAR_OFF_DIAG_LOCAL(NSPEC)
!
INTEGER :: IP1, IP2
INTEGER :: ITH, IK
INTEGER :: IE, POS, JSEA
INTEGER :: I, I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
INTEGER :: ISP, IP_glob, IPP1, IPP2
INTEGER :: counter
!/REF1 INTEGER :: eIOBPDR
INTEGER :: POS_TRICK(3,2)
REAL*8 :: DTK, TMP3
REAL*8 :: LAMBDA(2)
REAL*8 :: CRFS(3), K(3)
REAL*8 :: KP(3), UV_CUR(3,2)
REAL*8 :: KM(3), CSX(3), CSY(3)
REAL*8 :: K1, eSI, eVS, eVD
REAL*8 :: eVal1, eVal2, eVal3
REAL*8 :: ien_local(6)
REAL*8 :: DELTAL(3), K_X(3,NK), K_Y(3,NK), K_U(3)
REAL*8 :: CRFS_X(3,NK), CRFS_Y(3,NK), CRFS_U(3)
REAL*8 :: NM, CGFAK(3,NK)
REAL*8 :: TRIA03, SIDT, CCOS, CSIN
REAL*8 :: FL11_X, FL12_X, FL21_X, FL22_X, FL31_X, FL32_X
REAL*8 :: FL11_Y, FL12_Y, FL21_Y, FL22_Y, FL31_Y, FL32_Y
REAL*8 :: FL11_U, FL12_U, FL21_U, FL22_U, FL31_U, FL32_U
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
ASPAR_DIAG_LOCAL = ZERO
B_JAC_LOCAL = ZERO
ASPAR_OFF_DIAG_LOCAL = ZERO
IP_glob=iplg(IP)
DO I = 1, PDLIB_CCON(IP)
J = J + 1
IE = PDLIB_IE_CELL(J)
TRIA03 = ONETHIRD * PDLIB_TRIA(IE)
IEN_LOCAL = PDLIB_IEN(:,IE)
POS = PDLIB_POS_CELL(J)
IP1 = INE(POS_TRICK(POS,1),IE)
IP2 = INE(POS_TRICK(POS,2),IE)
IPP1 = POS_TRICK(POS,1)
IPP2 = POS_TRICK(POS,2)
NI = INE(:,IE)
NI_GLOB = iplg(NI)
NI_ISEA = MAPFS(1,NI_GLOB)
CRFS_U = ZERO
K_U = ZERO
IF (FLCUR) THEN
UV_CUR(:,1) = FACX * CX(NI_ISEA) / CLATS(NI_ISEA)
UV_CUR(:,2) = FACY * CY(NI_ISEA)
LAMBDA(1) = ONESIXTH*(UV_CUR(1,1)+UV_CUR(2,1)+UV_CUR(3,1))
LAMBDA(2) = ONESIXTH*(UV_CUR(1,2)+UV_CUR(2,2)+UV_CUR(3,2))
K_U(1) = LAMBDA(1) * IEN_LOCAL(1) + LAMBDA(2) * IEN_LOCAL(2)
K_U(2) = LAMBDA(1) * IEN_LOCAL(3) + LAMBDA(2) * IEN_LOCAL(4)
K_U(3) = LAMBDA(1) * IEN_LOCAL(5) + LAMBDA(2) * IEN_LOCAL(6)
FL11_U = UV_CUR(2,1)*IEN_LOCAL(1)+UV_CUR(2,2)*IEN_LOCAL(2)
FL12_U = UV_CUR(3,1)*IEN_LOCAL(1)+UV_CUR(3,2)*IEN_LOCAL(2)
FL21_U = UV_CUR(3,1)*IEN_LOCAL(3)+UV_CUR(3,2)*IEN_LOCAL(4)
FL22_U = UV_CUR(1,1)*IEN_LOCAL(3)+UV_CUR(1,2)*IEN_LOCAL(4)
FL31_U = UV_CUR(1,1)*IEN_LOCAL(5)+UV_CUR(1,2)*IEN_LOCAL(6)
FL32_U = UV_CUR(2,1)*IEN_LOCAL(5)+UV_CUR(2,2)*IEN_LOCAL(6)
CRFS_U(1) = - ONESIXTH*(2.d0 *FL31_U + FL32_U + FL21_U + 2.d0 * FL22_U)
CRFS_U(2) = - ONESIXTH*(2.d0 *FL32_U + 2.d0 * FL11_U + FL12_U + FL31_U)
CRFS_U(3) = - ONESIXTH*(2.d0 *FL12_U + 2.d0 * FL21_U + FL22_U + FL11_U)
ENDIF
DO IK = 1, NK
CSX = CG(IK,NI_ISEA) / CLATS(NI_ISEA)
CSY = CG(IK,NI_ISEA)
LAMBDA(1)=ONESIXTH * (CSX(1) + CSX(2) + CSX(3))
LAMBDA(2)=ONESIXTH * (CSY(1) + CSY(2) + CSY(3))
K_X(1,IK) = LAMBDA(1) * IEN_LOCAL(1)
K_X(2,IK) = LAMBDA(1) * IEN_LOCAL(3)
K_X(3,IK) = LAMBDA(1) * IEN_LOCAL(5)
K_Y(1,IK) = LAMBDA(2) * IEN_LOCAL(2)
K_Y(2,IK) = LAMBDA(2) * IEN_LOCAL(4)
K_Y(3,IK) = LAMBDA(2) * IEN_LOCAL(6)
FL11_X = CSX(2) * IEN_LOCAL(1)
FL12_X = CSX(3) * IEN_LOCAL(1)
FL21_X = CSX(3) * IEN_LOCAL(3)
FL22_X = CSX(1) * IEN_LOCAL(3)
FL31_X = CSX(1) * IEN_LOCAL(5)
FL32_X = CSX(2) * IEN_LOCAL(5)
FL11_Y = CSY(2) * IEN_LOCAL(2)
FL12_Y = CSY(3) * IEN_LOCAL(2)
FL21_Y = CSY(3) * IEN_LOCAL(4)
FL22_Y = CSY(1) * IEN_LOCAL(4)
FL31_Y = CSY(1) * IEN_LOCAL(6)
FL32_Y = CSY(2) * IEN_LOCAL(6)
CRFS_X(1,IK)= - ONESIXTH*(2.d0*FL31_X + FL32_X + FL21_X + 2.d0 * FL22_X)
CRFS_X(2,IK)= - ONESIXTH*(2.d0*FL32_X + 2.d0 * FL11_X + FL12_X + FL31_X)
CRFS_X(3,IK)= - ONESIXTH*(2.d0*FL12_X + 2.d0 * FL21_X + FL22_X + FL11_X)
CRFS_Y(1,IK)= - ONESIXTH*(2.d0*FL31_Y + FL32_Y + FL21_Y + 2.d0 * FL22_Y)
CRFS_Y(2,IK)= - ONESIXTH*(2.d0*FL32_Y + 2.d0 * FL11_Y + FL12_Y + FL31_Y)
CRFS_Y(3,IK)= - ONESIXTH*(2.d0*FL12_Y + 2.d0 * FL21_Y + FL22_Y + FL11_Y)
ENDDO
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
K = K_X(:,IK) * CCOS + K_Y(:,IK) * CSIN + K_U
CRFS = CRFS_X(:,IK) * CCOS + CRFS_Y(:,IK) * CSIN + CRFS_U
KM = MIN(ZERO,K)
KP = MAX(ZERO,K)
DELTAL = CRFS - KP
NM = 1.d0/MIN(-THR,SUM(KM))
DTK = KP(POS) * DBLE(DTG) * IOBDP(IP_glob) * IOBPD(ITH,IP_glob) * (1-IOBPA(IP_glob))
TMP3 = DTK * NM
!IF (IP == 224 .AND. ISP == 121) WRITE(10007,'(I10,20F20.15)') ISP, KP(POS), DTK, TMP3, DELTAL(POS)
IF (FSGEOADVECT) THEN
ASPAR_DIAG_LOCAL(ISP) = ASPAR_DIAG_LOCAL(ISP) + TRIA03 + DTK - TMP3*DELTAL(POS)
ASPAR_OFF_DIAG_LOCAL(ISP) = ASPAR_OFF_DIAG_LOCAL(ISP) - TMP3*DELTAL(IPP1)*VA(ISP,IP1)
ASPAR_OFF_DIAG_LOCAL(ISP) = ASPAR_OFF_DIAG_LOCAL(ISP) - TMP3*DELTAL(IPP2)*VA(ISP,IP2)
ELSE
ASPAR_DIAG_LOCAL(ISP) = ASPAR_DIAG_LOCAL(ISP) + TRIA03
END IF
B_JAC_LOCAL(ISP) = B_JAC_LOCAL(ISP) + TRIA03 * VAOLD(ISP,IP) * IOBDP(IP_glob) * IOBPD(ITH,IP_glob)
END DO
END DO
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI5(IE,DTG,FACX,FACY,VGX,VGY)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for advection part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, FACHFA, IOBPD, DMIN
USE W3GDATMD, ONLY: NSEAL, IOBDP, IOBPA, CLATS
USE W3GDATMD, ONLY: MAPSTA, NK
USE W3WDATMD, ONLY: VA, VAOLD
USE W3ADATMD, ONLY: CG, DW, WN, CX, CY
!/MEMCHECK USE W3ADATMD, ONLY: MALLINFOS
USE W3IDATMD, ONLY: FLCUR
USE W3GDATMD, ONLY: ECOS, ESIN, MAPFS
USE W3PARALL, ONLY : ONESIXTH, ZERO, THR, ONETHIRD
use yowElementpool, ONLY: ne, INE
USE YOWNODEPOOL, ONLY: PDLIB_IEN, PDLIB_TRIA, &
PDLIB_IE_CELL, PDLIB_POS_CELL, PDLIB_CCON, NP, NPA, &
PDLIB_IA_P, PDLIB_POSI, PDLIB_IA, PDLIB_NNZ, iplg, &
PDLIB_I_DIAG, PDLIB_JA
USE W3GDATMD, ONLY: IOBP
USE W3ODATMD, ONLY : IAPROC
!/MEMCHECK USE MallocInfo_m
!/DB1 USE W3SDB1MD
!/BT1 USE W3SBT1MD
!/BT4 USE W3SBT4MD
!/BT8 USE W3SBT8MD
!/BT9 USE W3SBT9MD
!/IC1 USE W3SIC1MD
!/IC2 USE W3SIC2MD
!/IC3 USE W3SIC3MD
!/TR1 USE W3STR1MD
implicit none
INTEGER, INTENT(IN) :: IE
REAL, INTENT(in) :: DTG, FACX, FACY, VGX, VGY
!
INTEGER :: IP, IP1, IP2
INTEGER :: ITH, IK
INTEGER :: POS, JSEA
INTEGER :: I, I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
INTEGER :: ISP, IP_glob, IPP1, IPP2
INTEGER :: counter
!/REF1 INTEGER :: eIOBPDR
INTEGER :: POS_TRICK(3,2)
REAL :: DTK(3), TMP3(NSPEC,3)
REAL :: LAMBDA(2)
REAL :: CRFS(3), K(3)
REAL :: KP(3), UV_CUR(3,2)
REAL :: KM(3), CSX(3), CSY(3)
REAL :: K1, eSI, eVS, eVD
REAL :: eVal1, eVal2, eVal3
REAL :: ien_local(6)
REAL :: DELTAL(NSPEC,3), K_X(3,NK), K_Y(3,NK), K_U(3)
REAL :: CRFS_X(3,NK), CRFS_Y(3,NK), CRFS_U(3)
REAL :: NM, CGFAK(3,NK)
REAL :: TRIA03, SIDT, CCOS, CSIN
REAL :: FL11_X, FL12_X, FL21_X, FL22_X, FL31_X, FL32_X
REAL :: FL11_Y, FL12_Y, FL21_Y, FL22_Y, FL31_Y, FL32_Y
REAL :: FL11_U, FL12_U, FL21_U, FL22_U, FL31_U, FL32_U
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
TRIA03 = ONETHIRD * PDLIB_TRIA(IE)
IEN_LOCAL = PDLIB_IEN(:,IE)
NI = INE(:,IE)
NI_GLOB = iplg(NI)
NI_ISEA = MAPFS(1,NI_GLOB)
CRFS_U = ZERO
K_U = ZERO
IF (FLCUR) THEN
UV_CUR(:,1) = CX(NI_ISEA) / CLATS(NI_ISEA)
UV_CUR(:,2) = CY(NI_ISEA)
LAMBDA(1)=ONESIXTH*(UV_CUR(1,1)+UV_CUR(2,1)+UV_CUR(3,1))
LAMBDA(2)=ONESIXTH*(UV_CUR(1,2)+UV_CUR(2,2)+UV_CUR(3,2))
K_U(1) = LAMBDA(1) * IEN_LOCAL(1) + LAMBDA(2) * IEN_LOCAL(2)
K_U(2) = LAMBDA(1) * IEN_LOCAL(3) + LAMBDA(2) * IEN_LOCAL(4)
K_U(3) = LAMBDA(1) * IEN_LOCAL(5) + LAMBDA(2) * IEN_LOCAL(6)
FL11_U = UV_CUR(2,1)*IEN_LOCAL(1)+UV_CUR(2,2)*IEN_LOCAL(2)
FL12_U = UV_CUR(3,1)*IEN_LOCAL(1)+UV_CUR(3,2)*IEN_LOCAL(2)
FL21_U = UV_CUR(3,1)*IEN_LOCAL(3)+UV_CUR(3,2)*IEN_LOCAL(4)
FL22_U = UV_CUR(1,1)*IEN_LOCAL(3)+UV_CUR(1,2)*IEN_LOCAL(4)
FL31_U = UV_CUR(1,1)*IEN_LOCAL(5)+UV_CUR(1,2)*IEN_LOCAL(6)
FL32_U = UV_CUR(2,1)*IEN_LOCAL(5)+UV_CUR(2,2)*IEN_LOCAL(6)
CRFS_U(1) = - ONESIXTH*(2.d0 *FL31_U + FL32_U + FL21_U + 2.d0 * FL22_U)
CRFS_U(2) = - ONESIXTH*(2.d0 *FL32_U + 2.d0 * FL11_U + FL12_U + FL31_U)
CRFS_U(3) = - ONESIXTH*(2.d0 *FL12_U + 2.d0 * FL21_U + FL22_U + FL11_U)
ENDIF
DO IK = 1, NK
CSX = CG(IK,NI_ISEA) / CLATS(NI_ISEA)
CSY = CG(IK,NI_ISEA)
LAMBDA(1) = ONESIXTH * (CSX(1) + CSX(2) + CSX(3))
LAMBDA(2) = ONESIXTH * (CSY(1) + CSY(2) + CSY(3))
K_X(1,IK) = LAMBDA(1) * IEN_LOCAL(1)
K_X(2,IK) = LAMBDA(1) * IEN_LOCAL(3)
K_X(3,IK) = LAMBDA(1) * IEN_LOCAL(5)
K_Y(1,IK) = LAMBDA(2) * IEN_LOCAL(2)
K_Y(2,IK) = LAMBDA(2) * IEN_LOCAL(4)
K_Y(3,IK) = LAMBDA(2) * IEN_LOCAL(6)
FL11_X = CSX(2) * IEN_LOCAL(1)
FL12_X = CSX(3) * IEN_LOCAL(1)
FL21_X = CSX(3) * IEN_LOCAL(3)
FL22_X = CSX(1) * IEN_LOCAL(3)
FL31_X = CSX(1) * IEN_LOCAL(5)
FL32_X = CSX(2) * IEN_LOCAL(5)
FL11_Y = CSY(2) * IEN_LOCAL(2)
FL12_Y = CSY(3) * IEN_LOCAL(2)
FL21_Y = CSY(3) * IEN_LOCAL(4)
FL22_Y = CSY(1) * IEN_LOCAL(4)
FL31_Y = CSY(1) * IEN_LOCAL(6)
FL32_Y = CSY(2) * IEN_LOCAL(6)
CRFS_X(1,IK) = - ONESIXTH * (2.d0*FL31_X + FL32_X + FL21_X + 2.d0 * FL22_X)
CRFS_X(2,IK) = - ONESIXTH * (2.d0*FL32_X + 2.d0 * FL11_X + FL12_X + FL31_X)
CRFS_X(3,IK) = - ONESIXTH * (2.d0*FL12_X + 2.d0 * FL21_X + FL22_X + FL11_X)
CRFS_Y(1,IK) = - ONESIXTH * (2.d0*FL31_Y + FL32_Y + FL21_Y + 2.d0 * FL22_Y)
CRFS_Y(2,IK) = - ONESIXTH * (2.d0*FL32_Y + 2.d0 * FL11_Y + FL12_Y + FL31_Y)
CRFS_Y(3,IK) = - ONESIXTH * (2.d0*FL12_Y + 2.d0 * FL21_Y + FL22_Y + FL11_Y)
ENDDO
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
K = K_X(:,IK) * CCOS + K_Y(:,IK) * CSIN + K_U
CRFS = CRFS_X(:,IK) * CCOS + CRFS_Y(:,IK) * CSIN + CRFS_U
KM = MIN(ZERO,K)
KP = MAX(ZERO,K)
DELTAL(ISP,:) = CRFS - KP
NM = 1.d0/MIN(-THR,SUM(KM))
DTK = KP * DTG * IOBDP(NI_ISEA) * IOBPD(ITH,NI_ISEA) * (1-IOBPA(NI_ISEA))
TMP3(ISP,:) = DTK * NM
ENDDO
DO I = 1, 3
IP = NI(I)
IP1 = INE(POS_TRICK(I,1),IE)
IP2 = INE(POS_TRICK(I,2),IE)
IPP1 = POS_TRICK(I,1)
IPP2 = POS_TRICK(I,2)
!ASPAR_DIAG(:,IP) = ASPAR_DIAG(:,IP) + TRIA03 + DTK(I) - TMP3(:,I) * DELTAL
!ASPAR_OFF_DIAG(:,IP1) = ASPAR_OFF_DIAG(:,IP1) - TMP3(:,IPP1) * DELTAL(:,IPP1) * VA(:,IP1)
!ASPAR_OFF_DIAG(:,IP2) = ASPAR_OFF_DIAG(:,IP2) - TMP3(:,IPP2) * DELTAL(:,IPP2) * VA(:,IP2)
ENDDO
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI_SPECTRAL(DTG,ASPAR_DIAG_LOCAL)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for spectral part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3GDATMD, ONLY: FSREFRACTION, FSFREQSHIFT
USE W3ODATMD, ONLY : IAPROC
USE YOWNODEPOOL, ONLY: np, iplg, PDLIB_SI, PDLIB_I_DIAG
USE W3GDATMD, ONLY: IOBP, MAPSTA, IOBPD, FACHFA, IOBPA, IOBDP
USE W3IDATMD, ONLY: FLLEV, FLCUR
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, MAPFS, DMIN, DSIP, NSEAL
USE W3PARALL, ONLY : PROP_REFRACTION_PR3, PROP_REFRACTION_PR1, PROP_FREQ_SHIFT, PROP_FREQ_SHIFT_M2, ZERO, IMEM
USE W3ADATMD, ONLY: CG, DW
IMPLICIT NONE
REAL, INTENT(in) :: DTG
REAL, INTENT(inout) :: ASPAR_DIAG_LOCAL(nspec,NSEAL)
INTEGER IP, IP_glob, ITH, IK
INTEGER ISEA, ISP
REAL :: eSI
REAL :: B_SIG(NSPEC), B_THE(NSPEC)
REAL :: CP_SIG(NSPEC), CM_SIG(NSPEC)
REAL :: CP_THE(NSPEC), CM_THE(NSPEC)
REAL :: CAD(NSPEC), CAS(NSPEC)
REAL :: DMM(0:NK2), eVal
REAL :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
LOGICAL :: DoLimiterRefraction = .FALSE.
LOGICAL :: DoLimiterFreqShit = .FALSE. !AR: This one is missing ...
INTEGER :: ITH0
LOGICAL :: LSIG = .FALSE.
!AR: this is missing in init ... but there is a design error in ww3_grid with FLCUR and FLLEV
LSIG = FLCUR .OR. FLLEV
DO IP = 1, np
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
eSI=PDLIB_SI(IP)
!
! The frequency shifting
!
IF (FSFREQSHIFT .AND. LSIG) THEN
IF (FreqShiftMethod .eq. 1) THEN
IF (MAPSTA(1,IP_glob) .eq. 1.and.IOBP(IP_glob).eq.1.and.IOBDP(IP_glob).eq.1.and.IOBPA(IP_glob).eq.0) THEN
CALL PROP_FREQ_SHIFT(IP, ISEA, CAS, DMM, DTG)
CP_SIG = MAX(ZERO,CAS)
CM_SIG = MIN(ZERO,CAS)
B_SIG=0
DO ITH=1,NTH
DO IK=1,NK
ISP=ITH + (IK-1)*NTH
B_SIG(ISP)= CP_SIG(ISP)/DMM(IK-1) - CM_SIG(ISP)/DMM(IK)
END DO
ISP = ITH + (NK-1)*NTH
B_SIG(ISP)= B_SIG(ISP) + CM_SIG(ISP)/DMM(NK) * FACHFA
END DO
IF (IMEM == 1) THEN
ASPAR_JAC(:,PDLIB_I_DIAG(IP))=ASPAR_JAC(:,PDLIB_I_DIAG(IP)) + B_SIG(:)*eSI
ELSE IF (IMEM == 2) THEN
ASPAR_DIAG_LOCAL(:,IP) = ASPAR_DIAG_LOCAL(:,IP) + B_SIG * eSI
ENDIF
ELSE
CAS=0
END IF
CAS_SIG(:,IP)=CAS
END IF
IF (FreqShiftMethod .eq. 2) THEN
IF (MAPSTA(1,IP_glob) .eq. 1.and.IOBP(IP_glob).eq.1) THEN
CALL PROP_FREQ_SHIFT_M2(IP, ISEA, CWNB_M2, DWNI_M2, DTG)
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'sum(CWNB_M2)=', sum(CWNB_M2)
DO ITH=1,NTH
DO IK=1,NK
ISP = ITH + (IK-1)*NTH
eVal = DWNI_M2(IK) * ( MIN(CWNB_M2(ISP - NTH), ZERO) - MAX(CWNB_M2(ISP),ZERO) )
IF (IMEM == 1) THEN
ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) - eSI * eVal
ELSE IF (IMEM == 2) THEN
ASPAR_DIAG_LOCAL(ISP,IP) = ASPAR_DIAG_LOCAL(ISP,IP) - eSI * eVal
ENDIF
END DO
eVal = DWNI_M2(NK) * MIN(CWNB_M2(ITH + (NK-1)*NTH), ZERO) * FACHFA
ITH0 = NSPEC - NTH
IF (IMEM == 1) THEN
ASPAR_JAC(ITH0 + ITH,PDLIB_I_DIAG(IP)) = ASPAR_JAC(ITH0 + ITH,PDLIB_I_DIAG(IP)) + eSI * eVal
ELSE IF (IMEM == 2) THEN
ASPAR_DIAG_LOCAL(ITH0 + ITH,IP) = ASPAR_DIAG_LOCAL(ITH0 + ITH,IP) + eSI * eVal
ENDIF
END DO
ELSE
CWNB_M2=0
END IF
CWNB_SIG_M2(:,IP)=CWNB_M2
END IF
END IF
!
! The refraction
!
IF (FSREFRACTION) THEN
!
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'refraction IP=', IP
!IF ((MAPSTA(1,IP_glob) .eq. 1).and.(SUM(IOBPD(:,IP_glob)) .EQ. NTH)) THEN
!IF (MAPSTA(1,IP_glob) .eq. 1) THEN
!IF (IOBP(IP_glob) .eq. 1) THEN
IF (MAPSTA(1,IP_glob) .eq. 1.and. IOBP(IP_glob) .eq. 1.and.IOBDP(IP_glob).eq.1.and.IOBPA(IP_glob).eq.0) THEN
!!/PR1 CALL PROP_REFRACTION_PR1(ISEA,DTG,CAD) !AR: Is this working?
!!/PR3 CALL PROP_REFRACTION_PR3(ISEA,DTG,CAD, DoLimiterRefraction)
CALL PROP_REFRACTION_PR3(IP,ISEA,DTG,CAD,DoLimiterRefraction)
ELSE
CAD=ZERO
END IF
!/DEBUGREFRACTION WRITE(740+IAPROC,*) 'refraction IP=', IP, ' ISEA=', ISEA
!/DEBUGREFRACTION WRITE(740+IAPROC,*) 'sum(abs(CAD))=', sum(abs(CAD))
CAD_THE(:,IP)=CAD
CP_THE = DTG*MAX(ZERO,CAD)
CM_THE = DTG*MIN(ZERO,CAD)
B_THE(:) = CP_THE(:) - CM_THE(:)
IF (IMEM == 1) THEN
ASPAR_JAC(:,PDLIB_I_DIAG(IP))=ASPAR_JAC(:,PDLIB_I_DIAG(IP)) + B_THE(:)*eSI
ELSEIF (IMEM ==2) THEN
ASPAR_DIAG_LOCAL(:,IP) = ASPAR_DIAG_LOCAL(:,IP) + B_THE(:)*eSI
ENDIF
END IF
END DO
!!/DEBUGSOLVER CALL PrintTotalOffContrib("Offdiag after the refraction")
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE CALCARRAY_JACOBI_SOURCE(DTG,ASPAR_DIAG_LOCAL)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for source part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3ODATMD, ONLY : IAPROC
USE YOWNODEPOOL, ONLY: iplg, PDLIB_SI, PDLIB_I_DIAG, NPA, NP
USE W3ADATMD, ONLY: CG, DW, WN
USE W3WDATMD, ONLY: UST, USTDIR
USE W3GDATMD, ONLY: NK, NTH, NSPEC, MAPFS, optionCall, DMIN
USE W3GDATMD, ONLY: IOBP, MAPSTA, FACP, SIG, IOBPD, IOBPA, IOBDP
USE W3PARALL, ONLY: IMEM
USE W3GDATMD, ONLY: NSEAL, CLATS
!/DB1 USE W3SDB1MD
USE W3WDATMD, ONLY: VA, VSTOT, VDTOT, SHAVETOT
USE constants, ONLY : TPI, TPIINV, GRAV
IMPLICIT NONE
REAL, INTENT(in) :: DTG
REAL, INTENT(inout) :: ASPAR_DIAG_LOCAL(:,:)!NSPEC,NSEAL)
REAL, PARAMETER :: COEF4 = 5.0E-07
REAL, PARAMETER :: FACDAM = 1
INTEGER JSEA, IP, IP_glob, ISEA
INTEGER IK, ITH, ISP, IS0
LOGICAL :: LBREAK
REAL :: eSI, eVS, eVD, SIDT
REAL :: DEPTH, DAM(NSPEC), RATIO, MAXDAC, VSDB(NSPEC), VDDB(NSPEC)
REAL :: PreVS, eDam, DVS, FREQ, EMEAN, FMEAN, WNMEAN, AMAX, CG1(NK),WN1(NK),SPEC_VA(NSPEC)
REAL TheFactor
DO JSEA = 1, NP
IP = JSEA
IP_glob = iplg(IP)
ISEA = MAPFS(1,IP_glob)
IF (IOBP(IP_glob).eq.1..and.IOBDP(IP_glob).eq.1.and.IOBPA(IP_glob).eq.0) THEN
DO IK=1, NK
DAM(1+(IK-1)*NTH) = FACP / ( SIG(IK) * WN(IK,ISEA)**3 )
END DO
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
DAM(ITH+IS0) = DAM(1+IS0)
END DO
END DO
eSI = PDLIB_SI(IP)
SIDT = eSI * DTG
DEPTH = DW(ISEA)
!/DB1 VSDB = 0.
!/DB1 VDDB = 0.
!/DB1 CG1 = CG(1:NK,ISEA)
!/DB1 WN1 = WN(1:NK,ISEA)
!/DB1 DO IK=1,NK
!/DB1 DO ITH=1,NTH
!/DB1 ISP=ITH + (IK-1)*NTH
!/DB1 SPEC_VA(ISP) = VA(ISP,JSEA) * CG(IK,ISEA) / CLATS(ISEA)
!/DB1 ENDDO
!/DB1 ENDDO
!/DB1 CALL COMPUTE_MEAN_PARAM(SPEC_VA, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX)
!/DB1 CALL W3SDB1 ( JSEA, SPEC_VA, DEPTH, EMEAN, FMEAN, WNMEAN, CG1, LBREAK, VSDB, VDDB )
! IF (JSEA == 10000) WRITE(*,'(2I20,10F20.10)') JSEA, ISEA, SUM(VSTOT(:,JSEA)), SUM(VDTOT(:,JSEA)), SUM(VSDB),SUM(VDDB), DEPTH, EMEAN, FMEAN, WNMEAN
DO IK=1,NK
DO ITH=1,NTH
ISP=ITH + (IK-1)*NTH
IF (SHAVETOT(JSEA)) THEN ! Limit ONLY the source term part ...
MAXDAC = FACDAM * DAM(ISP)
TheFactor = DTG / MAX ( 1. , (1.-DTG*VDTOT(ISP,JSEA)))
DVS = VSTOT(ISP,JSEA) * TheFactor
DVS = SIGN(MIN(MAXDAC,ABS(DVS)),DVS)
PreVS = DVS / TheFactor
ELSE
PreVS = VSTOT(ISP,JSEA)
END IF
eVS = PreVS / CG(IK,ISEA) * CLATS(ISEA)
eVD = DBLE(VDTOT(ISP,JSEA))
!/DB1 eVS = eVS + DBLE(VSDB(ISP)) / CG(IK,ISEA) * CLATS(ISEA)
!/DB1 eVD = evD + DBLE(VDDB(ISP))
B_JAC(ISP,IP) = B_JAC(ISP,IP) + SIDT * (eVS - eVD*VA(ISP,JSEA))
IF (IMEM == 1) THEN
ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) - SIDT * eVD
ELSE IF (IMEM == 2) THEN
ASPAR_DIAG_LOCAL(ISP,IP) = ASPAR_DIAG_LOCAL(ISP,IP) - SIDT * eVD
ENDIF
END DO
END DO
END IF
END DO
!!/DEBUGSOLVER CALL PrintTotalOffContrib("Offdiag after the source terms") !AR: Need to rewrite for IMEM == 2
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Before frequency shifting business'
!!/DEBUGSOLVER FLUSH(740+IAPROC)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI_SOURCE_LOCAL(DTG,ASPAR_DIAG_SOURCE, B_JAC_SOURCE)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Compute matrix coefficients for source part
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE W3ODATMD, ONLY : IAPROC
USE YOWNODEPOOL, ONLY: iplg, PDLIB_SI, PDLIB_I_DIAG
USE W3ADATMD, ONLY: CG, DW, WN
USE W3WDATMD, ONLY: UST, USTDIR
USE W3GDATMD, ONLY: NK, NTH, NSPEC, MAPFS, optionCall, DMIN
USE W3GDATMD, ONLY: IOBP, MAPSTA, FACP, SIG, IOBPD, IOBPA, IOBDP
USE W3GDATMD, ONLY: NSEAL, CLATS
USE W3WDATMD, ONLY: VA, VSTOT, VDTOT, SHAVETOT
!/DB1 USE W3SDB1MD
USE constants, ONLY : TPI, TPIINV, GRAV
IMPLICIT NONE
REAL, INTENT(in) :: DTG
REAL, INTENT(out) :: ASPAR_DIAG_SOURCE(NSPEC), B_JAC_SOURCE(NSPEC)
REAL, PARAMETER :: COEF4 = 5.0E-07
REAL, PARAMETER :: FACDAM = 1
INTEGER JSEA, IP, IP_glob, ISEA
LOGICAL :: LBREAK
INTEGER IK, ITH, ISP, IS0
REAL :: eSI, eVS, eVD, SIDT
REAL :: DEPTH, DAM(NSPEC), RATIO, MAXDAC, VSDB(NSPEC), VDDB(NSPEC)
REAL :: PreVS, eDam, DVS, FREQ, EMEAN, FMEAN, WNMEAN, AMAX, CG1(NSPEC),WN1(NSPEC),SPEC_VA(NSPEC)
REAL TheFactor
DO JSEA=1,NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA=MAPFS(1,IP_glob)
IF (IOBP(IP_glob).eq.1.and.MAPSTA(1,IP_glob).eq.1.and.IOBDP(IP_glob).eq.1.and.IOBPA(IP_glob).eq.0) THEN
DO IK=1, NK
DAM(1+(IK-1)*NTH) = FACP / ( SIG(IK) * WN(IK,ISEA)**3 )
END DO
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
DAM(ITH+IS0) = DAM(1+IS0)
END DO
END DO
eSI = PDLIB_SI(IP)
SIDT = eSI * DTG
DEPTH = DW(ISEA)
!/DB1 CG1 = CG(1:NK,ISEA)
!/DB1 WN1 = WN(1:NK,ISEA)
!/DB1 DO IK=1,NK
!/DB1 DO ITH=1,NTH
!/DB1 ISP=ITH + (IK-1)*NTH
!/DB1 SPEC_VA(ISP) = VA(ISP,JSEA) * CG(IK,ISEA) / CLATS(ISEA)
!/DB1 ENDDO
!/DB1 ENDDO
!/DB1 SPEC_VA = VA(:,JSEA)
!/DB1 CALL COMPUTE_MEAN_PARAM(SPEC_VA, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX)
!/DB1 CALL W3SDB1 ( JSEA, SPEC_VA, DEPTH, EMEAN, FMEAN, WNMEAN, CG1, LBREAK, VSDB, VDDB )
DO IK=1,NK
DO ITH=1,NTH
ISP=ITH + (IK-1)*NTH
IF (SHAVETOT(JSEA)) THEN ! Limit ONLY the source term part ...
MAXDAC = FACDAM * DAM(ISP)
TheFactor = DTG / MAX ( 1. , (1.-DTG*VDTOT(ISP,JSEA)))
DVS = VSTOT(ISP,JSEA) * TheFactor
DVS = SIGN(MIN(MAXDAC,ABS(DVS)),DVS)
PreVS = DVS / TheFactor
ELSE
PreVS=VSTOT(ISP,JSEA)
END IF
eVS = DBLE(PreVS) / CG(IK,ISEA) * CLATS(ISEA)
eVD = DBLE(VDTOT(ISP,JSEA))
!/DB1 eVS = eVS + DBLE(VSDB(ISP)) / CG(IK,ISEA) * CLATS(ISEA)
!/DB1 eVD = eVD + DBLE(VDDB(ISP))
B_JAC_SOURCE(ISP) = B_JAC_SOURCE(ISP) + SIDT * (eVS - eVD*VA(ISP,IP))
ASPAR_DIAG_SOURCE(ISP) = ASPAR_DIAG_SOURCE(ISP) - SIDT * eVD
END DO
END DO
END IF
END DO
!!/DEBUGSOLVER CALL PrintTotalOffContrib("Offdiag after the source terms")
!!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Before frequency shifting business'
!!/DEBUGSOLVER FLUSH(740+IAPROC)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE ADD_SOURCE_TERMS_NONLINEAR(DTG)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Add source terms nonlinera-
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3ODATMD, ONLY : IAPROC
USE YOWNODEPOOL, ONLY: iplg, PDLIB_SI, PDLIB_I_DIAG
USE W3ADATMD, ONLY: CG, DW, WN, BEDFORMS, TAUBBL
USE W3GDATMD, ONLY: NK, NTH, NSPEC, MAPFS, optionCall, DMIN
USE W3GDATMD, ONLY: IOBP, MAPSTA, IOBDP, IOBPA
!/BT4 USE W3GDATMD, ONLY: SED_D50, SED_PSIC
USE W3GDATMD, ONLY: NSEAL, CLATS
USE W3WDATMD, ONLY: VA, VSTOT, VDTOT, SHAVETOT
!/DB1 USE W3SDB1MD
!/DBX USE W3SDBXMD
!/TR1 USE W3STR1MD
!/TRX USE W3STRXMD
!/BT1 USE W3SBT1MD
!/BT4 USE W3SBT4MD
!/BT8 USE W3SBT8MD
!/BT9 USE W3SBT9MD
!/BTX USE W3SBTXMD
!/BS1 USE W3SBS1MD
!/BSX USE W3SBSXMD
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
REAL, INTENT(in) :: DTG
INTEGER IP, IP_glob, ISEA, IX, IY, JSEA
REAL :: SPEC_VA(NSPEC)
REAL :: CG1(NK), WN1(NK)
REAL :: eSI, eVS, eVD, SIDT
REAL :: DEPTH
INTEGER :: ITH, IK, ISP
REAL :: PreVS, AMAX, EMEAN, FMEAN, WNMEAN, D50, PSIC, TMP1(2), TMP2(3)
LOGICAL :: LBREAK
!/DB1 REAL :: VSDB(NSPEC), VDDB(NSPEC)
!/TR1 REAL :: VSTR(NSPEC), VDTR(NSPEC)
!/BT1 REAL :: VSBT(NSPEC), VDBT(NSPEC)
!/BT4 REAL :: VSBT(NSPEC), VDBT(NSPEC)
!/BS1 REAL :: VSBS(NSPEC), VDBS(NSPEC)
!/S CALL STRACE (IENT, 'ADD_SOURCE_TERMS_NONLINEAR')
DO JSEA=1,NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA=MAPFS(1,IP_glob)
eSI=PDLIB_SI(IP)
SIDT = eSI * DTG
DEPTH = DW(ISEA)
CG1 = CG(1:NK,ISEA)
WN1 = WN(1:NK,ISEA)
SPEC_VA = VA(:,JSEA)
CALL COMPUTE_MEAN_PARAM(SPEC_VA, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX)
!/DB1 VSDB = 0.
!/TR1 VSTR = 0.
!/BT1 VSBT = 0.
!/DB1 VDDB = 0.
!/TR1 VDTR = 0.
!/BT1 VDBT = 0.
!/TR1 CALL W3STR1 ( SPEC_VA, CG1, WN1, DEPTH, IX,VSTR, VDTR )
!/TRX CALL W3STRX
!/DB1 CALL W3SDB1 ( JSEA, SPEC_VA, DEPTH, EMEAN, FMEAN, WNMEAN, CG1, LBREAK, VSDB, VDDB )
!/BT1 CALL W3SBT1 ( SPEC_VA, CG1, WN1, DEPTH, VSBT, VDBT )
!/BT4 D50=SED_D50(ISEA)
!/BT4 PSIC=SED_PSIC(ISEA)
!/BT4 TMP1=TAUBBL(JSEA,1:2)
!/BT4 TMP2=BEDFORMS(JSEA,1:3)
!/BT4 CALL W3SBT4 ( SPEC_VA, CG1, WN1, DEPTH, D50, PSIC, TMP1,TMP2, VSBT, VDBT, IX, IY )
!/BT8 CALL W3SBT8 ( SPEC_VA, DEPTH, VSBT, VDBT, IX, IY )
!/BT9 CALL W3SBT9 ( SPEC_VA, DEPTH, VSBT, VDBT, IX, IY )
!/BTX CALL W3SBTX
!
!/BS1 CALL W3SBS1 ( SPEC_VA, CG1, WN1, DEPTH, CX, CY,TAUSCX, TAUSCY, VSBS, VDBS )
!/BSX CALL W3SBSX
DO IK=1,NK
DO ITH=1,NTH
ISP=ITH + (IK-1)*NTH
PreVS=0
eVD=0
!/DB1 PreVS = PreVS + VSDB(ISP)
!/TR1 PreVS = PreVS + VSTR(ISP)
!/BT1 PreVS = PreVS + VSBT(ISP)
!/BS1 PreVS = PreVS + VSBS(ISP)
eVS=DBLE(PreVS) / CG(IK,ISEA) * CLATS(ISEA)
!/DB1 eVD=eVD+DBLE(MIN(0., VDDB(ISP)))
!/TR1 eVD=eVD+DBLE(MIN(0., VDTR(ISP)))
!/BT1 eVD=eVD+DBLE(MIN(0., VDBT(ISP)))
!/BS1 eVD=eVD+DBLE(MIN(0., VDBS(ISP)))
IF (optionCall .eq. 1) THEN
B_JAC(ISP,IP) = B_JAC(ISP,IP) + SIDT * eVS
ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) - SIDT * eVD
ELSE IF (optionCall .eq. 2) THEN
B_JAC(ISP,IP) = B_JAC(ISP,IP) + SIDT * (eVS - eVD*VA(ISP,IP))
ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) - SIDT * eVD
ELSE IF (optionCall .eq. 3) THEN ! All source terms go with the REAL +- sign. E.g. dissipation is negative
B_JAC(ISP,IP) = B_JAC(ISP,IP) + SIDT * (eVS - eVD*VA(ISP,IP))
ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(IP)) - SIDT * eVD
END IF
ENDDO
ENDDO
ENDDO
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE APPLY_BOUNDARY_CONDITION_VA
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Boudary conditions on VA
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE yowRankModule, ONLY : IPGL_npa
USE W3GDATMD, ONLY: NSEAL, CLATS, GTYPE, UNGTYPE
USE W3WDATMD, ONLY: TIME
USE W3TIMEMD, ONLY: DSEC21
USE W3ADATMD, ONLY: CG, CX, CY
USE W3WDATMD, ONLY: VA
USE W3GDATMD, ONLY: NK, NK2, NTH, ECOS, ESIN, NSPEC
USE W3ODATMD, ONLY: TBPI0, TBPIN, FLBPI, IAPROC, NAPROC, BBPI0, BBPIN, ISBPI, NBI
USE W3PARALL, ONLY : ISEA_TO_JSEA
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
REAL :: RD1, RD2, RD10, RD20
REAL :: eVA, eAC
INTEGER :: IK, ITH, ISEA, JSEA
INTEGER :: IBI, ISP
!/S CALL STRACE (IENT, 'APPLY_BOUNDARY_CONDITION_VA')
IF (GTYPE .eq. UNGTYPE) THEN
IF ( FLBPI ) THEN
RD10 = DSEC21 ( TBPI0, TIME )
RD20 = DSEC21 ( TBPI0, TBPIN )
ELSE
RD10=1.
RD20=0.
END IF
IF (FLBPI .and. (IAPROC .le. NAPROC)) THEN
RD1=RD10 ! I am not completely sure about that
RD2=RD20
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
ELSE
RD1 = 0.
RD2 = 1.
END IF
DO IBI=1, NBI
ISEA=ISBPI(IBI)
JSEA=ISEA_TO_JSEA(ISEA)
IF (JSEA .gt. 0) THEN
DO ITH=1,NTH
DO IK=1,NK
ISP=ITH + (IK-1)*NTH
eAC = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) ) &
/ CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
eVA = MAX(0., CG(IK,ISEA)/CLATS(ISEA)*eAC)
VA(ISP,JSEA) = eVA
END DO
END DO
END IF
END DO
END IF
END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE APPLY_BOUNDARY_CONDITION
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Apply boundary conditions
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE YOWNODEPOOL, ONLY: npa
USE yowRankModule, ONLY : IPGL_npa
USE W3GDATMD, ONLY: NSEAL, CLATS, MAPSF
USE W3WDATMD, ONLY: TIME
USE W3TIMEMD, ONLY: DSEC21
USE W3WDATMD, ONLY : VA
USE W3ADATMD, ONLY: CG, CX, CY
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC
USE W3ODATMD, ONLY: TBPI0, TBPIN, FLBPI, IAPROC, BBPI0, BBPIN, ISBPI, NBI
!/DEBUGIOBC USE W3GDATMD, ONLY: DDEN
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/DEBUGSOLVER REAL*8 :: sumAC(NSPEC)
!/DEBUGSOLVER REAL :: sumBPI0(NSPEC), sumBPIN(NSPEC), sumCG, sumCLATS
!/DEBUGIOBC REAL :: ETOT, HSIG_bound, eVA, eAC, FACTOR
REAL :: RD1, RD2, RD10, RD20
INTEGER :: IK, ITH, ISEA
INTEGER :: IBI, IP_glob, ISP, JX
!/S CALL STRACE (IENT, 'APPLY_BOUNDARY_CONDITION')
IF ( FLBPI ) THEN
RD10 = DSEC21 ( TBPI0, TIME )
RD20 = DSEC21 ( TBPI0, TBPIN )
ELSE
RD10=1.
RD20=0.
END IF
IF ( FLBPI ) THEN
RD1=RD10 ! I am not completely sure about that
RD2=RD20
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
ELSE
RD1 = 0.
RD2 = 1.
END IF
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Begin of APPLY_BOUNDARY_CONDITION'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NBI=', NBI
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVER sumAC=0
!/DEBUGSOLVER sumBPI0=0
!/DEBUGSOLVER sumBPIN=0
!/DEBUGSOLVER sumCG=0
!/DEBUGSOLVER sumCLATS=0
DO IBI=1, NBI
ISEA=ISBPI(IBI)
IP_glob = MAPSF(ISEA,1)
JX=IPGL_npa(IP_glob)
IF (JX .gt. 0) THEN
DO ITH=1,NTH
DO IK=1,NK
ISP=ITH + (IK-1)*NTH
VA(ISP,JX) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) ) &
/ CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
END DO
END DO
!/DEBUGIOBC ETOT=0
!/DEBUGIOBC DO ITH=1,NTH
!/DEBUGIOBC DO IK=1,NK
!/DEBUGIOBC FACTOR = DDEN(IK)/CG(IK,ISEA)
!/DEBUGIOBC ISP=ITH + (IK-1)*NTH
!/DEBUGIOBC eAC=REAL(VA(ISP,JX))
!/DEBUGIOBC eVA=CG(IK,ISEA)/CLATS(ISEA)*eAC
!/DEBUGIOBC ETOT = ETOT + eVA*FACTOR
!/DEBUGIOBC END DO
!/DEBUGIOBC END DO
!/DEBUGIOBC HSIG_bound=4.*SQRT(ETOT)
!/DEBUGIOBC WRITE(740+IAPROC,*) 'IBI=', IBI, ' HSIG=', HSIG_bound
!/DEBUGSOLVER sumAC=sumAC + VA(:,JX)
!/DEBUGSOLVER sumBPI0=sumBPI0 + BBPI0(:,IBI)
!/DEBUGSOLVER sumBPIN=sumBPIN + BBPIN(:,IBI)
!/DEBUGSOLVER sumCG=sumCG + CG(IK,ISBPI(IBI))
!/DEBUGSOLVER sumCLATS=sumCLATS + CLATS(ISBPI(IBI))
END IF
ENDDO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'RD1=', RD1, ' RD2=', RD2
!/DEBUGSOLVER DO ISP=1,NSPEC
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'RD1=', RD1, ' RD2=', RD2
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumAC=', sumAC(ISP)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumBPI0=', sumBPI0(ISP)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumBPIN=', sumBPIN(ISP)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumCG=', sumCG
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'ISP=', ISP, 'sumCLATS=', sumCLATS
!/DEBUGSOLVER END DO
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Begin of APPLY_BOUNDARY_CONDITION'
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(npa) after boundary", npa)
END IF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE ACTION_LIMITER_LOCAL(IP,ACLOC,ACOLD, DTG)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Computation of the limiter function
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
use YOWNODEPOOL, ONLY: iplg
USE CONSTANTS, ONLY : GRAV, TPI
USE W3ADATMD, ONLY : WN, CG
USE W3GDATMD, ONLY : NTH, NK, NSPEC, MAPFS, SIG, FACP
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER, INTENT(in) :: IP
REAL, INTENT(in) :: ACOLD(NSPEC)
REAL, INTENT(inout) :: ACLOC(NSPEC)
REAL, INTENT(in) :: DTG
INTEGER :: MELIM = 1
REAL :: LIMFAK = 0.1
REAL :: CONST, SND, eWN, eWK, eWKpow
REAL :: eFact, eSPSIG
REAL :: NewVAL
REAL :: OLDAC, NEWAC, NEWDAC
REAL :: MAXDAC
REAL :: dac, limac, eDam
INTEGER IP_glob, ISEA
INTEGER :: IK, ITH, ISP
LOGICAL :: LLIMITER_WWM
!/S CALL STRACE (IENT, 'ACTION_LIMITER_LOCAL')!
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
eSPSIG=SIG(NK)
CONST = TPI**2*3.0*1.0E-7*DTG*eSPSIG
SND = TPI*5.6*1.0E-3
LLIMITER_WWM = .false.
IF (LLIMITER_WWM) THEN
MAXDAC = 0
DO IK=1,NK
IF (MELIM .eq. 1) THEN
eFact=2.*SIG(IK)
eWN=WN(IK,ISEA)
eWK=eWN
eWKpow=eWK**3
MAXDAC = DBLE(0.0081*LIMFAK/(eFact*eWKpow*CG(IK,ISEA)))
END IF
DO ITH=1,NTH
ISP=ITH + (IK-1)*NTH
NEWAC = ACLOC(ISP)
OLDAC = ACOLD(ISP)
NEWDAC = NEWAC - OLDAC
NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC)
NewVAL = MAX(0., OLDAC + NEWDAC )
ACLOC(ISP) = NewVAL
END DO
END DO
ELSE
DO IK = 1, NK
eDam=DBLE(FACP / (SIG(IK) * WN(IK,ISEA)**3))
DO ITH=1,NTH
isp = ITH + (IK-1)*NTH
dac = ACLOC(isp) - ACOLD(isp)
limac = SIGN (MIN(eDam,ABS(dac)),dac)
ACLOC(isp) = MAX(0., ACLOC(isp) + limac)
END DO
END DO
ENDIF
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(FACX, FACY, DTG, VGX, VGY)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Block Gauss Seidel and Jacobi solver
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!
!/S USE W3SERVMD, ONLY: STRACE
!/
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_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
!/DEBUGSRC USE W3GDATMD, ONLY: optionCall
!/DEBUGSRC USE W3WDATMD, ONLY: SHAVETOT
USE YOWNODEPOOL, ONLY: PDLIB_I_DIAG, PDLIB_IA_P, PDLIB_JA, np
USE YOWNODEPOOL, ONLY: PDLIB_SI, PDLIB_NNZ, PDLIB_CCON
use yowDatapool, ONLY: rtype
use YOWNODEPOOL, ONLY: npa, iplg
use yowExchangeModule, ONLY : PDLIB_exchange2DREAL
USE W3ADATMD, ONLY: WN
USE MPI, ONLY : MPI_SUM, MPI_INT
USE W3ADATMD, ONLY: MPI_COMM_WCMP
!/MEMCHECK USE W3ADATMD, ONLY: MALLINFOS
USE W3GDATMD, ONLY: IOBP, IOBPD, NSEA, SIG, IOBDP, IOBPA
USE W3GDATMD, ONLY: NK, NK2, NTH, ECOS, ESIN, NSPEC, MAPFS
USE W3WDATMD, ONLY: TIME
USE W3ODATMD, ONLY: NBI
USE W3TIMEMD, ONLY: DSEC21
USE W3GDATMD, ONLY: NSEAL, CLATS, FACHFA
USE W3IDATMD, ONLY: FLCUR
USE W3WDATMD, ONLY: VA, VAOLD, VSTOT, VDTOT
USE W3ADATMD, ONLY: CG, CX, CY
USE W3ODATMD, ONLY: TBPIN, FLBPI, IAPROC
USE W3PARALL, ONLY : IMEM
USE W3PARALL, ONLY : INIT_GET_JSEA_ISPROC, ZERO, THR8
USE W3PARALL, ONLY : ListISPprevDir, ListISPnextDir
USE W3PARALL, ONLY : JX_TO_JSEA
USE W3GDATMD, ONLY: B_JGS_NLEVEL, B_JGS_SOURCE_NONLINEAR
USE yowfunction, ONLY : pdlib_abort
!/MEMCHECK USE MallocInfo_m
implicit none
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
!
INTEGER :: IP, ISP, ITH, IK, JSEA, ISEA, IP_glob
INTEGER :: myrank
INTEGER :: nbIter, ISPnextDir, ISPprevDir
INTEGER :: ISPp1, ISPm1, JP, ICOUNT1, ICOUNT2
! for the exchange
REAL :: CCOS, CSIN, CCURX, CCURY
REAL :: eSum(NSPEC)
REAL :: eA_THE, eC_THE, eA_SIG, eC_SIG
REAL :: CAD(NSPEC), CAS(NSPEC), ACLOC(NSPEC)
REAL :: CP_SIG(NSPEC), CM_SIG(NSPEC)
REAL :: eFactM1, eFactP1
REAL :: Sum_Prev, Sum_New
REAL :: prop_conv, eSI, p_is_converged
REAL :: Sum_L2, Sum_L2_GL
REAL :: DMM(0:NK2)
REAL :: DiffNew
REAL :: eDiff(NSPEC), eProd(NSPEC)
REAL :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
REAL :: VAnew(NSPEC), VFLWN(1-NTH:NSPEC)
REAL :: VAAnew(1-NTH:NSPEC+NTH), VAAacloc(1-NTH:NSPEC+NTH)
REAL :: VAinput(NSPEC), VAacloc(NSPEC), eDiffB(NSPEC),ASPAR_DIAG(NSPEC)
REAL :: aspar_diag_local(nspec), aspar_off_diag_local(nspec), b_jac_local(nspec)
REAL :: eDiffSing, eSumPart
REAL :: eVal1, eVal2, extmp(nspec,nseal)
REAL :: eVA
LOGICAL :: LCONVERGED(NSEAL)
!/DEBUGSRC REAL :: IntDiff, eVA_w3srce, eVAsolve, SumACout
!/DEBUGSRC REAL :: SumVAin, SumVAout, SumVAw3srce, SumVS, SumVD, VS_w3srce
!/DEBUGSRC REAL :: VAsolve(NSPEC)
!/DEBUGSRC REAL*8 :: ACsolve
!/DEBUGSRC REAL :: eB
REAL :: ASPAR_DIAG_ALL(NSPEC,NSEAL)
!/DEBUGSOLVERCOH REAL :: TheARR(NSPEC, npa)
!/DEBUGSOLVERCOH REAL :: PRE_VA(NSPEC, npa)
!/DEBUGSOLVERCOH REAL :: OffDIAG(NSPEC, npa)
!/DEBUGSOLVERCOH REAL*8 :: eOff(NSPEC)
!/DEBUGSOLVERCOH REAL*8 :: eSum1(NSPEC), eSum2(NSPEC)
CHARACTER(len=128) eFile
INTEGER ierr, i
INTEGER JP_glob
INTEGER is_converged, itmp
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 0'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
CCURX = FACX
CCURY = FACY
CALL MPI_COMM_RANK(MPI_COMM_WCMP, myrank, ierr)
!
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK, begin'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NX=', NX
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NP=', NP
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NPA=', NPA
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NSEA=', NSEA
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NSEAL=', NSEAL
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NBI=', NBI
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_TERMINATE_NORM=', B_JGS_TERMINATE_NORM
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_TERMINATE_DIFFERENCE=', B_JGS_TERMINATE_DIFFERENCE
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_TERMINATE_MAXITER=', B_JGS_TERMINATE_MAXITER
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_MAXITER=', B_JGS_MAXITER
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_BLOCK_GAUSS_SEIDEL=', B_JGS_BLOCK_GAUSS_SEIDEL
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FSREFRACTION=', FSREFRACTION
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FSFREQSHIFT=', FSFREQSHIFT
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_LIMITER=', B_JGS_LIMITER
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_BLOCK_GAUSS_SEIDEL=', B_JGS_BLOCK_GAUSS_SEIDEL
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSRC WRITE(740+IAPROC,*) 'optionCall=', optionCall
!/DEBUGSRC FLUSH(740+IAPROC)
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 1'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
! 1.b Initialize arrays
!
CALL SETDEPTH_PDLIB
!
! 2. Convert to Wave Action ---------------- *
!
!/DEBUGSRC WRITE(740+IAPROC,*) 'NSEAL =', NSEAL, 'NP =', NP, 'NPA =', NPA
DO JSEA=1,NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA=MAPFS(1,IP_glob)
!!/DEBUGSRC WRITE(740+IAPROC,*) 'IP =', IP, 'IP_glob =', IP_glob, 'ISEA =', ISEA
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
VA(ISP,JSEA) = VA(ISP,JSEA) / CG(IK,ISEA) * CLATS(ISEA)
END DO
END DO
VAOLD = VA(:,1:NSEAL)
!/DEBUGSRC DO JSEA=1,NSEAL
!/DEBUGSRC WRITE(740+IAPROC,*) 'JSEA=', JSEA
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VA)=', minval(VA(:,JSEA)), maxval(VA(:,JSEA)), sum(VA(:,JSEA))
!/DEBUGSRC END DO
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(np) just defined", np)
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(npa) just defined", npa)
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'JACOBI_SOLVER, step 4'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FSSOURCE=', FSSOURCE
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FSREFRACTION=', FSREFRACTION
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FSFREQSHIFT=', FSFREQSHIFT
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FSGEOADVECT=', FSGEOADVECT
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'DTG=', DTG
!
! init matrix and right hand side
!
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 2'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
IF (IMEM == 1) THEN
ASPAR_JAC = ZERO
ELSE IF (IMEM == 2) THEN
ASPAR_DIAG_ALL = ZERO
ENDIF
B_JAC = ZERO
!/DEBUGSOLVER !WRITE(740+IAPROC,'(A20,20E20.10)') 'SUM BJAC INIT', sum(B_JAC), SUM(ASPAR_JAC)
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 3'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
! source terms
!
IF (FSSOURCE .and. B_JGS_NLEVEL .eq. 0) THEN
IF (B_JGS_SOURCE_NONLINEAR) THEN
call ADD_SOURCE_TERMS_NONLINEAR(DTG)
ELSE
call CALCARRAY_JACOBI_SOURCE(DTG,ASPAR_DIAG_ALL)
ENDIF
END IF
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 4'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
! geographical advection
!
IF (IMEM == 1) call calcARRAY_JACOBI2(DTG,FACX,FACY,VGX,VGY)
!/DEBUGSOLVER !WRITE(740+IAPROC,'(A20,20E20.10)') 'SUM BJAC 1', sum(B_JAC), SUM(ASPAR_JAC)
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 5'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
!/DEBUGSOLVER !WRITE(740+IAPROC,'(A20,20E20.10)') 'SUM BJAC 1', sum(B_JAC), SUM(ASPAR_JAC)
!
! spectral advection
!
IF (FSFREQSHIFT .or. FSREFRACTION) THEN
call calcARRAY_JACOBI_SPECTRAL(DTG,ASPAR_DIAG_ALL)
END IF
!/DEBUGSOLVER !WRITE(740+IAPROC,'(A20,20E20.10)') 'SUM BJAC 2', sum(B_JAC), SUM(ASPAR_JAC)
!
CALL APPLY_BOUNDARY_CONDITION
!/DEBUGSOLVER !WRITE(740+IAPROC,'(A20,20E20.10)') 'SUM BJAC 3', sum(B_JAC), SUM(ASPAR_JAC)
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 6'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(B_JAC, "B_JAC after calcARRAY", np)
!/DEBUGSOLVERCOH DO IP=1,npa
!/DEBUGSOLVERCOH TheArr(:, IP)=REAL(ASPAR_JAC(:, PDLIB_I_DIAG(IP)))
!/DEBUGSOLVERCOH END DO
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(TheArr, "ASPAR(:,I_DIAG) after calArr", np)
!/DEBUGSOLVER !WRITE(740+IAPROC,'(A20,20E20.10)') 'SUM BJAC 4', sum(B_JAC), SUM(ASPAR_JAC)
nbIter=0
Lconverged = .false.
!
DO
IF (B_JGS_TERMINATE_DIFFERENCE) is_converged=0
IF (FSSOURCE .and. (B_JGS_NLEVEL .eq. 1) .and. B_JGS_SOURCE_NONLINEAR) CALL ADD_SOURCE_TERMS_NONLINEAR(DTG)
ICOUNT1 = 0
ICOUNT2 = 0
!AR:
!ADD LOOP OVER ELEMENTS HERE FOR ALL ASPAR STUFF
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION SOLVER LOOP 1'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
DO IP=1,np
IP_glob=iplg(IP)
JSEA=JX_TO_JSEA(IP)
ISEA=MAPFS(1,IP_glob)
eSI=PDLIB_SI(IP)
ACLOC=VA(:,JSEA)
IF (.NOT. LCONVERGED(IP)) THEN
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'Begin loop'
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'IP/IP_glob/ISEA/JSEA=', IP, IP_glob, ISEA, JSEA
!/DEBUGSRC WRITE(740+IAPROC,*) 'IP=', IP, ' IP_glob=', IP_glob
!/DEBUGSRC WRITE(740+IAPROC,*) 'sum(VA)in=', sum(VA(:,IP))
!/DEBUGFREQSHIFT DO ISP=1,NSPEC
!/DEBUGFREQSHIFT VAold(ISP) = VA(ISP,JSEA)
!/DEBUGFREQSHIFT IK=MAPWN(ISP)
!/DEBUGFREQSHIFT VAinput(ISP) = DBLE(CG(IK,ISEA)/CLATS(ISEA)) * VA(ISP, IP)
!/DEBUGFREQSHIFT VAacloc(ISP) = DBLE(CG(IK,ISEA)/CLATS(ISEA)) * ACLOC(ISP)
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'sum(VAold/VAinput/VAacloc)=', sum(VAold), sum(VAinput), sum(VAacloc)
!!/DEBUGFREQSHIFT DO ISP=1,NSPEC
!!/DEBUGFREQSHIFT eVal1 = eSI * VA(ISP,IP)
!!/DEBUGFREQSHIFT eVal2 = B_JAC(ISP,IP)
!!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'eVal12=', eVal1, eVal2
!!/DEBUGFREQSHIFT END DO
Sum_Prev=sum(ACLOC)
IF (IMEM == 2) THEN
CALL calcARRAY_JACOBI4(IP,ICOUNT2,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
!WRITE(*,'(A10,10F20.10)') 'JAC4', SUM(ASPAR_DIAG_LOCAL), SUM(ASPAR_OFF_DIAG_LOCAL), SUM(B_JAC_LOCAL)
!CALL calcARRAY_JACOBI3(IP,ICOUNT1,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
!WRITE(*,'(A10,10F20.10)') 'JAC3', SUM(ASPAR_DIAG_LOCAL), SUM(ASPAR_OFF_DIAG_LOCAL), SUM(B_JAC_LOCAL)
ASPAR_DIAG = ASPAR_DIAG_LOCAL + ASPAR_DIAG_ALL(:,IP)
!IF (ANY(ABS(ASPAR_DIAG) .LT. TINY(1.))) THEN
! WRITE(*,'(8I10,4F20.10)') IP, JSEA, ISEA, NSEA, NSEAL, np, npa, IP_glob, SUM(ASPAR_DIAG), SUM(ASPAR_DIAG_LOCAL), SUM(ASPAR_DIAG_ALL(:,IP)), SUM(B_JAC(:,IP))
! CALL PDLIB_ABORT(25)
!ENDIF
esum = B_JAC_LOCAL - ASPAR_OFF_DIAG_LOCAL + B_JAC(:,IP)
ELSEIF (IMEM == 1) THEN
!CALL calcARRAY_JACOBI4(IP,ICOUNT2,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
!WRITE(*,'(A10,10F20.10)') 'JAC4', SUM(ASPAR_DIAG_LOCAL), SUM(ASPAR_OFF_DIAG_LOCAL), SUM(B_JAC_LOCAL)
!CALL calcARRAY_JACOBI3(IP,ICOUNT1,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
!WRITE(*,'(A10,10F20.10)') 'JAC3', SUM(ASPAR_DIAG_LOCAL), SUM(ASPAR_OFF_DIAG_LOCAL), SUM(B_JAC_LOCAL)
eSum = B_JAC(:,IP)
ASPAR_DIAG = ASPAR_JAC(:,PDLIB_I_DIAG(IP))
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'eSI=', eSI
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'sum(ASPAR_DIAG)=', sum(ASPAR_DIAG)
!/DEBUGSRC WRITE(740+IAPROC,*) 'Step 1: sum(eSum)=', sum(eSum)
!/DEBUGSOLVERCOH eOff=ZERO
DO i = PDLIB_IA_P(IP)+1, PDLIB_IA_P(IP+1)
JP=PDLIB_JA(I)
IF (JP .ne. IP) THEN
eProd = ASPAR_JAC(:,i)*VA(:,JP)
eSum = eSum - eProd
!/DEBUGSOLVER WRITE(740+IAPROC,'(A20,3I10,20E20.10)') 'OFF DIAGONAL', IP, i, jp, sum(B_JAC(:,IP)), sum(eSum), SUM(ASPAR_JAC(:,i)), SUM(VA(:,JP))
!/DEBUGSOLVERCOH eOff=eOff + abs(ASPAR_JAC(:,i))
END IF
END DO
ENDIF ! IMEM
!/DEBUGSOLVERCOH OffDiag(:, IP)=REAL(eOff)
!/DEBUGSOLVERCOH WRITE(740+IAPROC,*) 'Step 2: sum(eSum)=', sum(eSum), ' eOff=', sum(eOff)
IF (FSREFRACTION) THEN
!/DEBUGREFRACTION WRITE(740+IAPROC,*) 'Adding refraction terms to eSum'
CAD=CAD_THE(:,IP)
DO ISP=1,NSPEC
ISPprevDir=ListISPprevDir(ISP)
ISPnextDir=ListISPnextDir(ISP)
eA_THE = - DTG*eSI*MAX(ZERO,CAD(ISPprevDir))
eC_THE = DTG*eSI*MIN(ZERO,CAD(ISPnextDir))
eSum(ISP) = eSum(ISP) - eA_THE * VA(ISPprevDir,IP)
eSum(ISP) = eSum(ISP) - eC_THE * VA(ISPnextDir,IP)
END DO
END IF
!/DEBUGSRC WRITE(740+IAPROC,*) 'Step 3: sum(eSum)=', sum(eSum)
IF (FSFREQSHIFT) THEN
IF (FreqShiftMethod .eq. 1) THEN
CAS=CAS_SIG(:,IP)
CP_SIG = MAX(ZERO,CAS)
CM_SIG = MIN(ZERO,CAS)
DO IK=0, NK
DMM(IK+1) = DBLE(WN(IK+1,ISEA) - WN(IK,ISEA))
END DO
DMM(NK+2) = ZERO
DMM(0)=DMM(1)
DO ITH=1,NTH
DO IK=2,NK
ISP = ITH + (IK -1)*NTH
ISPm1 = ITH + (IK-1 -1)*NTH
eFactM1 = CG(IK-1,ISEA) / CG(IK,ISEA)
eA_SIG = - eSI * CP_SIG(ISPm1)/DMM(IK-1) * eFactM1
eSum(ISP) = eSum(ISP) - eA_SIG*VA(ISPm1,IP)
END DO
DO IK=1,NK-1
ISP = ITH + (IK -1)*NTH
ISPp1 = ITH + (IK+1 -1)*NTH
eFactP1 = CG(IK+1,ISEA) / CG(IK,ISEA)
eC_SIG = eSI * CM_SIG(ISPp1)/DMM(IK) * eFactP1
eSum(ISP) = eSum(ISP) - eC_SIG*VA(ISPp1,IP)
END DO
END DO
ELSE IF (FreqShiftMethod .eq. 2) THEN
CWNB_M2=CWNB_SIG_M2(:,IP)
DO IK=1, NK
DWNI_M2(IK) = DBLE( CG(IK,ISEA) / DSIP(IK) )
END DO
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'Before FreqShift oper eSum=', sum(abs(eSum))
DO ITH=1,NTH
DO IK=2,NK
ISP = ITH + (IK -1)*NTH
ISPm1 = ITH + (IK-1 -1)*NTH
eFactM1 = DBLE( CG(IK-1,ISEA) / CG(IK,ISEA) )
eA_SIG = - eSI * DWNI_M2(IK) * MAX(CWNB_M2(ISPm1),ZERO) *eFactM1
eSum(ISP) = eSum(ISP) - eA_SIG*VA(ISPm1,IP)
END DO
DO IK=1,NK-1
ISP = ITH + (IK -1)*NTH
ISPp1 = ITH + (IK+1 -1)*NTH
eFactP1 = DBLE( CG(IK+1,ISEA) / CG(IK,ISEA) )
eC_SIG = eSI * DWNI_M2(IK) * MIN(CWNB_M2(ISP),ZERO) * eFactP1
eSum(ISP) = eSum(ISP) - eC_SIG*VA(ISPp1,IP)
END DO
END DO
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) ' after FreqShift oper eSum=', sum(abs(eSum))
END IF
END IF
!/DEBUGSRC WRITE(740+IAPROC,*) 'Step 4: sum(eSum)=', sum(eSum)
!/DEBUGSOLVERCOH PRE_VA(:, IP)=REAL(eSum)
eSum=eSum/ASPAR_DIAG
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'JSEA=', JSEA, ' nbIter=', nbIter
!/DEBUGFREQSHIFT DO ISP=1,NSPEC
!/DEBUGFREQSHIFT IK=MAPWN(ISP)
!/DEBUGFREQSHIFT VAnew(ISP) = DBLE(CG(IK,ISEA)/CLATS(ISEA)) * eSum(ISP)
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT DO ISP=1,NSPEC
!/DEBUGFREQSHIFT VAAnew(ISP) = VAnew(ISP)
!/DEBUGFREQSHIFT VAAacloc(ISP) = VAacloc(ISP)
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT DO ITH=1,NTH
!/DEBUGFREQSHIFT VAAnew(ITH + NSPEC) = FACHFA * VAAnew(ITH + NSPEC - NTH)
!/DEBUGFREQSHIFT VAAnew(ITH - NTH ) = 0.
!/DEBUGFREQSHIFT VAAacloc(ITH + NSPEC) = FACHFA * VAAacloc(ITH + NSPEC - NTH)
!/DEBUGFREQSHIFT VAAacloc(ITH - NTH ) = 0.
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT DO ISP=1-NTH,NSPEC
!/DEBUGFREQSHIFT VFLWN(ISP) = MAX(CWNB_M2(ISP),0.) * VAAnew(ISP) + MIN(CWNB_M2(ISP),0.) * VAAnew(ISP + NTH)
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT DO ISP=1,NSPEC
!/DEBUGFREQSHIFT eDiff(ISP) = VAnew(ISP) - VAold(ISP) - DWNI_M2(MAPWN(ISP)) * (VFLWN(ISP-NTH) - VFLWN(ISP) )
!/DEBUGFREQSHIFT eVal1=MAX(CWNB_M2(ISP-NTH),0.) * VAAacloc(ISP-NTH) + MIN(CWNB_M2(ISP-NTH),0.) * VAAnew(ISP)
!/DEBUGFREQSHIFT eVal2=MAX(CWNB_M2(ISP),0.) * VAAnew(ISP) + MIN(CWNB_M2(ISP),0.) * VAAacloc(ISP + NTH)
!/DEBUGFREQSHIFT eDiffB(ISP) = VAnew(ISP) - VAold(ISP) - DWNI_M2(MAPWN(ISP)) * (eVal1 - eVal2)
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT IF (ISEA .eq. 190) THEN
!/DEBUGFREQSHIFT DO IK=1,NK
!/DEBUGFREQSHIFT DO ITH=1,NTH
!/DEBUGFREQSHIFT ISP = ITH + (IK-1)*NTH
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'ISP/ITH/IK=', ISP, ITH, IK
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'eDiff(A/B)=', eDiff(ISP), eDiffB(ISP)
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT END IF
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'NK=', NK, ' NTH=', NTH
!/DEBUGFREQSHIFT eSumPart=0
!/DEBUGFREQSHIFT DO IK=1,NK
!/DEBUGFREQSHIFT DO ITH=1,NTH
!/DEBUGFREQSHIFT ISP = ITH + (IK-1)*NTH
!/DEBUGFREQSHIFT eSumPart = eSumPart + abs(eDiff(ISP))
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT IF (ISEA .eq. 190) THEN
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'IK=', IK, ' eSumDiff=', eSumPart
!/DEBUGFREQSHIFT END IF
!/DEBUGFREQSHIFT END DO
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'sum(eDiff/VAnew/VAold)=', sum(abs(eDiff)), sum(abs(VAnew)), sum(abs(VAold))
Sum_New=sum(eSum)
IF (B_JGS_LIMITER) THEN
CALL ACTION_LIMITER_LOCAL(IP, eSum, ACLOC, DTG)
END IF
IF (B_JGS_BLOCK_GAUSS_SEIDEL) THEN
VA(:,IP)=eSum
ELSE
U_JAC(:,IP)=eSum
END IF
ELSE
esum = VA(:,IP)
DO I = 1, PDLIB_CCON(IP)
ICOUNT1 = ICOUNT1 + 1
ICOUNT2 = ICOUNT2 + 1
END DO
ENDIF ! LCONVERGED
!
IF (B_JGS_TERMINATE_DIFFERENCE) THEN
if (Sum_new .gt. thr8) then
DiffNew=sum(abs(ACLOC - eSum))
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'DiffNew=', DiffNew, ' Sum_new=', Sum_new
!DiffOld=abs(Sum_prev - Sum_new)
p_is_converged = DiffNew/Sum_new
!/DEBUGSOLVER write(740+IAPROC,'(10E20.10)') p_is_converged, sum(ASPAR_DIAG), SUM(B_JAC(:,IP)), DiffNew, Sum_new, sum(ACLOC), sum(esum)
else
p_is_converged = zero
endif
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'p_is_converged=', p_is_converged
IF (p_is_converged .lt. B_JGS_DIFF_THR) then
is_converged = is_converged + 1
lconverged(ip) = .true.
ENDIF
END IF
!IF (IP == 2) STOP
!/DEBUGSRC WRITE(740+IAPROC,*) 'sum(VA)out=', sum(VA(:,IP))
END DO ! IP
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION SOLVER LOOP 2'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!/DEBUGSOLVERCOH WRITE(740+IAPROC,*) 'sum(OffDiag)=', sum(OffDiag)
!/DEBUGSOLVERCOH WRITE (eFile,40) nbIter
!/DEBUGSOLVERCOH 40 FORMAT ('PRE_VA_',i4.4,'.txt')
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(OffDiag, "OffDiag(np) just check", np)
!/DEBUGSOLVERCOH ! CALL WRITE_VAR_TO_TEXT_FILE(PRE_VA, eFile)
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(PRE_VA, "PRE_VA(np) just check", np)
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(np) before exchanges", np)
IF (B_JGS_BLOCK_GAUSS_SEIDEL) THEN
extmp = VA(:,1:NPA)
CALL PDLIB_exchange2DREAL(extmp)
VA(:,1:NPA) = extmp
ELSE
CALL PDLIB_exchange2DREAL(U_JAC)
VA(:,1:NPA) = U_JAC
END IF
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION SOLVER LOOP 3'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(npa) after exchanges", npa)
!
! Terminate via number of iteration
!
IF (B_JGS_TERMINATE_MAXITER) THEN
nbIter=nbIter+1
IF (nbIter .gt. B_JGS_MAXITER) THEN
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Exiting by TERMINATE_MAXITER'
EXIT
END IF
END IF
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION SOLVER LOOP 4'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
! Terminate via differences
!
IF (B_JGS_TERMINATE_DIFFERENCE) THEN
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.
!write(*,*) prop_conv, nbIter, is_converged
if (myrank == 0) WRITE(*,*) 'solver', nbiter, is_converged, prop_conv, B_JGS_PMIN
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'solver', nbiter, is_converged, prop_conv, B_JGS_PMIN
!/DEBUGSOLVER FLUSH(740+IAPROC)
IF (prop_conv .le. B_JGS_PMIN) THEN
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'prop_conv=', prop_conv
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'NX=', NX
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'is_converged=', is_converged
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'Exiting by TERMINATE_DIFFERENCE'
EXIT
END IF
END IF
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION SOLVER LOOP 5'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!
! Terminate via norm
!
IF (B_JGS_TERMINATE_NORM) THEN
Sum_L2 =0
IF (FSSOURCE .and. (B_JGS_NLEVEL .eq. 1) .and. B_JGS_SOURCE_NONLINEAR) CALL ADD_SOURCE_TERMS_NONLINEAR(DTG)
DO IP = 1, np
IP_glob=iplg(IP)
IF ((MAPSTA(1,IP_glob).eq.1).and.(IOBP(IP_glob).eq.1)) THEN
JSEA=JX_TO_JSEA(IP)
eSI=PDLIB_SI(IP)
eSum=B_JAC(:,IP)
ACLOC=VA(:,IP)
eSum(:) = eSum(:) - ASPAR_DIAG(:)*ACLOC
DO I = PDLIB_IA_P(IP)+1, PDLIB_IA_P(IP+1)
JP=PDLIB_JA(I)
eSum(:) = eSum(:) - ASPAR_JAC(:,i)*VA(:,JP)
END DO
IF (FSREFRACTION) THEN
CAD=CAD_THE(:,IP)
DO ISP=1,NSPEC
ISPprevDir=ListISPprevDir(ISP)
ISPnextDir=ListISPnextDir(ISP)
eA_THE = - DTG*eSI*MAX(ZERO,CAD(ISPprevDir))
eC_THE = DTG*eSI*MIN(ZERO,CAD(ISPnextDir))
eSum(ISP) = eSum(ISP) - eA_THE*VA(ISPprevDir,IP)
eSum(ISP) = eSum(ISP) - eC_THE*VA(ISPnextDir,IP)
END DO
END IF
IF (FSFREQSHIFT) THEN
CAS=CAS_SIG(:,IP)
CP_SIG = MAX(ZERO,CAS)
CM_SIG = MIN(ZERO,CAS)
DO ITH=1,NTH
IF (IOBPD(ITH,IP_glob) .NE. 0) THEN
DO IK=2,NK
ISP =ITH + (IK -1)*NTH
ISPm1=ITH + (IK-1-1)*NTH
eFactM1=CG(IK-1,ISEA) / CG(IK,ISEA)
eA_SIG= - eSI*CP_SIG(ISPm1)/DMM(IK-1) * eFactM1
eSum(ISP) = eSum(ISP) - eA_SIG*VA(ISPm1,IP)
END DO
DO IK=1,NK-1
ISP =ITH + (IK -1)*NTH
ISPp1=ITH + (IK+1-1)*NTH
eFactP1=CG(IK+1,ISEA) / CG(IK,ISEA)
eC_SIG= eSI*CM_SIG(ISPp1)/DMM(IK) * eFactP1
eSum(ISP) = eSum(ISP) - eC_SIG*VA(ISPp1,IP)
END DO
END IF
END DO
END IF
Sum_L2 = Sum_L2 + sum(eSum*eSum)
END IF
END DO
CALL MPI_ALLREDUCE(Sum_L2, Sum_L2_GL, 1, rtype, MPI_SUM, MPI_COMM_WCMP, ierr)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'Sum_L2_gl=', Sum_L2_gl
!/DEBUGSOLVER FLUSH(740+IAPROC)
IF (Sum_L2_gl .le. B_JGS_NORM_THR) THEN
!/DEBUGFREQSHIFT WRITE(740+IAPROC,*) 'Exiting by TERMINATE_NORM'
EXIT
END IF
END IF
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION SOLVER LOOP 6'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
END DO ! Open Do Loop ... End of Time Interval
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'nbIter=', nbIter, ' B_JGS_MAXITER=', B_JGS_MAXITER
!/DEBUGSOLVER FLUSH(740+IAPROC)
! Tihs is below also goes into the matrix ... like the wave boundary ...
DO IP = 1, npa
IP_glob=iplg(IP)
!/DEBUGSRC WRITE(740+IAPROC,*) 'IOBPD loop, Before, sum(VA)=', sum(VA(:,IP))
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
VA(ISP,IP)=MAX(ZERO, VA(ISP,IP))*IOBDP(IP_glob)*DBLE(IOBPD(ITH,IP_glob))
END DO
!/DEBUGSRC WRITE(740+IAPROC,*) 'IOBPD loop, After, sum(VA)=', sum(VA(:,IP))
END DO
!!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(npa) after loop", npa)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FLBPI=', FLBPI
!/DEBUGSOLVER FLUSH(740+IAPROC)
DO JSEA=1, NSEAL
IP = JSEA
IP_glob = iplg(IP)
ISEA=MAPFS(1,IP_glob)
!
!/DEBUGSRC IntDiff=0
!/DEBUGSRC SumVS=0
!/DEBUGSRC SumVD=0
!/DEBUGSRC SumVAin=0
!/DEBUGSRC SumVAout=0
!/DEBUGSRC SumVAw3srce=0
!/DEBUGSRC SumACout=0
!
DO ISP=1,NSPEC
IK = 1 + (ISP-1)/NTH
eVA = MAX ( ZERO ,CG(IK,ISEA)/CLATS(ISEA)*REAL(VA(ISP,IP)) )
!/DEBUGSRC SumACout=SumACout + REAL(VA(ISP,IP))
!/DEBUGSRC VS_w3srce = VSTOT(ISP,JSEA) * DTG / MAX(1., (1. - DTG*VDTOT(ISP,JSEA)))
!/DEBUGSRC eVA_w3srce = MAX(0., VA(ISP,JSEA) + VS_w3srce)
!/DEBUGSRC IntDiff = IntDiff + abs(eVA - eVA_w3srce)
!/DEBUGSRC ACsolve=B_JAC(ISP,IP)/ASPAR_JAC(ISP,PDLIB_I_DIAG(IP))
!/DEBUGSRC eB=VA(ISP,JSEA) + DTG*(VSTOT(ISP,JSEA) - VDTOT(ISP,JSEA)*VA(ISP,JSEA))
!/DEBUGSRC eVAsolve=MAX(0., CG(IK,ISEA)/CLATS(ISEA)*ACsolve)
!/DEBUGSRC VAsolve(ISP)=eVAsolve
!/DEBUGSRC SumVS = SumVS + abs(VSTOT(ISP,JSEA))
!/DEBUGSRC SumVD = SumVD + abs(VDTOT(ISP,JSEA))
!/DEBUGSRC SumVAin = SumVAin + abs(VA(ISP,JSEA))
!/DEBUGSRC SumVAout = SumVAout + abs(eVA)
!/DEBUGSRC SumVAw3srce = SumVAw3srce + abs(eVA_w3srce)
VA(ISP,JSEA) = eVA
END DO
!/DEBUGSRC WRITE(740+IAPROC,*) 'ISEA=', ISEA, ' IntDiff=', IntDiff, ' DTG=', DTG
!/DEBUGSRC IF (ISEA .eq. TESTNODE) THEN
!/DEBUGSRC DO ISP=1,NSPEC
!/DEBUGSRC WRITE(740+IAPROC,*) 'ISP=', ISP, 'VA/VAsolve=', VA(ISP,JSEA), VAsolve(ISP)
!/DEBUGSRC END DO
!/DEBUGSRC END IF
!/DEBUGSRC WRITE(740+IAPROC,*) 'SHAVE=', SHAVETOT(JSEA)
!/DEBUGSRC WRITE(740+IAPROC,*) 'Sum(VS/VD)=', SumVS, SumVD
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VS)=', minval(VSTOT(:,JSEA)), maxval(VSTOT(:,JSEA)), sum(VSTOT(:,JSEA))
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VD)=', minval(VDTOT(:,JSEA)), maxval(VDTOT(:,JSEA)), sum(VDTOT(:,JSEA))
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VA)=', minval(VA(:,JSEA)), maxval(VA(:,JSEA)), sum(VA(:,JSEA))
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VAsolve)=', minval(VAsolve), maxval(VAsolve), sum(VAsolve)
!/DEBUGSRC WRITE(740+IAPROC,*) 'SumVA(in/out/w3srce)=', SumVAin, SumVAout, SumVAw3srce
!/DEBUGSRC WRITE(740+IAPROC,*) 'SumACout=', SumACout
END DO ! JSEA
!/MEMCHECK write(740+IAPROC,*) 'memcheck_____:', 'WW3_PROP SECTION 7'
!/MEMCHECK call getMallocInfo(mallinfos)
!/MEMCHECK call printMallInfo(IAPROC,mallInfos)
!/DEBUGSRC DO JSEA=1,NSEAL
!/DEBUGSRC WRITE(740+IAPROC,*) 'JSEA=', JSEA
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VA)=', minval(VA(:,JSEA)), maxval(VA(:,JSEA)), sum(VA(:,JSEA))
!/DEBUGSRC END DO
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VAtot)=', minval(VA), maxval(VA), sum(VA)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK, end'
!/DEBUGSOLVER FLUSH(740+IAPROC)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_EXPLICIT_BLOCK(FACX, FACY, DTG, VGX, VGY)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Explicit block solver
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!
!/S USE W3SERVMD, ONLY: STRACE
!/
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, NTRI
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
USE MPI, ONLY : MPI_MIN
use yowElementpool, ONLY: ne, INE
USE YOWNODEPOOL, ONLY: PDLIB_I_DIAG, PDLIB_IA_P, PDLIB_JA, np
USE YOWNODEPOOL, ONLY: PDLIB_SI, PDLIB_IEN, PDLIB_CCON, NPA, PDLIB_IE_CELL2, PDLIB_POS_CELL2
use yowDatapool, ONLY: rtype
use YOWNODEPOOL, ONLY: npa, iplg
use yowExchangeModule, ONLY : PDLIB_exchange2DREAL
USE W3ADATMD, ONLY: WN
USE MPI, ONLY : MPI_SUM, MPI_INT
USE W3ADATMD, ONLY: MPI_COMM_WCMP, CFLXYMAX
USE W3GDATMD, ONLY: IOBP, IOBPD, NSEA, SIG, IOBDP
USE W3GDATMD, ONLY: NK, NK2, NTH, ECOS, ESIN, NSPEC, MAPFS
USE W3WDATMD, ONLY: TIME
USE W3TIMEMD, ONLY: DSEC21
USE W3GDATMD, ONLY: NSEAL, CLATS, FACHFA
USE W3IDATMD, ONLY: FLCUR
!/DEBUGSRC USE W3WDATMD, ONLY: SHAVETOT
USE W3WDATMD, ONLY: VA, VSTOT, VDTOT
USE W3ADATMD, ONLY: CG, CX, CY, MPI_COMM_WCMP
USE W3ODATMD, ONLY: TBPIN, FLBPI, IAPROC
USE W3PARALL, ONLY : INIT_GET_JSEA_ISPROC, ZERO, THR8
USE W3PARALL, ONLY : ListISPprevDir, ListISPnextDir
USE W3PARALL, ONLY : JX_TO_JSEA
USE W3GDATMD, ONLY: B_JGS_NLEVEL
!
implicit none
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
!
INTEGER :: IP, ISP, ITH, IK, JSEA, ISEA, IP_glob, IE, IPOS
! for the exchange
REAL :: CCOS, CSIN, CCURX, CCURY
REAL :: eSum(NSPEC)
INTEGER :: ITER_EXP(nspec)
INTEGER :: ISPp1, ISPm1, JP, I1, I2, I3, NI(3), IT
INTEGER, SAVE :: ITER_MAX
REAL :: eFactM1, eFactP1
REAL :: Sum_Prev, Sum_New
REAL :: prop_conv, eSI, p_is_converged
REAL :: Sum_L2, Sum_L2_GL
REAL :: DMM(0:NK2)
REAL :: DiffNew, DTMAX_GLOBAL_EXP, DTMAX_EXP
REAL :: eDiff(NSPEC), eProd(NSPEC), u33(nspec,3)
REAL :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC), LAMBDA(NSPEC,3)
REAL :: VAnew(NSPEC), VAold(NSPEC), REST, CFLXY, DT4AI
REAL :: VAinput(NSPEC), VAacloc(NSPEC), eDiffB(NSPEC), KTMP(nspec,3), TMP(nspec)
REAL :: eDiffSing, eSumPart, N(nspec,ntri), kksum(nspec,npa), ST3(nspec), utilde33(nspec)
REAL :: FL11(NSPEC),FL12(NSPEC),FL21(NSPEC),FL22(NSPEC),FL31(NSPEC),FL32(NSPEC)
REAL :: FL111(NSPEC), FL112(NSPEC), FL211(NSPEC), FL212(NSPEC), FL311(NSPEC), FL312(NSPEC)
REAL :: KELEMGL(NSPEC,3,NTRI), FLALLGL(NSPEC,3,NTRI)
REAL :: eVal1, eVal2,thr
REAL :: eVA
!/DEBUGSRC REAL :: IntDiff, eVA_w3srce, eVAsolve, SumACout
!/DEBUGSRC REAL :: SumVAin, SumVAout, SumVAw3srce, SumVS, SumVD, VS_w3srce
!/DEBUGSRC REAL :: VAsolve(NSPEC)
!/DEBUGSRC REAL*8 :: ACsolve
!/DEBUGSRC REAL :: eB
REAL :: ASPAR_DIAG(NSPEC)
!/DEBUGSOLVERCOH REAL*8 :: PRE_VA(NSPEC, npa)
!/DEBUGSOLVERCOH REAL*8 :: OffDIAG(NSPEC, npa)
!/DEBUGSOLVERCOH REAL*8 :: eOff(NSPEC)
!/DEBUGSOLVERCOH REAL*8 :: eSum1(NSPEC), eSum2(NSPEC)
CHARACTER(len=128) eFile
INTEGER ierr, i
INTEGER JP_glob
INTEGER is_converged, itmp
thr = dble(tiny(1.))
CCURX = FACX
CCURY = FACY
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'EXPLICIT BLOCK SOLVER, begin'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NX=', NX
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NP=', NP
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NPA=', NPA
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NSEA=', NSEA
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'NSEAL=', NSEAL
!/DEBUGSOLVER FLUSH(740+IAPROC)
!
! 1.b Initialize arrays
!
CALL SETDEPTH_PDLIB
!
! 2. Calculate velocities ---------------- *
!
!/DEBUGSRC WRITE(740+IAPROC,*) 'NSEAL =', NSEAL
!/DEBUGSRC WRITE(740+IAPROC,*) 'NP =', NP
!/DEBUGSRC WRITE(740+IAPROC,*) 'NPA =', NPA
DO JSEA=1,NSEAL
IP = JSEA
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
!!/DEBUGSRC WRITE(740+IAPROC,*) 'IP =', IP
!!/DEBUGSRC WRITE(740+IAPROC,*) 'IP_glob =', IP_glob
!!/DEBUGSRC WRITE(740+IAPROC,*) 'ISEA =', ISEA
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
CCOS = FACX * ECOS(ITH)
CSIN = FACY * ESIN(ITH)
VA(ISP,IP) = DBLE(VA(ISP,JSEA) / CG(IK,ISEA) * CLATS(ISEA))
!/MGP VLCFLX(ISP,IP) = VLCFLX(ISP,IP) - CCURX*VGX/CLATS(ISEA)
!/MGP VLCFLY(ISP,IP) = VLCFLY(ISP,IP) - CCURY*VGY
END DO
END DO
!/DEBUGSRC DO IP=1,NP
!/DEBUGSRC WRITE(740+IAPROC,*) 'IP=', IP
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VA)=', minval(VA(:,IP)), maxval(VA(:,IP)), sum(VA(:,IP))
!/DEBUGSRC END DO
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(np) just defined", np)
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(npa) just defined", npa)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'JACOBI_SOLVER, step 3'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FLCUR=', FLCUR
!/DEBUGSOLVER FLUSH(740+IAPROC)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'EXPLICIT BLOCK SOLVER, step 4'
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'min,max(0)=', 0
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'min,max(0)=', 0
!/DEBUGSOLVER FLUSH(740+IAPROC)
DO IE = 1, NE
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
!LAMBDA(:,1) = 1./6. *(C(:,I1,1)+C(:,I2,1)+C(:,I3,1))
!LAMBDA(:,2) = 1./6. *(C(:,I1,2)+C(:,I2,2)+C(:,I3,2))
KELEMGL(:,1,IE) = LAMBDA(:,1) * PDLIB_IEN(1,IE) + LAMBDA(:,2) * PDLIB_IEN(2,IE)
KELEMGL(:,2,IE) = LAMBDA(:,1) * PDLIB_IEN(3,IE) + LAMBDA(:,2) * PDLIB_IEN(4,IE)
KELEMGL(:,3,IE) = LAMBDA(:,1) * PDLIB_IEN(5,IE) + LAMBDA(:,2) * PDLIB_IEN(6,IE)
KTMP(:,1) = KELEMGL(:,1,IE)
KTMP(:,2) = KELEMGL(:,2,IE)
KTMP(:,3) = KELEMGL(:,3,IE)
TMP(:) = SUM(MIN(ZERO,KTMP(:,:)),DIM=2)
N(:,IE) = -1.d0/MIN(-THR,TMP(:))
KELEMGL(:,1,IE) = MAX(ZERO,KTMP(:,1))
KELEMGL(:,2,IE) = MAX(ZERO,KTMP(:,2))
KELEMGL(:,3,IE) = MAX(ZERO,KTMP(:,3))
! FL11 = C(:,I2,1) * PDLIB_IEN(1,IE) + C(:,I2,2) * PDLIB_IEN(2,IE)
! FL12 = C(:,I3,1) * PDLIB_IEN(1,IE) + C(:,I3,2) * PDLIB_IEN(2,IE)
! FL21 = C(:,I3,1) * PDLIB_IEN(3,IE) + C(:,I3,2) * PDLIB_IEN(4,IE)
! FL22 = C(:,I1,1) * PDLIB_IEN(3,IE) + C(:,I1,2) * PDLIB_IEN(4,IE)
! FL31 = C(:,I1,1) * PDLIB_IEN(5,IE) + C(:,I1,2) * PDLIB_IEN(6,IE)
! FL32 = C(:,I2,1) * PDLIB_IEN(5,IE) + C(:,I2,2) * PDLIB_IEN(6,IE)
FL111 = 2*FL11+FL12
FL112 = 2*FL12+FL11
FL211 = 2*FL21+FL22
FL212 = 2*FL22+FL21
FL311 = 2*FL31+FL32
FL312 = 2*FL32+FL31
FLALLGL(:,1,IE) = (FL311 + FL212) * 1./6. + KELEMGL(:,1,IE)
FLALLGL(:,2,IE) = (FL111 + FL312) * 1./6. + KELEMGL(:,2,IE)
FLALLGL(:,3,IE) = (FL211 + FL112) * 1./6. + KELEMGL(:,3,IE)
END DO
KKSUM = 0.d0
DO IE = 1, NE
NI = INE(:,IE)
KKSUM(:,NI) = KKSUM(:,NI) + KELEMGL(:,:,IE)
END DO
DTMAX_GLOBAL_EXP = 1.d0/THR
DO IP = 1, NP
DTMAX_EXP = PDLIB_SI(IP)/MAX(THR,MAXVAL(KKSUM(:,IP)))
DTMAX_GLOBAL_EXP = MIN ( DTMAX_GLOBAL_EXP, DTMAX_EXP)
CFLXYMAX(IP) = DBLE(DTG)/DTMAX_EXP
END DO
rest = CFLXYMAX(1)
DO IP = 2, NP
if (rest .lt. CFLXYMAX(IP)) then
rest = CFLXYMAX(IP)
iter_max = ip
endif
END DO
DTMAX_EXP=DTMAX_GLOBAL_EXP
call mpi_allreduce(DTMAX_EXP,DTMAX_GLOBAL_EXP,1,rtype,MPI_MIN,MPI_COMM_WCMP,ierr)
CFLXY = DTG/DTMAX_GLOBAL_EXP
REST = ABS(MOD(CFLXY,1.d0))
IF (REST .LT. THR) THEN
ITER_MAX = ABS(NINT(CFLXY))
ELSE IF (REST .GT. THR .AND. REST .LT. 0.5d0) THEN
ITER_MAX = ABS(NINT(CFLXY)) + 1
ELSE
ITER_MAX = ABS(NINT(CFLXY))
END IF
DT4AI = DTG/ITER_MAX
DO IT = 1, ITER_MAX
DO IP = 1, NPA
ST3 = ZERO
DO I = 1, PDLIB_CCON(IP)
IE = PDLIB_IE_CELL2(IP,I)
U33 = VA(:,INE(:,IE))
UTILDE33 = N(:,IE)*(FLALLGL(:,1,IE)*U33(:,1)+FLALLGL(:,2,IE)*U33(:,2)+FLALLGL(:,3,IE)*U33(:,3))
IPOS = PDLIB_POS_CELL2(IP,I)
ST3 = ST3 + KELEMGL(:,IPOS,IE)*(U33(:,IPOS)-UTILDE33)
END DO
VA(:,IP) = MAX(ZERO,VA(:,IP)-DT4AI/PDLIB_SI(IP)*ST3)
END DO !IP
CALL PDLIB_exchange2DREAL(VA(:,1:NPA))
END DO !IT
DO IP = 1, npa
IP_glob=iplg(IP)
!/DEBUGSRC WRITE(740+IAPROC,*) 'IOBPD loop, Before, sum(VA)=', sum(VA(:,IP))
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
VA(ISP,IP)=MAX(ZERO, VA(ISP,IP))*IOBDP(IP_glob)*DBLE(IOBPD(ITH,IP_glob))
END DO
!/DEBUGSRC WRITE(740+IAPROC,*) 'IOBPD loop, After, sum(VA)=', sum(VA(:,IP))
END DO
!/DEBUGSOLVERCOH CALL CHECK_ARRAY_INTEGRAL_NX_R8(VA, "VA(npa) after loop", npa)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'FLBPI=', FLBPI
!/DEBUGSOLVER FLUSH(740+IAPROC)
DO JSEA=1, NSEAL
IP=JSEA
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
!
!/DEBUGSRC IntDiff=0
!/DEBUGSRC SumVS=0
!/DEBUGSRC SumVD=0
!/DEBUGSRC SumVAin=0
!/DEBUGSRC SumVAout=0
!/DEBUGSRC SumVAw3srce=0
!/DEBUGSRC SumACout=0
!
DO ISP=1,NSPEC
IK = 1 + (ISP-1)/NTH
eVA = MAX ( 0. ,CG(IK,ISEA)/CLATS(ISEA)*REAL(VA(ISP,IP)) )
!/DEBUGSRC SumACout=SumACout + REAL(VA(ISP,IP))
!/DEBUGSRC VS_w3srce = VSTOT(ISP,JSEA) * DTG / MAX(1., (1. - DTG*VDTOT(ISP,JSEA)))
!/DEBUGSRC eVA_w3srce = MAX(0., VA(ISP,JSEA) + VS_w3srce)
!/DEBUGSRC IntDiff = IntDiff + abs(eVA - eVA_w3srce)
!/DEBUGSRC ACsolve=B_JAC(ISP,IP)/ASPAR_JAC(ISP,PDLIB_I_DIAG(IP))
!/DEBUGSRC eB=VA(ISP,JSEA) + DTG*(VSTOT(ISP,JSEA) - VDTOT(ISP,JSEA)*VA(ISP,JSEA))
!/DEBUGSRC eVAsolve=MAX(0., CG(IK,ISEA)/CLATS(ISEA)*ACsolve)
!/DEBUGSRC VAsolve(ISP)=eVAsolve
!/DEBUGSRC SumVS = SumVS + abs(VSTOT(ISP,JSEA))
!/DEBUGSRC SumVD = SumVD + abs(VDTOT(ISP,JSEA))
!/DEBUGSRC SumVAin = SumVAin + abs(VA(ISP,JSEA))
!/DEBUGSRC SumVAout = SumVAout + abs(eVA)
!/DEBUGSRC SumVAw3srce = SumVAw3srce + abs(eVA_w3srce)
VA(ISP,JSEA) = eVA
END DO
!/DEBUGSRC WRITE(740+IAPROC,*) 'ISEA=', ISEA, ' IntDiff=', IntDiff, ' DTG=', DTG
!/DEBUGSRC IF (ISEA .eq. TESTNODE) THEN
!/DEBUGSRC DO ISP=1,NSPEC
!/DEBUGSRC WRITE(740+IAPROC,*) 'ISP=', ISP, 'VA/VAsolve=', VA(ISP,JSEA), VAsolve(ISP)
!/DEBUGSRC END DO
!/DEBUGSRC END IF
!/DEBUGSRC WRITE(740+IAPROC,*) 'SHAVE=', SHAVETOT(JSEA)
!/DEBUGSRC WRITE(740+IAPROC,*) 'Sum(VS/VD)=', SumVS, SumVD
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VS)=', minval(VSTOT(:,JSEA)), maxval(VSTOT(:,JSEA)), sum(VSTOT(:,JSEA))
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VD)=', minval(VDTOT(:,JSEA)), maxval(VDTOT(:,JSEA)), sum(VDTOT(:,JSEA))
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VA)=', minval(VA(:,JSEA)), maxval(VA(:,JSEA)), sum(VA(:,JSEA))
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VAsolve)=', minval(VAsolve), maxval(VAsolve), sum(VAsolve)
!/DEBUGSRC WRITE(740+IAPROC,*) 'SumVA(in/out/w3srce)=', SumVAin, SumVAout, SumVAw3srce
!/DEBUGSRC WRITE(740+IAPROC,*) 'SumACout=', SumACout
END DO
!/DEBUGSRC DO IP=1,NP
!/DEBUGSRC WRITE(740+IAPROC,*) 'IP=', IP
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VA)=', minval(VA(:,IP)), maxval(VA(:,IP)), sum(VA(:,IP))
!/DEBUGSRC END DO
!/DEBUGSRC WRITE(740+IAPROC,*) 'min/max/sum(VAtot)=', minval(VA), maxval(VA), sum(VA)
!/DEBUGSOLVER WRITE(740+IAPROC,*) 'PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK, end'
!/DEBUGSOLVER FLUSH(740+IAPROC)
END SUBROUTINE
!/ ------------------------------------------------------------------- /
SUBROUTINE BLOCK_SOLVER_INIT
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Initialization of the block solver
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
USE CONSTANTS, ONLY : LPDLIB
USE W3GDATMD, ONLY: MAPSF, NSEAL, DMIN, IOBDP, MAPSTA, IOBP, MAPFS, NX
USE W3ADATMD, ONLY: DW
USE W3PARALL, ONLY: INIT_GET_ISEA
USE YOWNODEPOOL, ONLY: iplg, np
USE yowfunction, ONLY: pdlib_abort
use YOWNODEPOOL, ONLY: npa
USE W3GDATMD, ONLY: B_JGS_USE_JACOBI
USE W3PARALL, ONLY : ListISPprevDir, ListISPnextDir
USE W3PARALL, ONLY : ListISPprevFreq, ListISPnextFreq
USE W3GDATMD, ONLY: NSPEC, NTH, NK
USE W3GDATMD, ONLY: FSTOTALIMP
USE W3ODATMD, ONLY: IAPROC
!/
IMPLICIT NONE
!
!/ ------------------------------------------------------------------- /
!/
INTEGER ISP, ITH, IK, ISPprevFreq, ISPnextFreq
INTEGER NewISP, JTH, istat
!/DEBUGINIT WRITE(740+IAPROC,*) 'BLOCK_SOLVER_INIT, step 1'
!/DEBUGINIT FLUSH(740+IAPROC)
ALLOCATE(ListISPnextDir(NSPEC), ListISPprevDir(NSPEC), ListISPnextFreq(NSPEC), ListISPprevFreq(NSPEC),stat=istat)
IF (istat /= 0) CALL PDLIB_ABORT(8)
!/DEBUGINIT WRITE(740+IAPROC,*) 'BLOCK_SOLVER_INIT, step 2'
!/DEBUGINIT FLUSH(740+IAPROC)
DO ISP=1,NSPEC
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
IF (IK .eq. 1) THEN
ISPprevFreq=-1
ELSE
ISPprevFreq=iTH + (IK-1 -1)*NTH
END IF
ListISPprevFreq(ISP)=ISPprevFreq
IF (IK .eq. NK) THEN
ISPnextFreq=-1
ELSE
ISPnextFreq=iTH + (IK+1 -1)*NTH
END IF
ListISPnextFreq(ISP)=ISPnextFreq
!
IF (ITH .eq. 1) THEN
JTH=NTH
ELSE
JTH=ITH-1
ENDIF
NewISP=JTH + (IK-1)*NTH
ListISPprevDir(ISP)=NewISP
IF (ITH .eq. NTH) THEN
JTH=1
ELSE
JTH=ITH+1
ENDIF
NewISP=JTH + (IK-1)*NTH
ListISPnextDir(ISP)=NewISP
END DO
!/DEBUGINIT WRITE(740+IAPROC,*) 'BLOCK_SOLVER_INIT, step 3'
!/DEBUGINIT FLUSH(740+IAPROC)
IF (FSTOTALIMP .and. B_JGS_USE_JACOBI) THEN
!/DEBUGINIT WRITE(740+IAPROC,*) 'BLOCK_SOLVER_INIT, step 4'
!/DEBUGINIT FLUSH(740+IAPROC)
CALL JACOBI_INIT
!/DEBUGINIT WRITE(740+IAPROC,*) 'BLOCK_SOLVER_INIT, step 5'
!/DEBUGINIT FLUSH(740+IAPROC)
END IF
!/DEBUGINIT WRITE(740+IAPROC,*) 'BLOCK_SOLVER_INIT, step 6'
!/DEBUGINIT FLUSH(740+IAPROC)
END SUBROUTINE
!/ ------------------------------------------------------------------ /
SUBROUTINE SETDEPTH_PDLIB
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Set depth pointer
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE CONSTANTS, ONLY : LPDLIB
USE W3GDATMD, ONLY: MAPSF, NSEAL, DMIN, IOBDP, MAPSTA, IOBP, MAPFS, NX
USE W3ADATMD, ONLY: DW
USE W3PARALL, ONLY: INIT_GET_ISEA
USE YOWNODEPOOL, ONLY: iplg, np
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!
INTEGER :: JSEA, ISEA, IX, IP, IP_glob
REAL*8, PARAMETER :: DTHR = 10E-6
!/S CALL STRACE (IENT, 'SETDEPTH_PDLIB')
IOBDP = 1
DO JSEA=1,NSEAL
IP=JSEA
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
IF (DW(ISEA) .LT. DMIN + DTHR) IOBDP(IP_glob) = 0
!WRITE(*,*) ip, ip_glob, MAPSTA(1,IP_glob), IOBP(IP_glob), DW(ISEA), DMIN
END DO
!/
!/ End of SETDEPTH_PDLIB --------------------------------------------- /
!/
END SUBROUTINE SETDEPTH_PDLIB
!/ ------------------------------------------------------------------- /
SUBROUTINE BLOCK_SOLVER_FINALIZE
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Finalize Solver
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY: B_JGS_USE_JACOBI
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!
!/S CALL STRACE (IENT, 'BLOCK_SOLVER_FINALIZE')
IF (B_JGS_USE_JACOBI) THEN
CALL JACOBI_FINALIZE
END IF
!/
!/ End of SETDEPTH_PDLIB --------------------------------------------- /
!/
END SUBROUTINE BLOCK_SOLVER_FINALIZE
!/ ------------------------------------------------------------------- /
SUBROUTINE JACOBI_INIT
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Init jacobi solver
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE W3GDATMD, ONLY: NSPEC, B_JGS_BLOCK_GAUSS_SEIDEL
use YOWNODEPOOL, ONLY: PDLIB_NNZ, npa, np
USE yowfunction, ONLY: pdlib_abort
USE W3GDATMD, ONLY: NTH, NK, NSEAL
USE W3PARALL, ONLY: IMEM
!/DEBUGINIT USE W3ODATMD, ONLY : IAPROC
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
INTEGER istat
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 1'
!/DEBUGINIT FLUSH(740+IAPROC)
IF (IMEM == 1) THEN
ALLOCATE(ASPAR_JAC(NSPEC, PDLIB_NNZ), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(9)
ENDIF
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 2'
!/DEBUGINIT FLUSH(740+IAPROC)
ALLOCATE(B_JAC(NSPEC,npa), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(10)
ALLOCATE(ASPAR_DIAG_SOURCES(NSPEC,npa), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(10)
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 3'
!/DEBUGINIT FLUSH(740+IAPROC)
ALLOCATE(CAD_THE(NSPEC,NSEAL), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(11)
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 4'
!/DEBUGINIT FLUSH(740+IAPROC)
CAD_THE = 0
IF (FreqShiftMethod .eq. 1) THEN
ALLOCATE(CAS_SIG(NSPEC,NSEAL), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(11)
CAS_SIG = 0
END IF
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 5, FreqShiftMethod=', FreqShiftMethod
!/DEBUGINIT FLUSH(740+IAPROC)
IF (FreqShiftMethod .eq. 2) THEN
!/DEBUGINIT WRITE(740+IAPROC,*) 'Before CWNB_SIG_M2 allocation, NTH=', NTH
!/DEBUGINIT FLUSH(740+IAPROC)
ALLOCATE(CWNB_SIG_M2(1-NTH:NSPEC,NSEAL), stat=istat)
!/DEBUGINIT WRITE(740+IAPROC,*) 'After CWNB_SIG_M2 allocation, istat=', istat
!/DEBUGINIT FLUSH(740+IAPROC)
if(istat /= 0) CALL PDLIB_ABORT(11)
!/DEBUGINIT WRITE(740+IAPROC,*) 'After istat test'
!/DEBUGINIT FLUSH(740+IAPROC)
CWNB_SIG_M2(:,:) = 0
!/DEBUGINIT WRITE(740+IAPROC,*) 'After CWNB_SIG_M2 setting to zero'
!/DEBUGINIT FLUSH(740+IAPROC)
END IF
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 6'
!/DEBUGINIT FLUSH(740+IAPROC)
IF (.NOT. B_JGS_BLOCK_GAUSS_SEIDEL) THEN
ALLOCATE(U_JAC(NSPEC,npa), stat=istat)
if(istat /= 0) CALL PDLIB_ABORT(12)
END IF
!/DEBUGINIT WRITE(740+IAPROC,*) 'JACOBI_INIT, step 7'
!/DEBUGINIT FLUSH(740+IAPROC)
!/
!/ End of JACOBI_INIT ------------------------------------------------ /
!/
END SUBROUTINE JACOBI_INIT
!/ ------------------------------------------------------------------- /
SUBROUTINE JACOBI_FINALIZE
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 01-June-2018 |
!/ +-----------------------------------+
!/
!/ 01-June-2018 : Origination. ( version 6.04 )
!/
! 1. Purpose : Finalize jacobi solver
! 2. Method :
! 3. Parameters :
!
! 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 :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETER
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'JACOBI_FINALIZE')
DEALLOCATE(ASPAR_JAC, B_JAC)
!/
!/ End of JACOBI_FINALIZE -------------------------------------------- /
!/
END SUBROUTINE JACOBI_FINALIZE
!/
!/ End of Module PDLIB_W3PROFSMD ------------------------------------- /
!/
END MODULE PDLIB_W3PROFSMD