w3profsmd Module Reference

Functions/Subroutines
subroutine	w3xypug (ISP, FACX, FACY, DTG, VQ, VGX, VGY, LCALC)
subroutine	w3cflug (ISEA, NKCFL, FACX, FACY, DT, MAPFS, CFLXYMAX, VGX, VGY)
subroutine	w3xypfsn2 (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3xypfspsi2 (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3xypfsnimp (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3xypfsfct2 (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	setdepth

Function/Subroutine Documentation

setdepth()

subroutine w3profsmd::setdepth

Definition at line 1587 of file w3profsmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | 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 :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    USE constants, ONLY : lpdlib
    USE w3gdatmd, ONLY: mapsf, nseal, dmin, iobdp, mapsta, iobp, mapfs, nx
    USE w3adatmd, ONLY: dw
 
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local PARAMETERs
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !
    INTEGER :: JSEA, ISEA, IX, IP
    real*8, PARAMETER :: dthr = 10e-6
#ifdef W3_S
    CALL strace (ient, 'SETDEPTH')
#endif
    iobdp = 1
    DO ip=1,nx
      IF (dw(ip) .LT. dmin + dthr) iobdp(ip) = 0
      !WRITE(*,*) ip, ip_glob, MAPSTA(1,IP_glob), IOBP(IP_glob), DW(ISEA), DMIN
    END DO
 

References w3gdatmd::dmin, w3adatmd::dw, w3gdatmd::iobdp, w3gdatmd::iobp, constants::lpdlib, w3gdatmd::mapfs, w3gdatmd::mapsf, w3gdatmd::mapsta, w3gdatmd::nseal, w3gdatmd::nx, and w3servmd::strace().

Referenced by w3xypug().

w3cflug()

subroutine w3profsmd::w3cflug	(	integer, intent(in)	ISEA,
		integer, intent(in)	NKCFL,
		real, intent(in)	FACX,
		real, intent(in)	FACY,
		real, intent(in)	DT,
		integer, dimension(ny*nx), intent(in)	MAPFS,
		real, intent(inout)	CFLXYMAX,
		real, intent(in)	VGX,
		real, intent(in)	VGY
	)

Definition at line 272 of file w3profsmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           Fabrice Ardhuin         |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jun-2013 |
    !/                  +-----------------------------------+
    !/
    !/    01-Mar-2011 : Origination.                        ( version 3.14 )
    !/    20-Jun-2013 : Computes only up to NKCFL for tests ( version 4.10 )
    !/     1-Jun-2017 : Rewrite routine for performance     ( version 5.xx )
    !/
    !  1. Purpose :
    !
    !     Computes the max CFL number for output purposes
    !
    !  2. Method :
    !
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISEA     Int.   I   Index of sea point
    !       NKCFL    Int.   I   Maximum frequency index
    !       FACX/Y   Real   I   Factor in propagation velocity.
    !                          ( 1 or 0 * DT / DX )
    !       DT       Real   I   Time step.
    !       MAPFS    I.A.   I   Storage map.
    !       CFLXYMAX Real       Maxmimum CFL for spatial advection
    !       VGX/Y    Real   I   Speed of grid.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       VCFL0X  R.A.  Local courant numbers for absolute group vel.
    !                     using local X-grid step.
    !       VCFL0Y  R.A.  Id. in Y.
    !     ----------------------------------------------------------------
    !
    !  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, nsea, mapsf, dtcfl, clats,          &
         flcx, flcy, nk, nth, dth, xfr,              &
         ecos, esin, sig,  pfmove,ien, index_cell,   &
         ntri, trigp, ccon ,                         &
         ie_cell, pos_cell, countri, si, iobp
 
    USE w3adatmd, ONLY: cg, cx, cy, atrnx, atrny, itime, dw
    USE w3idatmd, ONLY: flcur
#ifdef W3_T
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin,       &
         isbpi, bbpi0, bbpin
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
 
    IMPLICIT NONE
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISEA, NKCFL, MAPFS(NY*NX)
    REAL, INTENT(IN)        :: FACX, FACY, DT, VGX, VGY
    REAL, INTENT(INOUT)     :: CFLXYMAX
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK
    INTEGER                 :: IP, IP2, ISEA2, I, J, IE, IV, I1, I2, I3
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CCOS, CSIN, CCURX, CCURY
    REAL                    :: C(NX,2)
    real*8                  :: kelem(3), ktmp(3), lambda(2)
    real*8                  :: kksum, dtmaxexp
    !/
    !/  Velocities
    !/
    real*8,  PARAMETER :: onesixth  = 1.0d0/6.0d0
    real*8,  PARAMETER :: thr8      = tiny(1.0d0)
    REAL,    PARAMETER :: THR       = tiny(1.0)
 
    !/ ------------------------------------------------------------------- /
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Set constants
    !
 
#ifdef W3_S
    CALL strace (ient, 'W3CFLUG')
#endif
    cflxymax=1e-10
    ip = mapsf(isea,3)
    !
    ! CFL no important on boundary
    !
    IF (iobp(ip).EQ.1) THEN
      ccurx  = facx
      ccury  = facy
      !
      ! Loop over spectral components
      !
      DO ik=1,nkcfl
        DO ith=1,nth
          ccos   = facx * ecos(ith)
          csin   = facy * esin(ith)
          c(:,:)=0.
 
          !
          ! 2.  Calculate advection velocities: group speed ---------------- *
          !
          !AR: needs to be rewritten for speed ... single node computation is costly ...
          !MA: you are right but now it is only called if CFX and UNST for CFL profiling
 
          DO i = index_cell(ip), index_cell(ip+1)-1
            ie=ie_cell(i)       ! TRIGP(IV,IE)=IP  with IV=POS_CELL(I)
            DO j=1,3
              ip2=trigp(j,ie)
              isea2=mapfs(ip2)
              c(ip2,1) = ccos * cg(ik,isea2) / clats(isea2)
              c(ip2,2) = csin * cg(ik,isea2)
#ifdef W3_MGP
              c(ip2,1) = c(ip2,1) - ccurx*vgx/clats(isea2)
              c(ip2,2) = c(ip2,2) - ccury*vgy
#endif
              IF ( flcur ) THEN
                IF (iobp(ip2) .EQ. 1) THEN
                  c(ip2,1) = c(ip2,1) + ccurx*cx(isea2)/clats(isea2)
                  c(ip2,2) = c(ip2,2) + ccury*cy(isea2)
                END IF
              END IF  ! end of ( FLCUR )
            END DO
          END DO
          !
          !3.     Calculate K-Values and contour based quantities ...
          !
          kksum = 0.d0
          DO i = index_cell(ip), index_cell(ip+1)-1
            ie=ie_cell(i)       ! TRIGP(IV,IE)=IP
            iv=pos_cell(i)
            i1 = trigp(1,ie)
            i2 = trigp(2,ie)
            i3 = trigp(3,ie)
            lambda(1) = onesixth *(c(i1,1)+c(i2,1)+c(i3,1)) ! Advection speed in X direction
            lambda(2) = onesixth *(c(i1,2)+c(i2,2)+c(i3,2)) ! Advection speed in Y direction
            kelem(1) = lambda(1) * ien(ie,1) + lambda(2) * ien(ie,2) ! K-Values - so called Flux Jacobians
            kelem(2) = lambda(1) * ien(ie,3) + lambda(2) * ien(ie,4) ! K-Values - so called Flux Jacobians
            kelem(3) = lambda(1) * ien(ie,5) + lambda(2) * ien(ie,6) ! K-Values - so called Flux Jacobians
 
            ktmp        = kelem ! Copy
            kelem       = max(0.d0,ktmp)
 
            kksum = kksum  +  kelem(iv)
          END DO ! COUNTRI
          !
          dtmaxexp = si(ip)/max(dble(10.e-10),kksum)
          cflxymax = max(dble(dt)/dtmaxexp,dble(cflxymax))
        END DO
      END DO
    END IF
    !
    RETURN

References w3adatmd::atrnx, w3adatmd::atrny, w3odatmd::bbpi0, w3odatmd::bbpin, w3gdatmd::ccon, w3adatmd::cg, w3gdatmd::clats, w3gdatmd::countri, w3adatmd::cx, w3adatmd::cy, w3timemd::dsec21(), w3gdatmd::dtcfl, w3gdatmd::dth, w3adatmd::dw, w3gdatmd::ecos, w3gdatmd::esin, w3odatmd::flbpi, w3idatmd::flcur, w3gdatmd::flcx, w3gdatmd::flcy, w3gdatmd::ie_cell, w3gdatmd::ien, w3gdatmd::index_cell, w3gdatmd::iobp, w3odatmd::isbpi, w3adatmd::itime, w3gdatmd::mapsf, w3odatmd::nbi, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nsea, w3gdatmd::nth, w3gdatmd::ntri, w3gdatmd::nx, w3gdatmd::ny, w3gdatmd::pfmove, w3gdatmd::pos_cell, w3gdatmd::si, w3gdatmd::sig, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3gdatmd::trigp, and w3gdatmd::xfr.

Referenced by w3wavemd::w3wave().

w3xypfsfct2()

subroutine w3profsmd::w3xypfsfct2	(	integer, intent(in)	ISP,
		real, dimension(:,:), intent(in)	C,
		logical, intent(in)	LCALC,
		real, intent(in)	RD10,
		real, intent(in)	RD20,
		real, intent(in)	DT,
		double precision, dimension(:), intent(inout)	AC
	)

Definition at line 1282 of file w3profsmd.F90.

    !/
    !/
    !/                  +-----------------------------------+
    !/                  |  WWIII Version of the WWM FS Code |
    !/                  |  by Aron Roland                   |
    !/                  |  for use in WWIII                 |
    !/                  |  GPL License                      |
    !/                  |  Last update :        19-Dec-2007 |
    !/                  +-----------------------------------+
    !/
    !  1. Purpose :
    !     Advection of a scalar in a arbitary velocity field on unstructured meshes
    !     for the conservative hyperbolic equation N,t + (c*N),xy = 0 in spatial space
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Subroutine tracing (!/S switch)
    !
    !  5. Called by :
    !
    !       W3XYPUG   Routine for advection on unstructured grid
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE w3gdatmd, ONLY : nk, nth, ntri, nx, ccon, ie_cell,pos_cell, si, &
         ien, trigp, clats, mapsf, iobpd, iobpa, tria, iobdp
#ifdef W3_REF1
    USE w3gdatmd, ONLY : refpars
#endif
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, iter
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin, isbpi, bbpi0, bbpin
    USE w3timemd, ONLY: dsec21
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    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(:,:)              ! Velocity field in it's X- and Y- Components,
    DOUBLE PRECISION,    INTENT(INOUT) :: AC(:)               ! 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
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
 
    real*8,  PARAMETER :: onesixth  = 1.0d0/6.0d0
    real*8,  PARAMETER :: onethird  = 1.0d0/3.0d0
    real*8,  PARAMETER :: thr8      = tiny(1.0d0)
    REAL,    PARAMETER :: THR       = tiny(1.0)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    !
    ! local integer
    !
    INTEGER :: IP, IE, IT, I1, I2, I3, I, ITH, IK
    INTEGER :: IBI, NI(3)
    !
    ! local real
    !
    REAL    :: RD1, RD2
    !:
    ! local double
    !
    real*8  :: utilde, boundary_forcing
    real*8  :: ft, cflxy
    real*8  :: fl11, fl12, fl21, fl22, fl31, fl32
    real*8  :: fl111, fl112, fl211, fl212, fl311, fl312
    real*8  :: dtsi(nx), u(nx), dt4ai
    real*8  :: dtmaxgl, dtmaxexp, rest
    real*8  :: lambda(2), ktmp(3), tmp(3)
    real*8  :: bet1(3), betahat(3)
    real*8  :: theta_l(3,ntri), theta_h(3,ntri), theta_ace(3,ntri), utmp(3)
    real*8  :: wii(2,nx), ul(nx), ustari(2,nx)
 
    real*8  :: pm(nx), pp(nx), uim(nx), uip(nx)
 
    real*8  :: kelem(3,ntri), flall(3,ntri)
    real*8  :: kksum(nx), st(nx), beta
    real*8  :: nm(ntri)
 
#ifdef W3_S
    CALL strace (ient, 'W3XYPFSFCT2')
#endif
 
    ! 1. initialisation
 
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    dtmaxgl = dble(10.e10)
    !
    !2       Propagation
    !2.a     Calculate K-Values and contour based quantities ...
    !
    DO ie = 1, ntri ! I precacalculate this arrays below as I assume that current velocity  changes continusly ...
      i1 = trigp(1,ie) ! Index of the Element Nodes (TRIGP)
      i2 = trigp(2,ie)
      i3 = trigp(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) * ien(ie,1) + lambda(2) * ien(ie,2) ! K-Values - so called Flux Jacobians
      kelem(2,ie) = lambda(1) * ien(ie,3) + lambda(2) * ien(ie,4)
      kelem(3,ie) = lambda(1) * ien(ie,5) + lambda(2) * ien(ie,6)
      ktmp        = kelem(:,ie) ! Copy
      nm(ie)      = - 1.d0/min(-thr8,sum(min(0.d0,ktmp))) ! N-Values
      fl11  = c(i2,1) * ien(ie,1) + c(i2,2) * ien(ie,2) ! Weights for Simpson Integration
      fl12  = c(i3,1) * ien(ie,1) + c(i3,2) * ien(ie,2)
      fl21  = c(i3,1) * ien(ie,3) + c(i3,2) * ien(ie,4)
      fl22  = c(i1,1) * ien(ie,3) + c(i1,2) * ien(ie,4)
      fl31  = c(i1,1) * ien(ie,5) + c(i1,2) * ien(ie,6)
      fl32  = c(i2,1) * ien(ie,5) + c(i2,2) * ien(ie,6)
      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 ! NTRI
 
    IF (lcalc) THEN ! 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
      kksum = 0.d0
      DO ie = 1, ntri
        ni = trigp(:,ie)
        kksum(ni) = kksum(ni) + kelem(:,ie)
      END DO ! IE
      dtmaxexp = 1e10 ! initialize to large number
      DO ip = 1, nx
        dtmaxexp = si(ip)/max(dble(10.e-10),kksum(ip)*iobdp(ip))
        dtmaxgl  = min( dtmaxgl, dtmaxexp)
      END DO ! IP
      cflxy = dble(dt)/dtmaxgl
      rest  = abs(mod(cflxy,1.0d0))
      IF (rest .LT. thr8) THEN
        iter(ik,ith) = abs(nint(cflxy))
      ELSE IF (rest .GT. thr8 .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/si(:) ! Some precalculations for the time integration.
 
    u = ac
    ul = u
    !
    !  Now loop on sub-timesteps
    !
    DO it = 1, iter(ik,ith)
 
      st = 0.d0
 
      DO ie = 1, ntri
        ni      = trigp(:,ie)
        utmp    = u(ni)
        ft      =  - onesixth*dot_product(utmp,flall(:,ie))
        tmp     =  max(0.d0,kelem(:,ie))
        utilde  =  nm(ie) * ( dot_product(tmp,utmp) - ft )
        theta_l(:,ie) =  tmp * ( utmp - utilde )
        IF (abs(ft) .GT. dble(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(0.d0,min(betahat(1),1.d0-betahat(2),1.d0))
            bet1(2)       = max(0.d0,min(betahat(2),1.d0-betahat(3),1.d0))
            bet1(3)       = max(0.d0,min(betahat(3),1.d0-betahat(1),1.d0))
            theta_l(:,ie) = ft * bet1
          END IF
        ELSE
          theta_l(:,ie) = 0.d0
        END IF
        !              THETA_H(:,IE) = (ONETHIRD+DT4AI/(2.d0*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/si(ni)
        st(ni)          = st(ni) + theta_l(:,ie)*dt4ai/si(ni)
      END DO
 
      !            UL          = MAX(0.d0,U-ST)*DBLE(IOBPD(ITH,:))!*DBLE(IOBDP(:)) ... add for IOBDP dry/wet flag.
 
      DO ip = 1,nx
        ul(ip) = max(0.d0,u(ip)-st(ip))*dble(iobpd(ith,ip))
      END DO
 
      ustari(1,:) = max(ul,u)
      ustari(2,:) = min(ul,u)
 
      uip = -thr8
      uim =  thr8
      pp  = 0.d0
      pm  = 0.d0
      DO ie = 1, ntri
        ni = trigp(:,ie)
        pp(ni)  = pp(ni) + max(  thr8, -theta_ace(:,ie))
        pm(ni)  = pm(ni) + min( -thr8, -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)/max( thr8,pp))
      wii(2,:) = min(1.0d0,(uim-ul)/min(-thr8,pm))
 
      st = 0.d0
      DO ie = 1, ntri
        DO i = 1, 3
          ip = trigp(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 = trigp(:,ie)
        st(ni) = st(ni) + beta * theta_ace(:,ie)
      END DO
 
      DO ip = 1,nx
        !
        ! IOBPD is the switch for removing energy coming from the shoreline
        !
        u(ip) = max(0.d0,ul(ip)-st(ip))*dble(iobpd(ith,ip))
#ifdef W3_REF1
        IF (refpars(3).LT.0.5.AND.iobpd(ith,ip).EQ.0.AND.iobpa(ip).EQ.0) u(ip)= ac(ip) ! restores reflected boundary values
#endif
      END DO
      !
      ! update spectrum
      ac = u
      !
      ! 4 Update boundaries: performs interpolation in time as done in rect grids (e.g. w3pro1md.ftn)
      !
      IF ( flbpi ) THEN
        !
        ! 4.1 In this case the boundary is read from the nest.ww3 file
        !
        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
        !
        ! Overwrites only the incoming energy ( IOBPD(ITH,IP) = 0)
        !
        DO ibi=1, nbi
          ip = mapsf(isbpi(ibi),1)
          ac(ip) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
               / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
        END DO
 
      ENDIF
      !
    END DO ! End of loop on time steps
    !      CALL EXTCDE ( 99 )
    !/
    !/ End of W3XYPFSN ----------------------------------------------------- /
    !/

References w3odatmd::bbpi0, w3odatmd::bbpin, w3gdatmd::ccon, w3adatmd::cg, w3gdatmd::clats, w3timemd::dsec21(), w3odatmd::flbpi, w3gdatmd::ie_cell, w3gdatmd::ien, w3gdatmd::iobdp, w3gdatmd::iobpa, w3gdatmd::iobpd, w3odatmd::isbpi, w3adatmd::iter, w3gdatmd::mapsf, w3odatmd::nbi, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nth, w3gdatmd::ntri, w3gdatmd::nx, w3gdatmd::pos_cell, w3gdatmd::refpars, w3gdatmd::si, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, w3gdatmd::tria, and w3gdatmd::trigp.

Referenced by w3xypug().

w3xypfsn2()

subroutine w3profsmd::w3xypfsn2	(	integer, intent(in)	ISP,
		real, dimension(:,:), intent(in)	C,
		logical, intent(in)	LCALC,
		real, intent(in)	RD10,
		real, intent(in)	RD20,
		real, intent(in)	DT,
		double precision, dimension(:), intent(inout)	AC
	)

Definition at line 470 of file w3profsmd.F90.

 
    !/
    !/
    !/                  +-----------------------------------+
    !/                  |  WWIII Version of the WWM FS Code |
    !/                  |  by Aron Roland                   |
    !/                  |     and Fabrice  Ardhuin          |
    !/                  |  for use in WWIII                 |
    !/                  |  GPL License                      |
    !/                  |  Last update :        15-Apr-2020 |
    !/                  +-----------------------------------+
    !/
    !/    19-Dec-2007 : Origination.                        ( version 3.13 )
    !/    25-Aug-2011 : Change of method for IOBPD          ( version 4.04 )
    !/    03-Nov-2011 : Addition of shoreline reflection    ( version 4.04 )
    !/    15-Apr-2020 : Adds optional opt-out for CFL on BC ( version 7.08 )
    !/
    !/
    !  1. Purpose :
    !     Advection of a scalar in a arbitary velocity field on unstructured meshes
    !     for the conservative hyperbolic equation N,t + (c*N),xy = 0 in spatial space
    !     This is the standard explicit N-Scheme from Roe as formulated in Abgrall
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Subroutine tracing (!/S switch)
    !
    !  5. Called by :
    !
    !       W3XYPUG   Routine for advection on unstructured grid
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE w3gdatmd, ONLY : nk, nth, ntri, nx, ccon, ie_cell,pos_cell, si, &
         ien, trigp, clats, mapsf, iobpd, iobp, iobdp,  &
         iobpa, fsbccfl
#ifdef W3_REF1
    USE w3gdatmd, ONLY : refpars
#endif
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, iter, dw
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin, isbpi, bbpi0, bbpin
    USE w3timemd, ONLY: dsec21
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    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(:,:)                ! Velocity field in it's X- and Y- Components,
    DOUBLE PRECISION, INTENT(INOUT):: AC(:)         ! 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
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
 
    real*8,  PARAMETER :: onesixth  = 1.0d0/6.0d0
    real*8,  PARAMETER :: thr8      = tiny(1.0d0)
    REAL,    PARAMETER :: THR       = tiny(1.0)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    !
    ! local integer
    !
    INTEGER :: IP, IE, IT, I1, I2, I3, ITH, IK
    INTEGER :: IBI, NI(3)
    !
    ! local real
    !
    REAL    :: RD1, RD2
    !
    ! local double
    !
    real*8  :: utilde, boundary_forcing
    real*8  :: cflxy
    real*8  :: fl11, fl12, fl21, fl22, fl31, fl32
    real*8  :: fl111, fl112, fl211, fl212, fl311, fl312
    real*8  :: dtsi(nx), u(nx)
    real*8  :: dtmaxgl, dtmaxexp, rest
    real*8  :: lambda(2), ktmp(3), cloc(2,3)
    real*8  :: kelem(3,ntri), flall(3,ntri)
    real*8  :: kksum(nx), st(nx)
    real*8  :: nm(ntri)
 
#ifdef W3_S
    CALL strace (ient, 'W3XYPFSN')
#endif
 
    ! 1. initialisation
 
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    dtmaxgl = dble(10.e10)
    !
    !2       Propagation
    !2.a     Calculate K-Values and contour based quantities ...
    !
    DO ie = 1, ntri ! I precacalculate this arrays below as I assume that current velocity  changes continusly ...
 
      i1 = trigp(1,ie) ! Index of the Element Nodes (TRIGP)
      i2 = trigp(2,ie)
      i3 = trigp(3,ie)
      ni = trigp(:,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) * ien(ie,1) + lambda(2) * ien(ie,2) ! K-Values - so called Flux Jacobians
      kelem(2,ie) = lambda(1) * ien(ie,3) + lambda(2) * ien(ie,4)
      kelem(3,ie) = lambda(1) * ien(ie,5) + lambda(2) * ien(ie,6)
      !
      ktmp        = kelem(:,ie) ! Copy
      nm(ie)      = - 1.d0/min(-thr8,sum(min(0.d0,ktmp))) ! N-Values
      kelem(:,ie) = max(0.d0,ktmp)
      fl11  = c(i2,1) * ien(ie,1) + c(i2,2) * ien(ie,2) ! Weights for Simpson Integration
      fl12  = c(i3,1) * ien(ie,1) + c(i3,2) * ien(ie,2)
      fl21  = c(i3,1) * ien(ie,3) + c(i3,2) * ien(ie,4)
      fl22  = c(i1,1) * ien(ie,3) + c(i1,2) * ien(ie,4)
      fl31  = c(i1,1) * ien(ie,5) + c(i1,2) * ien(ie,6)
      fl32  = c(i2,1) * ien(ie,5) + c(i2,2) * ien(ie,6)
      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)
      ! IF (I1.EQ.1.OR.I2.EQ.1.OR.I3.EQ.1) WRITE(6,*) 'TEST N1 :',IK,ITH,IP,IE,KELEM(:,IE),'##',LAMBDA
    END DO ! NTRI
 
    IF (lcalc) THEN ! 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
      kksum = 0.d0
      DO ie = 1, ntri
        ni = trigp(:,ie)
        kksum(ni) = kksum(ni) + kelem(:,ie)
      END DO ! IE
      dtmaxexp = 1e10 ! initialize to large number
      DO ip = 1, nx
        IF (iobp(ip) .EQ. 1 .OR. fsbccfl) THEN
          dtmaxexp = si(ip)/max(dble(10.e-10),kksum(ip)*iobdp(ip))
          dtmaxgl  = min( dtmaxgl, dtmaxexp)
        END IF
      END DO ! IP
      cflxy = dble(dt)/dtmaxgl
      rest  = abs(mod(cflxy,1.0d0))
      IF (rest .LT. thr8) THEN
        iter(ik,ith) = abs(nint(cflxy))
      ELSE IF (rest .GT. thr8 .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
 
    DO ip = 1, nx
      dtsi(ip) = dble(dt)/dble(iter(ik,ith))/si(ip) ! Some precalculations for the time integration.
    END DO
 
    DO it = 1, iter(ik,ith)
      u = ac
      st = 0.d0
      DO ie = 1, ntri
        ni     = trigp(:,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
 
      DO ip = 1, nx
        !
        ! IOBPD=0  : waves coming from land (or outside grid)
        ! Possibly set flux to zero by multiplying ST by IOBPD but also in UTILDE multiply U(NI) by IOBPD ...
        !
        u(ip) = max(0.d0,u(ip)-dtsi(ip)*st(ip)*(1-iobpa(ip)))*dble(iobpd(ith,ip))
 
#ifdef W3_REF1
        WRITE(10111,*) refpars(3), iobpd(ith,ip), iobpa(ip), u(ip),  ac(ip)
        IF (refpars(3).LT.0.5.AND.iobpd(ith,ip).EQ.0.AND.iobpa(ip).EQ.0) THEN
          u(ip) = ac(ip) ! restores reflected boundary values
        ENDIF
#endif
      END DO
      ! update spectrum
      ac = u
      !
      ! 4 Update boundaries: performs interpolation in time as done in rect grids (e.g. w3pro1md.ftn)
      !
      IF ( flbpi ) THEN
        !
        ! 4.1 In this case the boundary is read from the nest.ww3 file
        !
        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
        !
        ! Overwrites only the incoming energy ( IOBPD(ITH,IP) = 0)
        !
        DO ibi=1, nbi
          ip = mapsf(isbpi(ibi),1)
          ac(ip) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
               / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
        END DO
 
      ENDIF
      !
    END DO ! End of loop on time steps
    !        CALL EXTCDE ( 99 )
    !/
    !/ End of W3XYPFSN ----------------------------------------------------- /
    !/

References w3odatmd::bbpi0, w3odatmd::bbpin, w3gdatmd::ccon, w3adatmd::cg, w3gdatmd::clats, w3timemd::dsec21(), w3adatmd::dw, w3odatmd::flbpi, w3gdatmd::fsbccfl, w3gdatmd::ie_cell, w3gdatmd::ien, w3gdatmd::iobdp, w3gdatmd::iobp, w3gdatmd::iobpa, w3gdatmd::iobpd, w3odatmd::isbpi, w3adatmd::iter, w3gdatmd::mapsf, w3odatmd::nbi, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nth, w3gdatmd::ntri, w3gdatmd::nx, w3gdatmd::pos_cell, w3gdatmd::refpars, w3gdatmd::si, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, and w3gdatmd::trigp.

Referenced by w3xypug().

w3xypfsnimp()

subroutine w3profsmd::w3xypfsnimp	(	integer, intent(in)	ISP,
		real, dimension(:,:), intent(in)	C,
		logical, intent(in)	LCALC,
		real, intent(in)	RD10,
		real, intent(in)	RD20,
		real, intent(in)	DT,
		double precision, dimension(:), intent(inout)	AC
	)

Definition at line 976 of file w3profsmd.F90.

 
    !/
    !/
    !/                  +-----------------------------------+
    !/                  |  WWIII Version of the WWM FS Code |
    !/                  |  by Aron Roland                   |
    !/                  |  for use in WWIII                 |
    !/                  |  GPL License                      |
    !/                  |  Last update :        15-Dec-2013 |
    !/                  +-----------------------------------+
    !/
    !/    15-Dec-2013 : Bug fix for implicit scheme         ( version 4.16 )
    !/
    !  1. Purpose :
    !     Advection of a scalar in a arbitary velocity field on unstructured meshes
    !     for the conservative hyperbolic equation N,t + (c*N),xy = 0 in spatial space
    !     This is the standard explicit N-Scheme from Roe as formulated in Abgrall
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Subroutine tracing (!/S switch)
    !
    !  5. Called by :
    !
    !       W3XYPUG   Routine for advection on unstructured grid
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE w3gdatmd, ONLY : nk, nth, ntri, nx, ccon, ie_cell,pos_cell, si, &
         ien, trigp, clats, mapsf, iobpd, iobpa, iobp, iaa, jaa, posi, &
         tria, nnz
#ifdef W3_REF1
    USE w3gdatmd, ONLY : refpars
#endif
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, iter
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin, isbpi, bbpi0, bbpin
    USE w3timemd, ONLY: dsec21
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    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(:,:)              ! Velocity field in it's X- and Y- Components,
    DOUBLE PRECISION,INTENT(INOUT) :: AC(:)         ! 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
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
 
    real*8,  PARAMETER :: onesixth  = 1.0d0/6.0d0
    real*8,  PARAMETER :: onethird  = 1.0d0/3.0d0
    real*8,  PARAMETER :: thr8      = tiny(1.0d0)
    REAL,    PARAMETER :: THR       = tiny(1.0)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    !
    ! local integer
    !
    INTEGER :: IP, IE, POS, I1, I2, I3, I, J, ITH, IK
    INTEGER :: IBI
    !
    ! local real
    !
    REAL    :: RD1, RD2
    !:
    ! local double
    !
    real*8  :: boundary_forcing
    real*8  :: fl11, fl12, fl21, fl22, fl31, fl32
    real*8  :: u(nx)
    real*8  :: dtmaxgl
    real*8  :: lambda(2)
    real*8  :: kp(3,ntri)
    real*8  :: k1, dtk, tmp3, km(3), k(3)
    real*8  :: nm(ntri), crfs(3), deltal(3,ntri)
    real*8  :: b(nx), x(nx)
    real*8  :: aspar(nnz)
 
    INTEGER :: IWKSP( 20*NX )
    INTEGER :: FLJU(NX)
    INTEGER :: FLJAU(NNZ+1)
 
 
    INTEGER :: POS_TRICK(3,2)
 
    INTEGER :: IPAR(16)
    INTEGER :: IERROR ! IWK                             ! ERROR Indicator and Work Array Size,
    INTEGER :: JAU(NNZ+1), JU(NX)
 
    real*8  :: fpar(16)  ! DROPTOL
    real*8  :: wksp( 8 * nx ) ! REAL WORKSPACES
    real*8  :: au(nnz+1)
    real*8  :: iniu(nx)
 
    external bcgstab
 
    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
 
#ifdef W3_S
    CALL strace (ient, 'W3XYPFSN')
#endif
 
    ! 1. initialisation
 
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    dtmaxgl = dble(10.e10)
 
    IF (.false.) THEN
      WRITE(*,*) 'NNZ', nnz
      WRITE(*,*) 'MINVAL IAA,JAA', minval(iaa), minval(jaa)
      WRITE(*,*) 'MINVAL IAA,JAA', maxval(iaa), maxval(jaa)
      WRITE(*,*) 'MAX/MIN POSI',   maxval(posi), minval(posi)
      WRITE(*,*) 'AC, AQ', sum(ac)
    END IF
    !
    !2       Propagation
    !2.a     Calculate K-Values and contour based quantities ...
    !
    DO ie = 1, ntri ! I precacalculate this arrays below as I assume that current velocity  changes continusly ...
      i1 = trigp(1,ie) ! Index of the Element Nodes (TRIGP)
      i2 = trigp(2,ie)
      i3 = trigp(3,ie)
      lambda(1) = onesixth * (c(i1,1)+c(i2,1)+c(i3,1))
      lambda(2) = onesixth * (c(i1,2)+c(i2,2)+c(i3,2))
      k(1)  = lambda(1) * ien(ie,1) + lambda(2) * ien(ie,2)
      k(2)  = lambda(1) * ien(ie,3) + lambda(2) * ien(ie,4)
      k(3)  = lambda(1) * ien(ie,5) + lambda(2) * ien(ie,6)
      kp(1,ie) = max(0.d0,k(1))
      kp(2,ie) = max(0.d0,k(2))
      kp(3,ie) = max(0.d0,k(3))
      km(1) = min(0.d0,k(1))
      km(2) = min(0.d0,k(2))
      km(3) = min(0.d0,k(3))
      fl11 = c(i2,1)*ien(ie,1)+c(i2,2)*ien(ie,2)
      fl12 = c(i3,1)*ien(ie,1)+c(i3,2)*ien(ie,2)
      fl21 = c(i3,1)*ien(ie,3)+c(i3,2)*ien(ie,4)
      fl22 = c(i1,1)*ien(ie,3)+c(i1,2)*ien(ie,4)
      fl31 = c(i1,1)*ien(ie,5)+c(i1,2)*ien(ie,6)
      fl32 = c(i2,1)*ien(ie,5)+c(i2,2)*ien(ie,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 )
      deltal(:,ie) = crfs(:)- kp(:,ie)
      nm(ie)       = 1.d0/min(dble(thr),sum(km(:)))
    END DO ! NTRI
 
    u = dble(ac)
    j = 0
    aspar = 0.d0
    b     = 0.d0
    DO ip = 1, nx
      DO i = 1, ccon(ip)
        j = j + 1
        ie    =  ie_cell(j)
        pos   =  pos_cell(j)
        k1    =  kp(pos,ie) * iobpd(ith,ip)
        IF (k1 > 0.) THEN
          dtk   =  k1 * dt
          tmp3  =  dtk * nm(ie)
          i1    =  posi(1,j)
          i2    =  posi(2,j)
          i3    =  posi(3,j)
          aspar(i1) =  onethird * tria(ie) + dtk - tmp3 * deltal(pos,ie)              + aspar(i1)
          aspar(i2) =                         - tmp3 * deltal(pos_trick(pos,1),ie) + aspar(i2)
          aspar(i3) =                         - tmp3 * deltal(pos_trick(pos,2),ie) + aspar(i3)
          b(ip)     =  b(ip) + onethird * tria(ie) * u(ip)
        ELSE
          i1    =  posi(1,j)
          aspar(i1) =  onethird * tria(ie) + aspar(i1)
          b(ip)     =  b(ip) + onethird * tria(ie) * u(ip)
        END IF
      END DO ! End loop over connected elements ...
    END DO
    !
    !2DO setup a semi-implicit integration scheme for source terms only ...
    !
    ipar(1) = 0       ! always 0 to start an iterative solver
    ipar(2) = 1       ! right preconditioning
    ipar(3) = 1       ! use convergence test scheme 1
    ipar(4) = 8*nx  !
    ipar(5) = 15
    ipar(6) = 1000    ! use at most 1000 matvec's
    fpar(1) = dble(1.0e-8)  ! relative tolerance 1.0E-6
    fpar(2) = dble(1.0e-10)  ! absolute tolerance 1.0E-10
    fpar(11) = 0.d0    ! clearing the FLOPS counter
 
    au    = 0.
    fljau = 0
    flju  = 0
    jau   = 0
    ju    = 0
 
    CALL ilu0  (nx, aspar, jaa, iaa, au, fljau, flju, iwksp, ierror)
 
    !         WRITE(*,*) 'PRECONDITIONER', IERROR
 
    !         IF (SUM(AC) .GT. 0.) THEN
    IF (.false.) THEN
      WRITE(*,*) sum(ac)
      WRITE(*,*) 'CALL SOLVER'
      WRITE(*,*) 'WRITE CG', sum(cg)
      WRITE(*,*) 'B, X, AC, U', sum(b), sum(x), sum(ac), sum(u)
      WRITE(*,*) 'IPAR, FPAR', sum(ipar), sum(fpar)
      WRITE(*,*) 'WKSP, INIU', sum(wksp), sum(iniu)
      WRITE(*,*) 'ASPAR, JAA, IAA',sum(aspar), sum(iaa), sum(jaa)
      WRITE(*,*) 'AU, FLJAU, FLJU',sum(au), sum(fljau), sum(flju)
    END IF
 
    iniu = u
    x    = 0.d0
 
    CALL runrc (nx, b, x, ipar, fpar, wksp, iniu, aspar, jaa, iaa, au, fljau, flju, bcgstab)
 
    IF (.false.) THEN
      WRITE(*,*) 'SOLUTION'
      WRITE(*,*) 'B, X, AC, U', sum(b), sum(x), sum(ac), sum(u)
      WRITE(*,*) 'IPAR, FPAR', sum(ipar), sum(fpar)
      WRITE(*,*) 'WKSP, INIU', sum(wksp), sum(iniu)
      WRITE(*,*) 'ASPAR, JAA, IAA', sum(aspar), sum(jaa), sum(iaa)
      WRITE(*,*) 'AU, FLJAU, FLJU', sum(au), sum(fljau), sum(flju)
    END IF
 
    DO ip = 1,nx
      u(ip) = max(0.d0,x(ip)*dble(iobpd(ith,ip)))
#ifdef W3_REF1
      IF (refpars(3).LT.0.5.AND.iobpd(ith,ip).EQ.0.AND.iobpa(ip).EQ.0) u(ip)= ac(ip) ! restores reflected boundary values
#endif
    END DO
    !
    ! update spectrum
    ac = real(u)
    !
    ! 4 Update boundaries: performs interpolation in time as done in rect grids (e.g. w3pro1md.ftn)
    !
    IF ( flbpi ) THEN
      rd1=rd10
      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
      !
      ! Time interpolation as done in rect grids
      !
      DO ibi=1, nbi
        ip    = mapsf(isbpi(ibi),1)
        ac(ip) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
             *iobpa(ip)*iobpd(ith,ip) / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
      END DO
    END IF
 
    !      CALL EXTCDE ( 99 )
    !/
    !/ End of W3XYPFSNIMP------------------------------------------------- /
    !/

References w3odatmd::bbpi0, w3odatmd::bbpin, bcgstab(), w3gdatmd::ccon, w3adatmd::cg, w3gdatmd::clats, w3timemd::dsec21(), w3odatmd::flbpi, w3gdatmd::iaa, w3gdatmd::ie_cell, w3gdatmd::ien, ilu0(), w3gdatmd::iobp, w3gdatmd::iobpa, w3gdatmd::iobpd, w3odatmd::isbpi, w3adatmd::iter, w3gdatmd::jaa, w3gdatmd::mapsf, w3odatmd::nbi, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nnz, w3gdatmd::nth, w3gdatmd::ntri, w3gdatmd::nx, w3gdatmd::pos_cell, w3gdatmd::posi, w3gdatmd::refpars, runrc(), w3gdatmd::si, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, w3gdatmd::tria, and w3gdatmd::trigp.

Referenced by w3xypug().

w3xypfspsi2()

subroutine w3profsmd::w3xypfspsi2	(	integer, intent(in)	ISP,
		real, dimension(:,:), intent(in)	C,
		logical, intent(in)	LCALC,
		real, intent(in)	RD10,
		real, intent(in)	RD20,
		real, intent(in)	DT,
		double precision, dimension(:), intent(inout)	AC
	)

Definition at line 724 of file w3profsmd.F90.

 
    !/
    !/
    !/                  +-----------------------------------+
    !/                  |  WWIII Version of the WWM FS Code |
    !/                  |  by Aron Roland                   |
    !/                  |  for use in WWIII                 |
    !/                  |  GPL License                      |
    !/                  |  Last update :        19-Dec-2007 |
    !/                  +-----------------------------------+
    !/
    !  1. Purpose :
    !     Advection of a scalar in a arbitary velocity field on unstructured meshes
    !     for the conservative hyperbolic equation N,t + (c*N),xy = 0 in spatial space
    !     This is the standard explicit N-Scheme from Roe as formulated in Abgrall
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Subroutine tracing (!/S switch)
    !
    !  5. Called by :
    !
    !       W3XYPUG   Routine for advection on unstructured grid
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE w3gdatmd, ONLY : nk, nth, ntri, nx, ccon, ie_cell,pos_cell, si, &
         ien, trigp, clats, mapsf, iobpa, iobpd, iobp, nnz, iobdp
#ifdef W3_REF1
    USE w3gdatmd, ONLY :  refpars
#endif
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, iter
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin, isbpi, bbpi0, bbpin
    USE w3timemd, ONLY: dsec21
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    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(:,:)              ! Velocity field in it's X- and Y- Components,
    DOUBLE PRECISION,INTENT(INOUT) :: AC(:)         ! 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
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
 
    real*8,  PARAMETER :: onesixth  = 1.0d0/6.0d0
    real*8,  PARAMETER :: thr8      = tiny(1.0d0)
    REAL,    PARAMETER :: THR       = tiny(1.0)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    !
    ! local integer
    !
    INTEGER :: IP, IE, IT, I1, I2, I3, ITH, IK
    INTEGER :: IBI, NI(3)
    !
    ! local real
    !
    REAL    :: RD1, RD2
    !:
    ! local double
    !
    real*8  :: utilde, boundary_forcing
    real*8  :: ft, cflxy
    real*8  :: fl11, fl12, fl21, fl22, fl31, fl32
    real*8  :: fl111, fl112, fl211, fl212, fl311, fl312
    real*8  :: dtsi(nx), u(nx)
    real*8  :: dtmaxgl, dtmaxexp, rest
    real*8  :: lambda(2), ktmp(3), tmp(3)
    real*8  :: theta_l(3), bet1(3), betahat(3)
    real*8  :: kelem(3,ntri), flall(3,ntri)
    real*8  :: kksum(nx), st(nx)
    real*8  :: nm(ntri)
 
#ifdef W3_S
    CALL strace (ient, 'W3XYPFSN')
#endif
 
    ! 1. initialisation
 
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    dtmaxgl = dble(10.e10)
    !
    !2       Propagation
    !2.a     Calculate K-Values and contour based quantities ...
    !
    DO ie = 1, ntri ! I precacalculate this arrays below as I assume that current velocity  changes continusly ...
      i1 = trigp(1,ie) ! Index of the Element Nodes (TRIGP)
      i2 = trigp(2,ie)
      i3 = trigp(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) * ien(ie,1) + lambda(2) * ien(ie,2) ! K-Values - so called Flux Jacobians
      kelem(2,ie) = lambda(1) * ien(ie,3) + lambda(2) * ien(ie,4)
      kelem(3,ie) = lambda(1) * ien(ie,5) + lambda(2) * ien(ie,6)
      ktmp        = kelem(:,ie) ! Copy
      nm(ie)      = - 1.d0/min(-thr8,sum(min(0.d0,ktmp))) ! N-Values
      kelem(:,ie) = max(0.d0,ktmp)
      fl11  = c(i2,1) * ien(ie,1) + c(i2,2) * ien(ie,2) ! Weights for Simpson Integration
      fl12  = c(i3,1) * ien(ie,1) + c(i3,2) * ien(ie,2)
      fl21  = c(i3,1) * ien(ie,3) + c(i3,2) * ien(ie,4)
      fl22  = c(i1,1) * ien(ie,3) + c(i1,2) * ien(ie,4)
      fl31  = c(i1,1) * ien(ie,5) + c(i1,2) * ien(ie,6)
      fl32  = c(i2,1) * ien(ie,5) + c(i2,2) * ien(ie,6)
      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 ! NTRI
 
    IF (lcalc) THEN ! 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
      kksum = 0.d0
      DO ie = 1, ntri
        ni = trigp(:,ie)
        kksum(ni) = kksum(ni) + kelem(:,ie)
      END DO ! IE
      dtmaxexp = 1e10 ! initialize to large number
      DO ip = 1, nx
        dtmaxexp = si(ip)/max(dble(10.e-10),kksum(ip)*iobdp(ip))
        dtmaxgl  = min( dtmaxgl, dtmaxexp)
      END DO ! IP
      cflxy = dble(dt)/dtmaxgl
      rest  = abs(mod(cflxy,1.0d0))
      IF (rest .LT. thr8) THEN
        iter(ik,ith) = abs(nint(cflxy))
      ELSE IF (rest .GT. thr8 .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
 
    DO ip = 1, nx
      dtsi(ip) = dble(dt)/dble(iter(ik,ith))/si(ip) ! Some precalculations for the time integration.
    END DO
 
    DO it = 1, iter(ik,ith)
      u = ac
 
      st = 0.d0
 
      DO ie = 1, ntri
        ni   = trigp(:,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(0.d0,min(betahat(1),1.d0-betahat(2),1.d0))
            bet1(2)       = max(0.d0,min(betahat(2),1.d0-betahat(3),1.d0))
            bet1(3)       = max(0.d0,min(betahat(3),1.d0-betahat(1),1.d0))
            theta_l(:) = ft * bet1
          END IF
        ELSE
          theta_l(:) = 0.d0
        END IF
        st(ni) = st(ni) + theta_l ! the 2nd term are the theta values of each node ...
      END DO
 
      DO ip = 1, nx
        u(ip) = max(0.d0,u(ip)-dtsi(ip)*st(ip)*(1-iobpa(ip)))*dble(iobpd(ith,ip))
#ifdef W3_REF1
        IF (refpars(3).LT.0.5.AND.iobpd(ith,ip).EQ.0.AND.iobpa(ip).EQ.0) u(ip)= ac(ip) ! restores reflected boundary values
#endif
      END DO
 
      ! update spectrum
      ac = u
      !
      ! 4 Update boundaries: performs interpolation in time as done in rect grids (e.g. w3pro1md.ftn)
      !
      IF ( flbpi ) THEN
        !
        ! 4.1 In this case the boundary is read from the nest.ww3 file
        !
        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
        !
        ! Overwrites only the incoming energy ( IOBPD(ITH,IP) = 0)
        !
        DO ibi=1, nbi
          ip = mapsf(isbpi(ibi),1)
          ac(ip) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
               / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
        END DO
 
      ENDIF
      !
    END DO ! End of loop on time steps
    !      CALL EXTCDE ( 99 )
    !/
    !/ End of W3XYPFSN ----------------------------------------------------- /
    !/

References w3odatmd::bbpi0, w3odatmd::bbpin, w3gdatmd::ccon, w3adatmd::cg, w3gdatmd::clats, w3timemd::dsec21(), w3odatmd::flbpi, w3gdatmd::ie_cell, w3gdatmd::ien, w3gdatmd::iobdp, w3gdatmd::iobp, w3gdatmd::iobpa, w3gdatmd::iobpd, w3odatmd::isbpi, w3adatmd::iter, w3gdatmd::mapsf, w3odatmd::nbi, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nnz, w3gdatmd::nth, w3gdatmd::ntri, w3gdatmd::nx, w3gdatmd::pos_cell, w3gdatmd::refpars, w3gdatmd::si, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, and w3gdatmd::trigp.

Referenced by w3xypug().

w3xypug()

subroutine w3profsmd::w3xypug	(	integer, intent(in)	ISP,
		real, intent(in)	FACX,
		real, intent(in)	FACY,
		real, intent(in)	DTG,
		real, dimension(1-ny:ny*(nx+2)), intent(inout)	VQ,
		real, intent(in)	VGX,
		real, intent(in)	VGY,
		logical, intent(in)	LCALC
	)

Definition at line 63 of file w3profsmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            Aron Roland            |
    !/                  |                        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 )
    !/
    !  1. Purpose :
    !
    !     Propagation in physical space for a given spectral component.
    !     Gives the choice of scheme on unstructured grid
    !
    !  2. Method :
    !
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISP     Int.   I   Number of spectral bin (IK-1)*NTH+ITH
    !       FACX/Y  Real   I   Factor in propagation velocity.
    !                          ( 1 or 0 * DT / DX )
    !       DTG     Real   I   Total time step.
    !       VQ      R.A.  I/O  Field to propagate.
    !       VGX/Y   Real   I   Speed of grid.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       VCFL0X  R.A.  Local courant numbers for absolute group vel.
    !                     using local X-grid step.
    !       VCFL0Y  R.A.  Id. in Y.
    !     ----------------------------------------------------------------
    !
    !  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, nsea, mapsf, mapfs, dtcfl, clats,   &
         flcx, flcy, nk, nth, dth, xfr,              &
         ecos, esin, sig,  pfmove,ien,               &
         ntri, trigp, ccon ,                         &
         ie_cell, pos_cell, iobp, iobpd, iobdp,      &
         fsn, fspsi, fsfct, fsnimp, gtype, ungtype
 
    USE w3wdatmd, ONLY: time
    USE w3odatmd, ONLY: tbpi0, tbpin, flbpi
    USE w3adatmd, ONLY: cg, cx, cy, atrnx, atrny, itime, cflxymax, dw
    USE w3idatmd, ONLY: flcur
    !      USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPI0, TBPIN,       &
    !                          ISBPI, BBPI0, BBPIN
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
 
    IMPLICIT NONE
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISP
    REAL, INTENT(IN)        :: FACX, FACY, DTG, VGX, VGY
    REAL, INTENT(INOUT)     :: VQ(1-NY:NY*(NX+2))
    LOGICAL, INTENT(IN)     :: LCALC
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, ISEA, IXY
    INTEGER                 :: IX
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CCOS, CSIN, CCURX, CCURY
    REAL                    :: C(NX,2)
    REAL                    :: RD1, RD2
    !/
    !/ Automatic work arrays
    !/
    REAL                    :: VLCFLX((NX+1)*NY), VLCFLY(NX*NY)
    DOUBLE PRECISION        :: AC(NX)
    !/ ------------------------------------------------------------------- /
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Set constants
    !
 
#ifdef W3_S
    CALL strace (ient, 'W3XYPUG')
#endif
    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.
    !
    ! 1.c Set depth
    !
    CALL setdepth
    !
    !
    ! 2.  Calculate velocities ---------------- *
    !
    DO isea = 1, nsea
      ixy         = mapsf(isea,3)
      vq(ixy)     = vq(ixy) / cg(ik,isea) * clats(isea)
      vlcflx(ixy) = ccos * cg(ik,isea) / clats(isea)
      vlcfly(ixy) = csin * cg(ik,isea)
#ifdef W3_MGP
      vlcflx(ixy) = vlcflx(ixy) - ccurx*vgx/clats(isea)
      vlcfly(ixy) = vlcfly(ixy) - ccury*vgy
#endif
    END DO
 
    IF ( flcur ) THEN
      DO isea=1, nsea
        ixy =  mapsf(isea,3)
        !
        ! Currents are not included on coastal boundaries (IOBP(IXY).EQ.0)
        !
        IF (iobp(ixy) .EQ. 1) THEN
          vlcflx(ixy) = vlcflx(ixy) + ccurx*cx(isea)/clats(isea)
          vlcfly(ixy) = vlcfly(ixy) + ccury*cy(isea)
        END IF
      END DO
    END IF
 
    !
    ! 3. initialize fluctuation splitting arrays ( to fit with WWM notations)
    !
    ac(1:nx)   = dble(vq(1:nx)) * iobdp(1:nx)
    c(1:nx,1)  = vlcflx(1:nx)   * iobdp(1:nx)
    c(1:nx,2)  = vlcfly(1:nx)   * iobdp(1:nx)
 
    !
    ! 4. Prepares boundary update
    !
    IF ( flbpi ) THEN
      rd1    = dsec21( tbpi0, time )
      rd2    = dsec21( tbpi0, tbpin )
    ELSE
      rd1=1.
      rd2=0.
    END IF
    !
    ! 4. propagate using the selected scheme
    !
    IF (fsn) THEN
      CALL w3xypfsn2 (isp, c, lcalc, rd1, rd2, dtg, ac)
    ELSE IF (fspsi) THEN
      CALL w3xypfspsi2 (isp, c, lcalc, rd1, rd2, dtg, ac)
    ELSE IF (fsfct) THEN
      CALL w3xypfsfct2 (isp, c, lcalc, rd1, rd2, dtg, ac)
    ELSE IF (fsnimp) THEN
      CALL w3xypfsnimp(isp, c, lcalc, rd1, rd2, dtg, ac)
    ENDIF
    !
    DO ix=1,nx
      isea=mapfs(1,ix)
      vq(ix)=real(ac(ix))
    ENDDO
 
    ! 6.  Store results in VQ in proper format --------------------------- *
    !
    DO isea=1, nsea
      ixy   = mapsf(isea,3)
      vq(ixy) =  max( 0. , cg(ik,isea)/clats(isea)*vq(ixy) )
    END DO

References w3adatmd::atrnx, w3adatmd::atrny, w3gdatmd::ccon, w3adatmd::cflxymax, w3adatmd::cg, w3gdatmd::clats, w3adatmd::cx, w3adatmd::cy, w3timemd::dsec21(), w3gdatmd::dtcfl, w3gdatmd::dth, w3adatmd::dw, w3gdatmd::ecos, w3gdatmd::esin, w3odatmd::flbpi, w3idatmd::flcur, w3gdatmd::flcx, w3gdatmd::flcy, w3gdatmd::fsfct, w3gdatmd::fsn, w3gdatmd::fsnimp, w3gdatmd::fspsi, w3gdatmd::gtype, w3gdatmd::ie_cell, w3gdatmd::ien, w3gdatmd::iobdp, w3gdatmd::iobp, w3gdatmd::iobpd, w3adatmd::itime, w3gdatmd::mapfs, w3gdatmd::mapsf, w3gdatmd::nk, w3gdatmd::nsea, w3gdatmd::nth, w3gdatmd::ntri, w3gdatmd::nx, w3gdatmd::ny, w3gdatmd::pfmove, w3gdatmd::pos_cell, setdepth(), w3gdatmd::sig, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, w3gdatmd::trigp, w3gdatmd::ungtype, w3xypfsfct2(), w3xypfsn2(), w3xypfsnimp(), w3xypfspsi2(), and w3gdatmd::xfr.

Referenced by w3wavemd::w3wave().