text:
module stream_module

  use amrex_fort_module, only : amrex_real, dim=>amrex_spacedim

  implicit none

  public

contains

  subroutine vtrace(T, T_lo, T_hi, nT, loc, loc_lo, loc_hi, nl,&
     &     ids, n_ids, g, g_lo, g_hi, computeVec, strm, strm_lo, strm_hi,&
     &     ncs, dx, plo, phi, hRK, errFlag) bind(C,name="vtrace")
    implicit none
    integer, intent(in) ::  nT, nl, computeVec, ncs, n_ids
    integer, intent(in) :: T_lo(3),T_hi(3),g_lo(3),g_hi(3),loc_lo(3),loc_hi(3),strm_lo(3),strm_hi(3)
    real(amrex_real), intent(in) :: T(T_lo(1):T_hi(1),T_lo(2):T_hi(2),T_lo(3):T_hi(3),nT)
    real(amrex_real), intent(inout) :: g(g_lo(1):g_hi(1),g_lo(2):g_hi(2),g_lo(3):g_hi(3),3)
    real(amrex_real), intent(in) :: loc(loc_lo(1):loc_hi(1),loc_lo(2):loc_hi(2),loc_lo(3):loc_hi(3),nl)
    real(amrex_real), intent(inout) :: strm(strm_lo(1):strm_hi(1),strm_lo(2):strm_hi(2),strm_lo(3):strm_hi(3),ncs)
    integer, intent(in) :: ids(0:n_ids-1)
    real(amrex_real), intent(in) :: dx(3), plo(3), phi(3), hRK

    integer :: i,j,k,n,m
    real(amrex_real) :: x(3),xp(3),xm(3)

    logical :: ok
    integer :: myproc, errFlag

    call bl_pd_myproc(myproc)

    ! Compute gradient field assuming grow cells valid
    if (computeVec .eq. 1) then
       do k=g_lo(3),g_hi(3)
          do j=g_lo(2),g_hi(2)
             do i=g_lo(1),g_hi(1)
                  g(i,j,k,1) = T(i+1,j,k,1) - T(i-1,j,k,1)
                  g(i,j,k,2) = T(i,j+1,k,1) - T(i,j-1,k,1)
                  if (dim.eq.3) then
                     g(i,j,k,3) = T(i,j,k+1,1) - T(i,j,k-1,1)
                  endif
             enddo
          enddo
       enddo
    endif

    ! Loop over list of nodes in this box and do work
    !  loc contains all the nodes, ids is a list of those in this box
    !  strm holds path trace data: i=cnt of node in this box, j=index away from iso, k=unused, m=path comp
    !  Note that the ids come from the inside_nodes structure which is 1-based, but the
    !  indexing into strm is boxlib 0-based, subtract 1 when setting j here.
    errFlag = 0
    do i=0,n_ids-1
       j = ids(i) - 1
       x(1:dim) = loc(j,loc_lo(2),loc_lo(3),1:dim)
       strm(i,0,strm_lo(3),1:dim) = x(1:dim)
       do m=1,nT
          call ntrpv(x,T(T_lo(1),T_lo(2),T_lo(3),m),T_lo,T_hi,dx,plo,phi,&
               strm(i,0,strm_lo(3),dim+m),1,ok)
          if (ok .eqv. .false.) then
             errFlag = 1
             return
          endif
       enddo

       ! Integrate with RK4, interpolate input data at new positions
       xp(1:dim) = x(1:dim)
       xm(1:dim) = x(1:dim)
       do n=-1,strm_lo(2),-1
          strm(i,n,strm_lo(3),1:dim) = xm(1:dim)
          call RK4(xm,-hRK,g,g_lo,g_hi,dx,plo,phi,ok)
          if (ok) then
             strm(i,n,strm_lo(3),1:dim) = xm(1:dim)
          else
             errFlag = 2
          endif
          do m=1,nT
             call ntrpv(xm,T(T_lo(1),T_lo(2),T_lo(3),m),T_lo,T_hi,dx,plo,phi,&
                  strm(i,n,strm_lo(3),dim+m),1,ok)
             if (ok .eqv. .false.) then
                strm(i,n,strm_lo(3),dim+m) = strm(i,n+1,strm_lo(3),dim+m)
             endif
          enddo
       enddo

       do n=1,strm_hi(2)
          strm(i,n,strm_lo(3),1:dim) = xp(1:dim)
          call RK4(xp,+hRK,g,g_lo,g_hi,dx,plo,phi,ok)
          if (ok) then
             strm(i,n,strm_lo(3),1:dim) = xp(1:dim)
          else
             errFlag = 4
          endif
          do m=1,nT
             call ntrpv(xp,T(T_lo(1),T_lo(2),T_lo(3),m),T_lo,T_hi,dx,plo,phi,&
                  strm(i,n,strm_lo(3),dim+m),1,ok)
             if (ok .eqv. .false.) then
                strm(i,n,strm_lo(3),dim+m) = strm(i,n-1,strm_lo(3),dim+m)
             endif
          enddo
       enddo
    enddo
  end subroutine vtrace

  subroutine IsOK(x,plo,phi,ok)
    implicit none
    real(amrex_real), intent(in) :: x(3),plo(3),phi(3)
    logical, intent(out) :: ok
    integer :: i
    ok = .true.
    do i=1,dim
       if (x(i).lt.plo(i) .or. x(i).gt.phi(i)) ok = .false.
    enddo
  end subroutine IsOK

  subroutine RK4(x,h,g,g_lo,g_hi,dx,plo,phi,ok)
    implicit none
    real(amrex_real), intent(inout) :: x(3)
    real(amrex_real), intent(in) :: h,dx(3),plo(3),phi(3)
    integer, intent(in) :: g_lo(3),g_hi(3)
    real(amrex_real), intent(in) :: g(g_lo(1):g_hi(1),g_lo(2):g_hi(2),g_lo(3):g_hi(3),3)

    real(amrex_real) :: xx(3),vec(3),k1(3),k2(3),k3(3),k4(3)
    logical ok

    xx(1:dim) = x(1:dim)
    call ntrpv(xx,g,g_lo,g_hi,dx,plo,phi,vec,dim,ok)
    if ( ok .eqv. .false. ) return
    call vnrml(vec)

    k1(1:dim) = vec(1:dim)*h
    xx(1:dim) = x(1:dim) + k1(1:dim)*.5d0
    call ntrpv(xx,g,g_lo,g_hi,dx,plo,phi,vec,dim,ok)
    if ( ok .eqv. .false. ) return
    call vnrml(vec)

    k2(1:dim) = vec(1:dim)*h
    xx(1:dim) = x(1:dim) + k2(1:dim)*.5d0
    call ntrpv(xx,g,g_lo,g_hi,dx,plo,phi,vec,dim,ok)
    if ( ok .eqv. .false. ) return
    call vnrml(vec)

    k3(1:dim) = vec(1:dim)*h
    xx(1:dim) = x(1:dim) + k3(1:dim)
    call ntrpv(xx,g,g_lo,g_hi,dx,plo,phi,vec,dim,ok)
    if ( ok .eqv. .false. ) return
    call vnrml(vec)

    k4(1:dim) = vec(1:dim)*h
    x(1:dim) = x(1:dim) + (k1(1:dim)+k4(1:dim))/6.d0 + (k2(1:dim)+k3(1:dim))/3.d0
  end subroutine RK4

  subroutine ntrpv(x,g,g_lo,g_hi,dx,plo,phi,u,nc,ok)
    implicit none
    integer, intent(in) :: nc
    real(amrex_real), intent(in) :: x(3), dx(3), plo(3), phi(3)
    real(amrex_real), intent(inout) :: u(nc)
    integer, intent(in) :: g_lo(3), g_hi(3)
    real(amrex_real), intent(in) :: g(g_lo(1):g_hi(1),g_lo(2):g_hi(2),g_lo(3):g_hi(3),nc)
    logical, intent(inout) :: ok

    integer :: b(3), i
    real(amrex_real) :: n(3), tmp

    call IsOK(x,plo,phi,ok)
    if (ok .eqv. .false.) return

    b = 0
    n = 0
    do i=1,dim
       tmp = (x(i) - plo(i)) / dx(i) - 0.5d0
       b(i) = FLOOR( tmp )
       n(i) = ( x(i) - ( (b(i)+0.5d0)*dx(i) + plo(i) ) )/dx(i)
       n(i) = MAX(0.d0,MIN(1.d0,n(i)))
    enddo

    ok = .true.
    do i=1,dim
       if (b(i).lt.g_lo(i) .or. b(i).gt.g_hi(i)-1) ok = .false.
    enddo
    if (.not. ok) then
       !print *,'b:',b
       !print *,'DIMS:',g_lo,g_hi
       return
    endif

    do i=1,nc
       if (dim.eq.2) then
          u(i) = &
               +       n(1)  *    n(2)   * g(b(1)+1,b(2)+1,g_lo(3),i) &
               +       n(1)  *(1.d0-n(2))* g(b(1)+1,b(2)  ,g_lo(3),i) &
               +  (1.d0-n(1))*    n(2)   * g(b(1)  ,b(2)+1,g_lo(3),i) &
               +  (1.d0-n(1))*(1.d0-n(2))* g(b(1)  ,b(2)  ,g_lo(3),i)
       else
          u(i) = &
               +       n(1)  *    n(2)   *    n(3)    * g(b(1)+1,b(2)+1,b(3)+1,i) &
               +       n(1)  *(1.d0-n(2))*    n(3)    * g(b(1)+1,b(2)  ,b(3)+1,i) &
               +       n(1)  *    n(2)   *(1.d0-n(3)) * g(b(1)+1,b(2)+1,b(3)  ,i) &
               +       n(1)  *(1.d0-n(2))*(1.d0-n(3)) * g(b(1)+1,b(2)  ,b(3)  ,i) &
               +  (1.d0-n(1))*    n(2)   *    n(3)    * g(b(1)  ,b(2)+1,b(3)+1,i) &
               +  (1.d0-n(1))*(1.d0-n(2))*    n(3)    * g(b(1)  ,b(2)  ,b(3)+1,i) &
               +  (1.d0-n(1))*    n(2)   *(1.d0-n(3)) * g(b(1)  ,b(2)+1,b(3)  ,i) &
               +  (1.d0-n(1))*(1.d0-n(2))*(1.d0-n(3)) * g(b(1)  ,b(2)  ,b(3)  ,i)
       endif
    enddo
  end subroutine ntrpv

  subroutine vnrml(vec)
    implicit none
    real(amrex_real), intent(inout) :: vec(3)
    real(amrex_real) :: eps, sum
    parameter (eps=1.e-12)
    integer :: i

    sum = 0._amrex_real
    do i=1,dim
       sum = sum + vec(i)*vec(i)
    enddo
    if (sum .gt. eps) then
       sum = sqrt(sum)
       do i=1,dim
          vec(i) = vec(i) / sum
       enddo
    endif
  end subroutine vnrml

end module stream_module

max_stars_repo_name: AMReX-Combustion/PeleAnalysis

max_issues_repo_name: AMReX-Combustion/PeleAnalysis

max_forks_repo_name: AMReX-Combustion/PeleAnalysis


text:
PROGRAM test_comptonprofiles

USE, INTRINSIC :: ISO_C_BINDING
USE, INTRINSIC :: ISO_FORTRAN_ENV
USE :: xraylib
USE :: libtest
IMPLICIT NONE

TYPE(xrl_error), POINTER :: error => NULL()
REAL (C_DOUBLE) :: profile, profile1, profile2

! pz == 0.0
profile = ComptonProfile(26, 0.0_C_DOUBLE, error)
CALL assert(abs(profile - 7.060_C_DOUBLE) < 1E-6)
CALL assert(.NOT. ASSOCIATED(error))

profile = ComptonProfile_Partial(26, N1_SHELL, 0.0_C_DOUBLE, error)
CALL assert(abs(profile - 1.550_C_DOUBLE) < 1E-6)
CALL assert(.NOT. ASSOCIATED(error))

profile1 = ComptonProfile_Partial(26, L2_SHELL, 0.0_C_DOUBLE, error)
profile2 = ComptonProfile_Partial(26, L3_SHELL, 0.0_C_DOUBLE, error)
CALL assert(abs(profile1 - profile2) < 1E-6)
CALL assert(abs(profile1 - 0.065_C_DOUBLE) < 1E-6)
CALL assert(.NOT. ASSOCIATED(error))

! pz == 100.0
profile = ComptonProfile(26, 100.0_C_DOUBLE, error)
CALL assert(abs(profile - 1.800E-05_C_DOUBLE) < 1E-8)
CALL assert(.NOT. ASSOCIATED(error))

profile = ComptonProfile_Partial(26, N1_SHELL, 100.0_C_DOUBLE, error)
CALL assert(abs(profile - 5.100E-09_C_DOUBLE) < 1E-12)
CALL assert(.NOT. ASSOCIATED(error))

profile1 = ComptonProfile_Partial(26, L2_SHELL, 100.0_C_DOUBLE, error)
profile2 = ComptonProfile_Partial(26, L3_SHELL, 100.0_C_DOUBLE, error)
CALL assert(abs(profile1 - profile2) < 1E-10)
CALL assert(abs(profile1 - 1.100E-08_C_DOUBLE) < 1E-10)
CALL assert(.NOT. ASSOCIATED(error))

! pz == 50.0 -> interpolated!
profile = ComptonProfile(26, 50.0_C_DOUBLE, error)
CALL assert(abs(profile - 0.0006843950273082384_C_DOUBLE) < 1E-8)
CALL assert(.NOT. ASSOCIATED(error))

profile = ComptonProfile_Partial(26, N1_SHELL, 50.0_C_DOUBLE, error)
CALL assert(abs(profile - 2.4322755767709126e-07_C_DOUBLE) < 1E-10)
CALL assert(.NOT. ASSOCIATED(error))

profile1 = ComptonProfile_Partial(26, L2_SHELL, 50.0_C_DOUBLE, error)
profile2 = ComptonProfile_Partial(26, L3_SHELL, 50.0_C_DOUBLE, error)
CALL assert(abs(profile1 - profile2) < 1E-10)
CALL assert(abs(profile1 - 2.026953933016568e-06_C_DOUBLE) < 1E-10)
CALL assert(.NOT. ASSOCIATED(error))

! bad input
profile = ComptonProfile(0, 0.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile(0, 0.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile(102, 0.0_C_DOUBLE, error)
CALL assert(.NOT. ASSOCIATED(error))

profile = ComptonProfile(103, 0.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile(103, 0.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile(26, -1.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile(26, -1.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile_Partial(0, K_SHELL, 0.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile_Partial(0, K_SHELL, 0.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile_Partial(102, K_SHELL, 0.0_C_DOUBLE, error)
CALL assert(.NOT. ASSOCIATED(error))

profile = ComptonProfile_Partial(103, K_SHELL, 0.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile_Partial(103, K_SHELL, 0.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile_Partial(26, K_SHELL, -1.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile_Partial(26, K_SHELL, -1.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile_Partial(26, -1, 0.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile_Partial(26, -1, 0.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

profile = ComptonProfile_Partial(26, N2_SHELL, 0.0_C_DOUBLE)
CALL assert(profile == 0.0)
profile = ComptonProfile_Partial(26, N2_SHELL, 0.0_C_DOUBLE, error)
CALL assert(profile == 0.0)
CALL assert(ASSOCIATED(error))
CALL assert(error%code == XRL_ERROR_INVALID_ARGUMENT)
DEALLOCATE(error)

ENDPROGRAM test_comptonprofiles

max_stars_repo_name: mschollmeier/xraylib

max_issues_repo_name: mschollmeier/xraylib

max_forks_repo_name: mschollmeier/xraylib


text:
      PROGRAM W02
      INTEGER I, N
      REAL T(10)

      I = 1
      N = 10

      DO WHILE (I.LT.N)
         PRINT *, T(I)
         I = I+1
      ENDDO

      PRINT *, I

      END

ext: f


subroutine getvort()
!Calculates the vorticity on the edge of the object.
!The vorticity is taken to be a constant on the surface:
! vortcoeff on the top edge, and -vortcoeff on the bottom edge
    use vars
    use parallel

    implicit none

    integer :: i, j

    call getv() !get the velocity

    do j=1,ny
        do i=1,nx
            if (boundary(i,j) .ne. 0) vort(i,j) = real(boundary(i,j))*vortcoeff
        enddo
    enddo

    vort(nx+1,:) = 0.

    !transfer halo values for the vorticity
    call haloswap(vort)

end subroutine


text:
	SUBROUTINE GDPINP  ( gdfile, gdatim, gvcord, gfunc, gpoint,
     +			     ptype,  xaxis,  yaxis, border, line,
     +			     marker, title,  clear, device, scale,
     +			     panel,  gvect,  wind,  refvec, winpos,
     +                       filter, text, output, thtaln, thteln,
     +                       mixrln, iret )
C************************************************************************
C* GDPINP								*
C*									*
C* This subroutine gets the input parameters for GDPROF.		*
C*									*
C* GDPINP  ( GDFILE, GDATIM, GVCORD, GFUNC, GPOINT, PTYPE, XAXIS, 	*
C*           YAXIS,  BORDER, LINE,  MARKER, TITLE,  CLEAR, DEVICE,	*
C*           SCALE,  PANEL,  GVECT, WIND, REFVEC,  WINPOS, FILTER,	*
C*           TEXT, OUTPUT, THTALN, THTELN, MIXRLN, IRET )		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 2/84						*
C* M. desJardins/GSFC	 9/88	GEMPAK4					*
C* K. Brill/GSC         12/89   Added wind input			*
C* K. Brill/GSC          1/90   Added TEXT				*
C* K. Brill/GSC          5/90   Added OUTPUT				*
C* K. Brill/GSC          5/90   Added THTALN, THTELN, MIXRLN		*
C* J. Nielsen/TAMU	11/91	Made FILTER a string			*
C* L. Sager/NMC		 7/93	Added refvec				*
C************************************************************************
	CHARACTER*(*)	gdfile, gdatim, gvcord, gfunc, gpoint, ptype,
     +			xaxis,  yaxis, border, line, marker, title,
     +			device, scale, panel, gvect, wind, refvec,
     +                  winpos, text, output, thtaln, thteln, mixrln,
     +			filter
	LOGICAL		clear
C------------------------------------------------------------------------
	CALL IP_STR  ( 'GDFILE',  gdfile,  ier1 )
	CALL IP_STR  ( 'GDATTIM', gdatim,  ier2 )
	CALL IP_STR  ( 'GVCORD',  gvcord,  ier3 )
	CALL IP_STR  ( 'GFUNC',   gfunc,   ier4 )
	CALL IP_STR  ( 'GPOINT',  gpoint,  ier5 )
	CALL IP_STR  ( 'PTYPE',   ptype,   ier6 )
	CALL IP_STR  ( 'XAXIS',   xaxis,   ier7 )
	CALL IP_STR  ( 'YAXIS',   yaxis,   ier8 )
	CALL IP_STR  ( 'BORDER',  border,  ier9 )
	CALL IP_STR  ( 'LINE',    line,    ier10 )
	CALL IP_STR  ( 'MARKER',  marker,  ier11 )
	CALL IP_STR  ( 'TITLE',   title,   ier12 )
	CALL IP_LOG  ( 'CLEAR',   clear,   ier13 )
	CALL IP_STR  ( 'DEVICE',  device,  ier14 )
	CALL IP_STR  ( 'SCALE',   scale,   ier15 )
	CALL IP_STR  ( 'PANEL',   panel,   ier16 )
	CALL IP_STR  ( 'GVECT',   gvect,   ier17 )
	CALL IP_STR  ( 'WIND',    wind,    ier18 )
	CALL IP_STR  ( 'REFVEC',  refvec,  ier19 )
        CALL IP_STR  ( 'WINPOS',  winpos,  ier20 )
	CALL IP_STR  ( 'FILTER',  filter,  ier21 )
	CALL IP_STR  ( 'TEXT',    text,    ier22 )
	CALL IP_STR  ( 'OUTPUT',  output,  ier23 )
	CALL IP_STR  ( 'THTALN',  thtaln,  ier24 )
	CALL IP_STR  ( 'THTELN',  thteln,  ier25 )
	CALL IP_STR  ( 'MIXRLN',  mixrln,  ier26 )
C*
	iret = ier1 + ier2 + ier3 + ier4 + ier5 + ier6 + ier7 + ier8 +
     +	       ier9 + ier10 + ier11 + ier12 + ier13 + ier14 + ier15 +
     +	       ier16 + ier17 + ier18 + ier19 + ier20 + ier21 + ier22 +
     +         ier23 + ier24 + ier25 + ier26
	IF  ( iret .ne. 0 )  iret = -2
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      MODULE MYSTRAN_Version

! Data used for writing the MYSTRAN version and related data to the F06 file

      USE PENTIUM_II_KIND, ONLY      :  BYTE, LONG, DOUBLE

      IMPLICIT NONE

      SAVE

      CHARACTER(256*BYTE)            :: MYSTRAN_COMMENT  = '*** Please report any problems to the author at dbcase29@gmail.com ***'
      CHARACTER(  8*BYTE), PARAMETER :: MYSTRAN_VER_NUM  = '13.3a'
      CHARACTER(  3*BYTE), PARAMETER :: MYSTRAN_VER_MONTH= 'Oct'
      CHARACTER(  2*BYTE), PARAMETER :: MYSTRAN_VER_DAY  = '31'
      CHARACTER(  4*BYTE), PARAMETER :: MYSTRAN_VER_YEAR = '2021'
      CHARACTER( 33*BYTE), PARAMETER :: MYSTRAN_AUTHOR   = 'MYSTRAN developed by Dr Bill Case'

      END MODULE MYSTRAN_Version


text:
Knit Bombing (also called graffiti knitting) has become increasingly popular in Davis and the effects of these random acts of knitting can be seen all over town.  The Davis Food Coops Knitting and Craft Circle http://www.youtube.com/watch?vj_oTGSzmGc&featureshare made a video of their knit bombings,

There are a couple of knit bombs near the coop on the South side of the parking lot. Looks good ~SD

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE GP1
C
C     GP1  BUILDS THE FOLLOWING DATA BLOCKS--
C       1. GRID POINT LIST (GPL)
C       2. EXTERNAL INTERNAL GRID POINT EQUIVALENCE TABLE (EQEXIN)
C       3. GRID POINT DEFINITION TABLE (GPDT)
C       4. COORDINATE SYSTEM TRANSFORMATION MATRICES (CSTM)
C       5. BASIC GRID POINT DEFINITION TABLE (BGPDT)
C       6. SCALAR INDEX LIST (SIL)
C
C     THE FOLLOWING CARDS ARE READ BY GP1--
C       1. GRID
C       2. CELASI, CDAMPI, CMASSI  (I=1,2,3,4)
C       3. SPOINT
C       4. SEQGP   (SEQEP IS PROCESSED IN DPD1)
C       5. CORDIJ  (I=1,2,  J=R,S,C)
C
C     IMPORTANT
C     =========
C     REVISED  7/89 BY G.CHAN/UNISYS, TO ALLOW GRID, SCALAR AND EXTRA
C     POINT EXTERNAL ID UP TO 8 DIGITS FOR ALL 32-BIT MACHINES
C     PREVIOUSLY, ID OF 2000000 IS THE UPPER LIMIT FOR IBM AND VAX
C
C     REVISED  8/89 BY G.CHAN/UNISYS, AS PART OF THE EFFORT TO ALLOW A
C     NASTRAN JOB TO EXCEED 65535 LIMIT.
C     NORMALLY, IF GRID POINTS OR SCALAR POINTS DO NOT HAVE VERY LARGE
C     EXTERNAL ID NUMBERS, THEIR ID NOS. ARE MULTIPLIED BY 1000, SO THAT
C     999 ADDITIONAL POINTS CAN SQUEEZE IN VIA SEQGP CARDS. (NOTE - A
C     7- OR 8-DIGIT ID NO., TIMES 1000, EXCEEDS A 32-BIT WORD COMPUTER
C     HARDWARE LIMIT). THIS MULTIPLY FACTOR IS NOW ADJUSTABLE, 1000,100,
C     OR 10, SO THAT ADDITIONAL DIGITS CAN BE USED FOR THE EXTERNAL GRID
C     OR SCALAR POINTS IN CASE THERE ARE LIMITTED SEQGP CARDS PRESENT.
C     THIS VARIABLE MULTIPLIER (10,100, OR 1000) IS ALSO RECORDED IN THE
C     3RD WORD OF THE HEADER RECORD OF THE GPL DATA BLOCK FOR LATER USE.
C     THE ACTUAL FACTOR OF THE MULTIPLIER IS ALSO MACHINE DEPENDENT.
C     UNIVAC, A 36-BIT MACHINE, CAN HAVE A MULTIPLIER OF 100 OR 1000.
C     OTHER 60- OR 64- BIT MACHINES, THE MULTIPLIER REMAINS AT 1000
C     IF THE MULTIPLIER IS 1000, THE SEQGP AND SEQEP CARDS, AS BEFORE,
C     CAN HAVE 4 SEQID LEVELS, SUCH AS XXX.X.X.X
C     IF THE MULTIPLIER IS 100, SEQGP AND SEQEP CARDS ARE LIMITED TO
C     3 SEQID LEVELS, XXX.X.X
C     FINALLY, IF MULTIPLIER IS 10, SEQGP AND SEQEP ARE LIMITED TO XXX.X
C
C     SPECIAL CONSIDERATION FOR THE AXISYM. AND HYDROELAS. PROBLEMS - 10
C     IS USED FOR THE MULTIPLIER, AND THEREFOR A ONE SEQID LEVEL IS
C     AVAILABLE. PREVIOUSLY, SEQGP CARDS WERE NOT USED IN AXISYM. AND
C     HYDROELAS. PROBLEMS, AND NO USER WARNING MESSAGE PRINTED
C
C     NO ADJUSTABLE MULTIPLY FACTOR FOR SUBRSTRUCTURING (MULT=1000,
C     SEE ALSO SGEN)
C
C     THE 65535 LIMITATION INVOLVES ONLY A SAMLL CHANGE IN STA 973
C
      EXTERNAL        RSHIFT
      INTEGER         RD,WRT,CLS,FILE,ELEM,AXIC,Z,SYSBUF,BUF1,BUF2,
     1                BUF3,GEOMP,GPL,EQEXIN,GPDT,CSTM,BGPDT,SIL,SCR1,
     2                SCR2,WRTREW,RDREW,A,SPOINT,FLAG,GRID,CLSREW,
     3                SEQGP,GPFL,CORD,CORDIJ,GP1AH,GEOM1,GEOM2,PTR,
     4                SOLV,SOLVP,SCALPT,TYPE,OFFSET,RSHIFT
      REAL            LENGTH
      DIMENSION       A(34),AA(34),AB(3),AC(3),AI(3),AJ(3),AK(3),AX(3),
     1                AR(3),SPOINT(2),GRID(2),SEQGP(2),CORDIJ(12),
     2                CORD(6),GP1AH(2),SCALPT(2),ZZ(1),MCB(7)
      CHARACTER*29    LVL1,LVL2
      CHARACTER       UFM*23,UWM*25,UIM*29
      COMMON /XMSSG / UFM,UWM,UIM
      COMMON /BLANK / LUSET,NOGPDT,NOCSTM
      COMMON /CONDAS/ PI,TWOPI,RADEG,DEGRA,S4PISQ
      COMMON /ZZZZZZ/ Z(1)
      COMMON /SYSTEM/ KSYSTM(100)
      COMMON /SETUP / NFILE(6),PTR
      COMMON /GPTA1 / NELEM,LASTX,INCRX,ELEM(1)
      COMMON /NAMES / RD,RDREW,WRT,WRTREW,CLSREW,CLS
      EQUIVALENCE     (KSYSTM( 1),SYSBUF), (KSYSTM( 2),IOUT ),
     1                (KSYSTM(24),ICFIAT), (KSYSTM(27),AXIC ),
     2                (KSYSTM(38),IAXIF ), (KSYSTM(40),NBPW ),
     3                (KSYSTM(56),ITHERM), (KSYSTM(69),ISUBS)
      EQUIVALENCE     (Z( 1),ZZ(1)), (A( 1),AA(1)), (A( 4),AB(1)),
     1                (A( 7),AC(1)), (A(10),AI(1)), (A(13),AJ(1)),
     2                (A(16),AK(1)), (A(19),AX(1)), (A(22),AR(1)),
     3                (NOCSTM,IFL ), (GEOMP,GEOM1), (MCB(2),KN  )
      EQUIVALENCE     (IGPDT,ICSDT)
      DATA    GEOM1 / 101/, GEOM2 / 102/,
     1        GPL   / 201/, EQEXIN/ 202/, GPDT  / 203/,
     2        CSTM  / 204/, BGPDT / 205/, SIL   / 206/,
     3        SCR1  / 301/, SCR2  / 302/
      DATA    GP1AH / 4HGP1 , 4H    /,
     1        CORD  / 6,6,6,13,13,13/,
     2        GRID  / 4501,45  /,
     3        SEQGP / 5301,53  /,
     4        CORDIJ/ 1701,17,1801,18,1901,19,2001,20,2101,21,2201,22/,
     5        SCALPT/ 5551,49  /
      DATA    MCB   / 7*0      /,
     1        LARGE / 100000000/,
     2        LVL1  / '3  I.E.  XXX.X.X.X TO XXX.X.X' /,
     3        LVL2  / '2  I.E.  XXX.X.X.X TO XXX.X  ' /
C
C     PERFORM GENERAL INITIALIZATION
C
      CALL DELSET
      NZ     = KORSZ(Z)
      BUF1   = NZ - SYSBUF - 2
      BUF2   = BUF1 - SYSBUF
      BUF3   = BUF2 - SYSBUF
      NOGO   = 0
      NOCSTM = 0
      NOGPDT =-1
      NOGMP1 = 1
      MAXA1  = 0
      MULT   = 1000
      AXI    = 0
      IF (AXIC.NE.0 .OR. IAXIF.NE.0) AXI = 1
      IF (AXI   .NE. 0) MULT = 10
      IF (ISUBS .NE. 0) MULT = 1000
      IMAX = LARGE
      IF (NBPW .EQ. 32) IMAX =  2147483
      IF (NBPW .EQ. 36) IMAX = 34359738
C         2147483=2**31/1000   34359738=2**35/1000
C
C     READ SCALAR ELEMENT CONNECTION CARDS (IF PRESENT).
C     EXTRACT SCALAR POINTS AND WRITE THEM ON SCR2.
C
      FILE = SCR2
      CALL OPEN (*1170,SCR2,Z(BUF2),WRTREW)
      NOSCLR= 0
      M8    =-8
      A(11) =-1
      DO 30 K = 12,16
   30 A(K) = 0
      CALL PRELOC (*80,Z(BUF1),GEOM2)
      I = 1
      DO 60 I = 1,LASTX,INCRX
      KK = ELEM(I+10)
      IF (KK .EQ. 0) GO TO 60
      CALL LOCATE (*60,Z(BUF1),ELEM(I+3),FLAG)
      NN = ELEM(I+5)
   40 CALL READ (*1180,*60,GEOM2,A,NN,0,FLAG)
      DO 50 K = 3,4
      IF (A(K).EQ.0 .OR. (KK.EQ.1 .AND. A(K+2).NE.0)) GO TO 50
      A(10)  = A(K)
      NOSCLR = 1
      CALL WRITE (SCR2,A(10),1,0)
   50 CONTINUE
      GO TO 40
   60 CONTINUE
C
C     COPY SCALAR POINTS DEFINED ON SPOINT CARDS (IF PRESENT) ONTO SCR2.
C
      CALL LOCATE (*80,Z(BUF1),SCALPT,FLAG)
      NOSCLR = 1
      CALL READ   (*1180,*70,GEOM2,Z,BUF2-1,1,N)
      CALL MESAGE (M8,0,GP1AH)
   70 CALL WRITE  (SCR2,Z,N,0)
C
C     CLOSE FILES. IF SCALAR POINTS PRESENT, SORT LIST.
C     THEN DISCARD DUPLICATES AND WRITE UNIQUE LIST ON SCR2.
C
   80 CALL WRITE (SCR2,0,0,1)
      CALL CLOSE (SCR2,CLSREW)
      CALL CLOSE (GEOM2,CLSREW)
      IF (NOSCLR .EQ. 0) GO TO 110
      NFILE(1) = GPDT
      NFILE(2) = BGPDT
      NFILE(3) = SIL
      CALL OPEN  (*1170,SCR2,Z(BUF1),RDREW)
      CALL SORTI (SCR2,0,1,1,Z,BUF1-1)
      CALL CLOSE (SCR2,CLSREW)
      FILE = NFILE(6)
      CALL OPEN (*1170,FILE,Z(BUF1),RDREW)
      CALL OPEN (*1170,SCR2,Z(BUF2),WRTREW)
      LAST = -1
   90 CALL READ (*1180,*100,FILE,A(10),1,0,FLAG)
      IF (A(10) .EQ. LAST) GO TO 90
      CALL WRITE (SCR2,A(10),1,0)
      LAST = A(10)
      GO TO 90
  100 CALL WRITE (SCR2,0,0,1)
      CALL CLOSE (SCR2,CLSREW)
      CALL CLOSE (FILE,CLSREW)
      CALL OPEN  (*1170,SCR2,Z(BUF3),RDREW)
C
C     READ GRID ENTRIES (IF PRESENT).
C     MERGE GRID AND SCALAR NOS.
C     CREATING LIST IN CORE OF EXTERNAL NO., MULT * EXTERNAL NO.
C     WRITE 7-WORD GRID AND SCALAR ENTRIES ON SCR1.
C
  110 A(1)  = LARGE
      A(10) = LARGE
      FILE  = SCR1
      IF (MAXA1 .EQ. 0) CALL OPEN (*1170,SCR1,Z(BUF2),WRTREW)
      I = -1
      NOGRID = 0
      IF (MAXA1 .EQ. 0) CALL PRELOC (*190,Z(BUF1),GEOM1)
      CALL LOCATE (*200,Z(BUF1),GRID,FLAG)
      NOGRID = 1
      CALL READ  (*1180,*1200,GEOM1,A,8,0,FLAG)
      CALL WRITE (SCR1,A,7,0)
  120 IF (NOSCLR .EQ. 0) GO TO 140
      CALL READ  (*1180,*1200,SCR2,A(10),1,0,FLAG)
      CALL WRITE (SCR1,A(10),7,0)
  130 IF (NOGRID .EQ. 0) GO TO 160
      IF (NOSCLR .EQ. 0) GO TO 140
      IF (A(1) -  A(10)) 140,1250,160
C
C     GRID NO. .LT. SCALAR NO.
C
  140 I = I + 2
      Z(I) = A(1)
C
C     GRID POINT EXTERNAL ID * MULT IS LIMITED TO COMPUTER MAXIMUM
C     INTEGER SIZE
C
      IF (A(1).LE.IMAX .OR. AXI.NE.0) GO TO 142
      IF (A(1) .GT. MAXA1) MAXA1 = A(1)
      GO TO 146
  142 Z(I+1) = MULT*A(1)
  146 CALL READ  (*1180,*150,GEOM1,A,8,0,FLAG)
      CALL WRITE (SCR1,A,7,0)
      GO TO 130
  150 NOGRID = 0
      A(1) = LARGE
      IF (NOSCLR .EQ. 0) GO TO 180
C
C     SCALAR NO. .LT. GRID NO.
C
  160 I = I + 2
      Z(I) = A(10)
C
C     SCALAR POINT EXTERNAL ID * MULT IS LIMITED TO COMPUTER MAXIMUM
C     INTEGER SIZE
C
      IF (A(10).LE.IMAX .OR. AXI.NE.0) GO TO 162
      IF (A(10) .GT. MAXA1) MAXA1 = A(10)
      GO TO 166
  162 Z(I+1) = MULT*A(10)
  166 CALL READ  (*1180,*170,SCR2,A(10),1,0,FLAG)
      CALL WRITE (SCR1,A(10),7,0)
      GO TO 130
  170 NOSCLR = 0
      A(10)  = LARGE
      IF (NOGRID .EQ. 0) GO TO 180
      GO TO 140
C
C     LIST COMPLETE ONLY IF MAXA1 .LE. ZERO
C
C     IF MAXA1 IS .GT. ZERO, SOME LARGE GRID OR SCALAR POINTS HAD BEEN
C     LEFT OUT IN LIST. MAXA1 IS THE LARGEST GRID OR SCALAR POINT
C     EXTERNAL ID.  RESET MULT AND REPEAT COMPILING LIST
C
  180 IF (MAXA1 .LE. 0) GO TO 185
      IF (ISUBS .NE. 0) GO TO 183
      CALL REWIND (SCR1)
      CALL REWIND (GEOM1)
      IF (NOSCLR .NE. 0) CALL REWIND (SCR2)
      MULT = 100
      IF (MAXA1 .GT. IMAX*10) MULT = 10
      IMAX  = (IMAX/MULT)*1000
      MAXA1 = -1
CWKBR CALL PAGE (-3)
      CALL PAGE2(-3)
      IF (MULT .EQ. 100) WRITE (IOUT,182) UWM,LVL1
      IF (MULT .EQ.  10) WRITE (IOUT,182) UWM,LVL2
  182 FORMAT (A25,' 2140A, DUE TO THE PRESENCE OF ONE OR MORE GRID OR ',
     1       'SCALAR POINTS WITH VERY LARGE EXTERNAL ID''S, THE SEQGP' ,
     2       /5X,'AND SEQEP CARDS, IF USED, ARE FORCED TO REDUCE FROM ',
     3       'ALLOWABLE 4 SEQID LEVELS TO ',A29,/)
      GO TO 110
C
  183 WRITE  (IOUT,184) UFM
  184 FORMAT (A23,' 2140B, EXTERNAL GRID OR SCALAR POINT ID TOO BIG')
      CALL MESAGE (-61,0,0)
C
  185 N    = I
      NEQEX= N
      N1   = N + 1
      N2   = N + 2
      IGPDT= N2
      ILIST= N2
      KN   = N1/2
      CALL CLOSE (SCR1,CLSREW)
      CALL CLOSE (SCR2,CLSREW)
      GO TO 210
C
C     NO GRID CARDS PRESENT-- TEST FOR ANY SCALAR PTS.
C
  190 NOGMP1 = 0
  200 IF (NOSCLR .EQ. 0) GO TO 980
      GO TO 120
C
C     READ THE SEQGP TABLE (IF PRESENT)
C     FOR EACH ENTRY, FIND MATCH IN THE SORTED EXTERNAL GRID POINTS
C     AND REPLACE SEQUENCE NO. WITH SEQGP NO.
C
  210 NOSEQ  = 0
      NOGPDT = 1
      IF (NOGMP1 .EQ. 0) GO TO 260
      ASSIGN 230 TO NDX
      SPOINT(2) = 0
      IERR = 1
      ASSIGN 220 TO NERR
      CALL LOCATE (*250,Z(BUF1),SEQGP,FLAG)
      NOSEQ = 1
      IFAIL = 0
 2010 CALL READ (*1180,*2020,GEOMP,Z(N2),BUF1-1,1,FLAG)
      IFAIL = IFAIL + 1
      GO TO 2010
 2020 IF (IFAIL .EQ. 0) GO TO 2060
      NWDS = (IFAIL-1)*(BUF1-1) + FLAG
      WRITE  (IOUT,2040) UFM,NWDS
 2040 FORMAT (A23,' 3135, UNABLE TO PROCESS SEQGP DATA IN SUBROUTINE ',
     1       'GP1 DUE TO INSUFFICIENT CORE.', //5X,
     2       'ADDITIONAL CORE REQUIRED =',I10,7H  WORDS)
      CALL MESAGE (-61,0,0)
C
C     CHECK FOR MULTIPLE REFERENCES TO GRID (OR SCALAR) POINT ID NOS.
C     AND SEQUENCE ID NOS. ON SEQGP CARDS
C
 2060 K  = N2
      KK = N2 + FLAG - 1
      JJ = KK - 2
 2080 DO 2285 I = K,JJ,2
      IF (Z(I).LT.0 .OR. I.GE.KK) GO TO 2275
      II = I + 2
      IFAIL = 0
      DO 2270 J = II,KK,2
      IF (Z(I) .NE. Z(J)) GO TO 2270
      IF (IFAIL .NE.   0) GO TO 2260
      IFAIL = 1
      NOGO  = 1
      IF (K .NE. N2) GO TO 2110
      WRITE  (IOUT,2100) UFM,Z(I)
 2100 FORMAT (A23,' 3136, MULTIPLE REFERENCES TO GRID (OR SCALAR) POINT'
     1,      ' ID NO.',I9,'  ON SEQGP CARDS.')
      GO TO 2260
 2110 IDSEQ1 = Z(I)/1000
      IRMNDR = Z(I) - 1000*IDSEQ1
      IF (IRMNDR.NE.0 .AND. MULT.GE.10) GO TO 2140
      IF (AXI .NE. 0) GO TO 2130
      WRITE  (IOUT,2120) UFM,IDSEQ1
 2120 FORMAT (A23,' 3137, MULTIPLE REFERENCES TO SEQUENCE ID NO.',I6,6X,
     1       ' ON SEQGP CARDS.')
      GO TO 2260
 2130 IF (AXI .EQ. 1) WRITE (IOUT,2135) UFM
 2135 FORMAT (A23,' 3137A, SEQGP CARDS WITH MORE THAN ONE SEQID LEVEL ',
     1       'ARE ILLEGAL FOR AXISYSM. OR HYDROELAS. PROBLEM')
      AXI  = 2
      NOGO = 1
      GO TO 2260
 2140 IDSEQ2 = IRMNDR/100
      IRMNDR = IRMNDR - 100*IDSEQ2
      IF (IRMNDR.NE.0 .AND. MULT.GE.100) GO TO 2180
      WRITE  (IOUT,2160) UFM,IDSEQ1,IDSEQ2
 2160 FORMAT (A23,' 3137, MULTIPLE REFERENCES TO SEQUENCE ID NO.',I6,
     1        1H.,I1,5X,'ON SEQGP CARDS.')
      GO TO 2260
 2180 IDSEQ3 = IRMNDR/10
      IRMNDR = IRMNDR - 10*IDSEQ3
      IF (IRMNDR .NE. 0) GO TO 2220
      WRITE  (IOUT,2200) UFM,IDSEQ1,IDSEQ2,IDSEQ3
 2200 FORMAT (A23,' 3137, MULTIPLE REFERENCES TO SEQUENCE ID NO.',I6,
     1       1H.,I1,1H.,I1,4X,'ON SEQGP CARDS.')
      GO TO 2260
 2220 WRITE  (IOUT,2240) UFM,IDSEQ1,IDSEQ2,IDSEQ3,IRMNDR
 2240 FORMAT (A23,' 3137, MULTIPLE REFERENCES TO SEQUENCE ID NO.',I6,
     1        1H.,I1,1H.,I1,1H.,I1,'  ON SEQGP CARDS.')
 2260 Z(J) = -Z(J)
 2270 CONTINUE
C
 2275 IF (JJ.LT.KK .OR. MULT.EQ.1000) GO TO 2285
      L = Z(I)
      IF (MULT   .LE.   10) GO TO 2280
      IF (MOD(L,10) .NE. 0) GO TO 2276
      Z(I) = L/10
      GO TO 2285
 2276 IF (MAXA1 .EQ. 0) GO TO 2285
      MAXA1 = 0
      NOGO  = 1
      WRITE  (IOUT,2277) UFM
 2277 FORMAT (A23,' 2140B, ILLEGAL DATA IN SEQGP CARD, POSSIBLY CAUSED',
     1       ' BY LARGE GRID OR SCALAR POINTS')
      GO TO 2285
 2280 IF (MULT .EQ. 1) GO TO 2282
      IF (MOD(L,100) .NE. 0) GO TO 2276
      Z(I) = L/100
      GO TO 2285
 2282 IF (AXI .EQ. 0) CALL MESAGE (-37,0,NAM)
      IF (MOD(L,1000) .EQ. 0) GO TO 2285
      IF (AXI .EQ. 1) WRITE (IOUT,2135) UFM
      AXI  = 2
      NOGO = 1
 2285 CONTINUE
C
      IF (K .NE. N2) GO TO 2290
      JJ = KK
      K  = K + 1
      GO TO 2080
C
 2290 DO 2300 I = N2,KK,2
      IF (Z(I) .LT. 0) Z(I) = -Z(I)
 2300 CONTINUE
      IF (NOGO .EQ. 1) GO TO 2400
C
C     CHECK TO SEE IF ANY SEQUENCE ID NO. ON SEQGP CARDS IS THE SAME
C     AS A GRID (OR SCALAR) POINT ID NO. THAT HAS NOT BEEN RESEQUENCED
C
      DO 2390 I = K,KK,2
      IF (Z(I) .LT. 0) GO TO 2390
      IDSEQ1 = Z(I)/MULT
      IRMNDR = Z(I) - MULT*IDSEQ1
      IF (IRMNDR .NE. 0) GO TO 2390
      DO 2320 J = N2,KK,2
      IF (IDSEQ1 .EQ. Z(J)) GO TO 2390
 2320 CONTINUE
      DO 2340 J = 1,N1,2
      IF (IDSEQ1 .EQ. Z(J)) GO TO 2360
 2340 CONTINUE
      GO TO 2390
 2360 NOGO = 1
      WRITE  (IOUT,2380) UFM,IDSEQ1
 2380 FORMAT (A23,' 3138, SEQUENCE ID NO.',I6,' ON SEQGP CARDS IS THE ',
     1       'SAME AS A', /5X,'GRID (OR SCALAR) POINT ID NO. THAT HAS ',
     2       'NOT BEEN RESEQUENCED.')
 2390 CONTINUE
 2400 CONTINUE
      I = -1
  220 I = I + 2
      IF (I .GT. FLAG) GO TO 240
      A(1) = Z(N2+I-1)
      A(2) = Z(N2+I  )
      GO TO 1060
  230 Z(2*K) = A(2)
      GO TO 220
C
C     SORT THE CORE TABLE BY INTERNAL GRID PT NO
C     THUS FORMING THE GPL (EXTERNAL GRID PT NOS IN SORT BY INTERNAL NO)
C
  240 IF (NOGO .NE. 0) GO TO 1165
      CALL SORTI (0,0,2,2,Z,N1)
C
C     CLOSE GEOM1. WRITE THE GPL. FIRST RECORD IS A SINGE ENTRIED LIST
C     OF EXTERNAL GRID NOS. IN INTERNAL SORT. SECOND RECORD IS A DOUBLE
C     ENTRIED LIST OF EXTERAL GRID NO., SEQUENCE NO. (SORT IS INTERNAL).
C     ADD THE MULTIPLIER, MULT, TO THE 3RD WORD OF GPL HEADER RECORD
C
  250 IF (NOGMP1 .NE. 0) CALL CLOSE (GEOM1,CLSREW)
  260 CALL FNAME (GPL,A)
      FILE = GPL
      CALL OPEN (*1170,GPL,Z(BUF1),WRTREW)
      A(3) = MULT
      CALL WRITE (GPL,A,3,1)
      DO 270 I = 1,N,2
  270 CALL WRITE (GPL,Z(I),1,0)
      CALL WRITE (GPL,0,0,1)
      CALL WRITE (GPL,Z,N1,1)
      CALL CLOSE (GPL,CLSREW)
      MCB(1) = GPL
      CALL WRTTRL (MCB)
C
C     FORM INTERNAL INDEX FOR EACH EXTERNAL GRID PT. NO.
C
      I = 2
      Z(I) = 1
      IF (N .EQ. 1) GO TO 310
      DO 290 I = 3,N,2
  290 Z(I+1) = Z(I-1) + 1
C
C     TEST TO SEE IF EXTERNAL GRID PT NOS ARE STILL IN EXTERNAL SORT
C     I.E., IF NO SEQGP TABLE, THEN SORT IS MAINTAINED
C     OTHERWISE, SORT ON EXTERNAL GRID NO.
C
      IF (NOSEQ .NE. 0) CALL SORTI (0,0,2,1,Z,N1)
C
C     DETERMINE IF THE GPDT CAN BE HELD IN CORE
C     NWDS= TOTAL NO OF WORDS IN THE GPDT
C     M= MAX NO OF ENTRIES CORE CAN HOLD WITH ONE BUFFER OPEN
C     IF NWDS/7.LE.M,CORE WILL HOLD THE GPDT
C     OTHERWISE THE FILE SORT ROUTINE WILL BE USED
C
  310 NWDS = 7*KN
      M    = (BUF1-N1)/7
      GPFL = 0
      IF (KN .GT. M) GPFL = 7
      FILE = SCR1
C
C     READ THE GRID AND SPOINT TABLES FROM SCR1
C     REPLACE THE EXTERNAL GRID PT NO WITH THE INTERNAL INDEX
C     IF CORE WILL HOLD THE GPDT, USE THE INTERNAL INDEX AS A POINTER
C     OTHERWISE, WRITE THE UNSORTED GPDT ON SCR2
C
      CALL OPEN (*1170,SCR1,Z(BUF1),RDREW)
      FILE = SCR2
      IF (GPFL .NE. 0) CALL OPEN (*1170,SCR2,Z(BUF2),WRTREW)
      FILE = SCR1
      ASSIGN 340 TO NDX
      IERR = 2
      ASSIGN 330 TO NERR
  330 CALL READ (*1180,*370,SCR1,A,7,0,FLAG)
      GO TO 1060
  340 IF (GPFL .NE. 0) GO TO 360
      J = N1 + 7*(A(1)-1)
      DO 350 K = 1,7
      I = J+K
  350 Z(I) = A(K)
      GO TO 330
  360 CALL WRITE (SCR2,A,7,0)
      GO TO 330
  370 IF (NOGO .NE. 0) GO TO 1165
      CALL CLOSE (SCR1,CLSREW)
C
C     OPEN OUTPUT FILE FOR GPDT AND WRITE HEADER DATA
C     IF GPDT IS IN CORE, WRITE IT OUT
C
      FILE = GPDT
      CALL FNAME (GPDT,A)
      CALL OPEN  (*1170,GPDT,Z(BUF1),WRTREW)
      CALL WRITE (GPDT,A,2,1)
      IF (GPFL .NE. 0) GO TO 390
      CALL WRITE (GPDT,Z(IGPDT),NWDS,1)
      GO TO 400
C
C     IF GPDT NOT IN CORE, CALL SORT
C
  390 NFILE(1) = SCR1
      NFILE(2) = CSTM
      NFILE(3) = BGPDT
      CALL CLOSE (SCR2,CLSREW)
      FILE = SCR2
      CALL OPEN  (*1170,SCR2,Z(BUF2),RDREW)
      CALL SORTI (SCR2,GPDT,7,1,Z(IGPDT),BUF2-IGPDT)
      CALL CLOSE (SCR2,CLSREW)
  400 CALL CLOSE (GPDT,CLSREW)
      MCB(1) = GPDT
      CALL WRTTRL (MCB)
C
C     READ THE CORDIJ TABLES INTO CORE (IF PRESENT)
C
      IFL = -1
      M   = ICSDT
      NOLIST = 0
      IF (NOGMP1 .EQ. 0) GO TO 810
      NDX   = BUF1 - 15
      NCORE = BUF1 - 15
      DO 420 I = ICSDT,BUF1
  420 Z(I) = 0
      FILE = GEOMP
      CALL PRELOC (*1170,Z(BUF1),GEOMP)
      DO 440 I = 1,6
      IJ = I + I - 1
      CALL LOCATE (*440,Z(BUF1),CORDIJ(IJ),FLAG)
      IFL = 1
  430 CALL READ (*1180,*440,GEOMP,Z(M),CORD(I),0,FLAG)
      M = M + 16
      IF (M .GT. NCORE) CALL MESAGE (-8,0,GP1AH)
      GO TO 430
  440 CONTINUE
      CALL CLOSE (GEOMP,CLSREW)
      M = M - 16
      NCSDT = M
C
C     TEST FOR PRESENCE OF ANY CORDIJ TABLES
C
      IF (IFL .EQ. -1) GO TO 810
C
C     REPLACE EXTERNAL GRID PT NO IN CORD1J ENTRIES (IF ANY)
C     WITH CORRESPONDING INTERNAL INDEX
C     SAVE A TABLE OF GRID PTS REFERENCED ON CORD1J ENTRIES
C
      JJ   = ICSDT
      ILIST= NCSDT + 16
      II   = ILIST - 1
      NCORE= BUF1  - 3
      IERR = 3
  470 IF (Z(JJ+2).NE.1) GO TO 510
      NOLIST = 1
      ASSIGN 480 TO NDX
      ASSIGN 485 TO NERR
      A(1) = Z(JJ+3)
      SPOINT(2) =  Z(JJ+1)
      GO TO 1060
  480 Z(JJ+3) = A(1)
      Z(II+1) = A(1)
  485 ASSIGN 490 TO NDX
      ASSIGN 495 TO NERR
      A(1) = Z(JJ+4)
      GO TO 1060
  490 Z(JJ+4) = A(1)
      Z(II+2) = A(1)
  495 ASSIGN 500 TO NDX
      ASSIGN 505 TO NERR
      A(1) = Z(JJ+5)
      GO TO 1060
  500 Z(JJ+5) = A(1)
      Z(II+3) = A(1)
  505 II = II+3
      IF (II .GT. NCORE) CALL MESAGE (-8,0,GP1AH)
  510 JJ = JJ + 16
      IF (JJ .LE. NCSDT) GO TO 470
      IF (NOGO .NE.   0) GO TO 1165
C
C     IF ANY CORD1J ENTRIES, PASS THE GPDT AND CREATE A TABLE OF THE
C     REFERENCED GRID PTS. THIS TABLE IS CALLED CSGP
C
      IF (NOLIST .EQ. 0) GO TO 550
      NLIST = II
      ICSGP = NLIST + 1
      CALL SORTI (0,0,1,1,Z(ILIST),ICSGP-ILIST)
      Z(ICSGP) = 0
      JJ = ILIST
      DO 530 KK = ILIST,NLIST
      IF (Z(KK+1) .EQ. Z(KK)) GO TO 530
      Z(JJ) = Z(KK)
      JJ = JJ + 1
  530 CONTINUE
      NLIST = JJ - 1
      ICSGP = JJ
      FILE  = GPDT
      CALL OPEN   (*1170,GPDT,Z(BUF1),RDREW)
      CALL FWDREC (*1180,GPDT)
      NCORE = BUF1 - 5
      I = ILIST
  540 CALL READ (*1180,*1200,GPDT,Z(JJ),7,0,FLAG)
      IF (Z(JJ) .NE. Z(I)) GO TO 540
      JJ = JJ + 5
      IF (JJ .GT. NCORE) CALL MESAGE (-8,0,GP1AH)
      I  = I + 1
      IF (I .LE. NLIST) GO TO 540
      NCSGP = JJ - 5
      CALL CLOSE (GPDT,CLSREW)
C
C     LOOP THRU THE CSDT SOLVING AS MANY COORDINATE SYSTEMS AS POSSIBLE
C     ON EACH PASS.
C
  550 NN = (NCSDT-ICSDT)/16 + 1
      SOLV  = 0
      SOLVP = 0
  560 II = ICSDT
  570 IF (Z(II+2)-2) 580,620,690
C
C     *****  TYPE = 1 *****
C     CHECK TO SEE IF EACH OF THE 3 REFERENCE GRID PTS IS IN BASIC SYS
C     IF SO,CALCULATE THE TRANSFORMATION TO BASIC AND SET COORD SYSTEM
C     AS SOLVED, IF NOT CONTINUE TO NEXT COORDINATE SYSTEM
C
  580 I = 0
  590 K = II + I
      J = ICSGP - 1
  600 IF (Z(J+1) .EQ. Z(K+3)) GO TO 610
      J = J + 5
      IF (J .LT. NCSGP) GO TO 600
      GO TO 1220
  610 IF (Z(J+2).NE.0) GO TO 700
      K = I*3
      AA(K+1) = ZZ(J+3)
      AA(K+2) = ZZ(J+4)
      AA(K+3) = ZZ(J+5)
      I = I+1
      IF (I .LE. 2) GO TO 590
      GO TO 1020
C
C     ***** TYPE = 2 *****
C     CHECK THE DEFINING LOCAL COORDINATE SYSTEM
C     IF BASIC, SOLVE AS IN TYPE=1
C     IF NOT BASIC, FIND THE REFERENCED COORD SYSTEM AND TEST IF THAT
C     SYSTEM IS SOLVED. IF YES, CALCULATE THE TRANSFORMATION TO BASIC
C     IF NO, CONTINUE THRU THE CSDT
C
  620 IF (Z(II+3) .NE. 0) GO TO 640
      DO 630 I = 1,9
      K = II + I
  630 AA(I) = ZZ(K+3)
      GO TO 1020
  640 I = ICSDT
  650 IF (Z(I) .EQ. Z(II+3)) GO TO 660
      I = I + 16
      IF (I .LE. NCSDT) GO TO 650
      GO TO 1230
  660 IF (Z(I+2).NE.3 .OR. Z(I+3).NE.0) GO TO 700
      K = 0
      ASSIGN 680 TO NDX
  670 L = K + II
      AX(1) = ZZ(L+4)
      AX(2) = ZZ(L+5)
      AX(3) = ZZ(L+6)
      IF (Z(I+1)-2) 990,1000,1010
  680 AA(K+1) = AR(1)
      AA(K+2) = AR(2)
      AA(K+3) = AR(3)
      K = K + 3
      IF (K .LE. 6) GO TO 670
      GO TO 1020
C
C     ***** TYPE = 3 *****
C     CHECK THE DEFINING LOCAL COORDINATE SYSTEM
C     IF BASIC, CONTINUE THRU CSDT
C     IF NOT BASIC, ERROR CONDITION
C
  690 IF (Z(II+3) .NE. 0) GO TO 1190
C
C     TEST FOR COMPLETION OF PASS THRU CSDT
C
  700 II = II + 16
      IF (II .LE. NCSDT) GO TO 570
C
C     LOOP THRU THE CSGP (IFPRESENT) AND TRANSFORM ALL
C     POSSIBLE GRID PTS TO BASIC
C
      IF (NOLIST .EQ. 0) GO TO 770
      JJ = ICSGP
  720 IF (Z(JJ+1) .EQ. 0) GO TO 760
      I = ICSDT
  730 IF (Z(JJ+1) .EQ. Z(I)) GO TO 740
      I = I + 16
      IF (I .LE. NCSDT) GO TO 730
      IERR = 6
      SPOINT(1) = Z(JJ  )
      SPOINT(2) = Z(JJ+1)
      GO TO 1190
  740 IF (Z(I+2).NE.3 .OR. Z(I+3).NE.0) GO TO 760
      AX(1) = ZZ(JJ+2)
      AX(2) = ZZ(JJ+3)
      AX(3) = ZZ(JJ+4)
      ASSIGN 750 TO NDX
      IF (Z(I+1)-2) 990,1000,1010
  750 ZZ(JJ+2) = AR(1)
      ZZ(JJ+3) = AR(2)
      ZZ(JJ+4) = AR(3)
      ZZ(JJ+1) = 0
  760 JJ = JJ + 5
      IF (JJ .LE. NCSGP) GO TO 720
C
C     TEST TO SEE IF ALL COORDINATE SYSTEMS SOLVED
C     IF NOT, TEST TO SEE IF ANY NEW SOLUTIONS ON LAST PASS
C     IF NONE, INCONSISTANT DEFINITION OF COORDINATE SYSTEMS
C     OTHERWISE LOOP BACK THRU THE CSDT
C
  770 IF (SOLV .EQ.    NN) GO TO 780
      IF (SOLV .EQ. SOLVP) GO TO 1240
      SOLVP = SOLV
      GO TO 560
C
C     WRITE THE CSTM
C
  780 CALL FNAME (CSTM,A)
      FILE = CSTM
      CALL OPEN  (*1170,CSTM,Z(BUF1),WRTREW)
      CALL WRITE (CSTM,A,2,1)
      DO 800 II = ICSDT,NCSDT,16
      CALL WRITE (CSTM,Z(II),2,0)
      CALL WRITE (CSTM,Z(II+4),12,0)
  800 CONTINUE
      CALL CLOSE (CSTM,CLSREW)
      NOCSTM = NN
      MCB(3) = NN
      MCB(1) = CSTM
      CALL WRTTRL (MCB)
C
C     OPEN EQEXIN AND WRITE HEADER RECORD.
C     THEN WRITE FIRST RECORD (PAIRS OF EXTERNAL GRID NO., INTERNAL NO.
C     IN EXTERNAL SORT).
C
  810 FILE = EQEXIN
      CALL OPEN  (*1170,EQEXIN,Z(BUF1),WRTREW)
      CALL FNAME (EQEXIN,A)
      CALL WRITE (EQEXIN,A,2,1)
      CALL WRITE (EQEXIN,Z,N1,1)
      CALL CLOSE (EQEXIN,CLS)
C
C     A LIST OF DEGREES OF FREEDOM FOR EACH GRID OR SCALAR POINT IS
C     FORMED BEGINNING AT Z(ILIST) BY READING GEOM2 AND USING THE
C     CONNECTION INFORMATION IN CONJUNCTION WITH THE ELEM TABLE IN
C     /GPTA1/.
C
      FILE   = GEOM2
      ILIST0 = ILIST - 1
      NLIST  = ILIST + (NEQEX+1)/2
      IF (NLIST .GE. BUF3) CALL MESAGE (-8,0,GP1AH)
      DO 8102 I = ILIST,NLIST
      Z(I) = 0
 8102 CONTINUE
      JERR = 0
      CALL OPEN   (*8130,GEOM2,Z(BUF1),RDREW)
 8103 CALL FWDREC (*1180,GEOM2)
 8104 CALL ECTLOC (*8130,GEOM2,A,I)
C
C     ELEMENT TYPE LOCATED--PREPARE TO PROCESS EACH ELEMENT
C
      IF (ELEM(I+9) .EQ. 0) GO TO 8103
      J1 = ELEM(I+12)
      NREAD = J1 + ELEM(I+9) - 1
      NSKIP =-(ELEM(I+5 ) - NREAD)
      MAXDOF=  ELEM(I+24)
      ITYPE =  ELEM(I+2 )
C
C     READ CONNECTION DATA FOR ELEMENT AND LOCATE EXT. GRID NBR IN
C     EQEXIN UPDATE DOF LIST FOR EACH GRID NBR
C
 8110 CALL READ (*1180,*8104,GEOM2,A,NREAD,0,M)
      DO 8128 I = J1,NREAD
      IF (A(I) .EQ. 0) GO TO 8128
      CALL BISLOC (*8122,A(I),Z,2,KN,K)
      J = ILIST0 + Z(K+1)
      IF (ITYPE.GE.76 .AND. ITYPE.LE.79) GO TO 8115
C
C     STRUCTURE ELEMENT AND OTHERS
C
      IF (Z(J) .LT. 0) GO TO 8124
      Z(J) = MAX0(Z(J),MAXDOF)
      GO TO 8128
C
C     FLUID ELEMENT (CFHEX1,CFHEX2,CFWEDGE,CFTETRA)
C
 8115 IF (Z(J) .GT. 0) GO TO 8124
C
      Z(J) = -1
      GO TO 8128
 8122 WRITE  (IOUT,8123) UFM,A(1),A(I)
 8123 FORMAT (A23,' 2007, ELEMENT',I8,' REFERENCES UNDEFINED GRID ',
     1       'POINT',I8)
      JERR = JERR + 1
      GO TO 8128
 8124 WRITE  (IOUT,8125) UFM,A(I)
 8125 FORMAT (A23,' 8011, GRID POINT',I8,' HAS BOTH STRUCTURE AND ',
     1       'FLUID ELEMENTS CONNECTED')
      JERR = JERR + 1
 8128 CONTINUE
      CALL READ (*1180,*8104,GEOM2,A,NSKIP,0,M)
      GO TO 8110
C
C     END-OF-FILE ON GEOM2---IF FATAL ERRORS, TERMINATE
C
 8130 CONTINUE
      IF (JERR .NE. 0) CALL MESAGE (-61,A,Z)
C
C     OPEN BGPDT AND SIL. WRITE HEADER RECORDS. OPEN GPDT. SKIP HEADER.
C
      OFFSET = RSHIFT(KN,5)
      CALL FNAME (BGPDT,A)
      CALL FNAME (SIL,A(3))
      FILE = BGPDT
      CALL OPEN (*1170,BGPDT,Z(BUF1),WRTREW)
      FILE = SIL
      CALL OPEN (*1170,SIL,Z(BUF2),WRTREW)
      FILE = GPDT
      CALL OPEN   (*1170,GPDT,Z(BUF3),RDREW)
      CALL FWDREC (*1180,GPDT)
      CALL WRITE  (BGPDT,A,2,1)
      CALL WRITE  (SIL,A(3),2,1)
      LUSET = 1
C
C     READ AN ENTRY FROM THE GPDT.
C     TEST FOR DEFINING COORDINATE SYSTEM.
C
  820 CALL READ (*1180,*970,GPDT,A,7,0,FLAG)
      IF (A(2)) 910,880,830
C
C     COORDINATE SYSTEM NOT BASIC--
C     USE CSDT IN CORE TO TRANSFORM TO BASIC.
C
  830 IF (NOCSTM .EQ. -1) GO TO 850
      I = ICSDT
  840 IF (Z(I).EQ.A(2)) GO TO 860
      I = I + 16
      IF (I .LE. NCSDT) GO TO 840
  850 IERR = 6
      SPOINT(1) = A(1)
      SPOINT(2) = A(2)
      GO TO 1190
  860 AX(1) = AA(3)
      AX(2) = AA(4)
      AX(3) = AA(5)
      ASSIGN 870 TO NDX
      IF (Z(I+1)-2) 990,1000,1010
  870 AA(3) = AR(1)
      AA(4) = AR(2)
      AA(5) = AR(3)
C
C     GRID POINT NOW BASIC--
C     STORE DISPLACEMENT SYSTEM COORD. SYSTEM ID AND SET TYPE.
C     MAKE SURE DISPLACEMENT COORD. SYSTEM IS DEFINED.
C
  880 A(2) = A(6)
      TYPE = 1
      KHR  = ILIST0 + A(1)
      INCR = Z(KHR)
C
C     IF INCR NEGATIVE - SPECIAL HYDROELASTIC GRID POINT WITH SINGLE
C     DEGREE OF FREEDOM
C
      IF (INCR .LT. 0) GO TO 905
      IF (INCR .EQ. 0) INCR = 6
C
C     ///////////////////////////////
C
C     TEMP PATCH
C
      INCR = MAX0(INCR,6)
C
C     ///////////////////////////////
C
      IF (A(2).EQ.0 .AND. ITHERM.EQ.0) GO TO 920
C
C     IF A(2) WHICH EQUALS A(6) IS EQUAL TO -1 THEN A FLUID GRID POINT
C     AS CREATED BY IFP4 IS AT HAND AND HAS ONLY 1 DEGREE OF FREEDOM
C     ..... IF -HEAT- PROBLEM THEN ALL GRIDS HAVE 1 DEGREE OF FREEDOM.
C
      IF (A(2).EQ.(-1) .OR. ITHERM.GT.0) GO TO 905
      IF (NOCSTM .EQ. -1) GO TO 900
      DO 890 IJK = ICSDT,NCSDT,16
      IF (A(2) .EQ. Z(IJK)) GO TO 920
  890 CONTINUE
  900 NOGO = 1
      CALL MESAGE (30,104,A(2))
      GO TO 920
C
C     SCALAR POINT-- SET TYPE.
C
  905 A(2) = 0
      A(6) = 0
  910 TYPE = 2
      INCR = 1
C
C     WRITE ENTRY ON BGPDT AND SIL.
C
  920 CALL WRITE (BGPDT,A(2),4,0)
      CALL WRITE (SIL, LUSET,1,0)
C
C     REPLACE INTERNAL NO. IN EQEXIN WITH CODED SIL NO.
C     THEN INCREMENT SIL NO.
C
      NCODE = 10*LUSET + TYPE
      IF (NOSEQ .NE. 0) GO TO 925
      K = 2*A(1)
      IF (Z(K) - A(1)) 950,960,950
  925 NCODE = -NCODE
      LMT1  = MAX0(2*(A(1)-OFFSET),2)
      DO 930 K = LMT1,N1,2
      IF (Z(K) .EQ. A(1)) GO TO 960
  930 CONTINUE
      DO 940 K = 2,LMT1,2
      IF (Z(K) .EQ. A(1)) GO TO 960
  940 CONTINUE
  950 CALL MESAGE (-30,2,A)
  960 Z(K)  = NCODE
      LUSET = LUSET + INCR
      GO TO 820
C
C     CLOSE BGPDT AND SIL. WRITE TRAILERS.
C
  970 CALL CLOSE (BGPDT,CLSREW)
      CALL CLOSE (SIL,CLSREW)
      CALL CLOSE (GPDT,CLSREW)
      LUSET = LUSET - 1
C    2147483647   = 2**31-1
      IF (LUSET .LE. 2147483647) GO TO 974
      WRITE (IOUT,972) UFM,LUSET
  972 FORMAT (A23,' 3175, TOTAL NUMBER OF DEGREES OF FREEDOM IN THE ',
     1       'PROBLEM (',I11,' ) EXCEEDS 2,147,483,647 (I.E., '
     2       '2**31 - 1)')
  973 FORMAT (A29,' 3175, PROBLEM SIZE,',I8,' DOF''S, EXCEEDS THE OLD ',
     1       'LIMIT OF 65535.', /5X,'GOOD NEWS, JOB WILL CONTINUE')
      CALL MESAGE (-61,0,0)
  974 MCB(1) = BGPDT
      MCB(3) = 0
      CALL WRTTRL (MCB)
      MCB(1) = SIL
      MCB(3) = LUSET
      CALL WRTTRL (MCB)
C
C     IF GRID NOS. ARE RESEQUENCED, SWITCH SIGN ON CODED SIL NO.
C     WRITE SECOND RECORD OF EQEXIN. CLOSE FILE AND WRITE TRAILER.
C
      IF (NOSEQ .EQ. 0) GO TO 978
      DO 976 K = 2,N1,2
      Z(K) = -Z(K)
  976 CONTINUE
  978 FILE = EQEXIN
      CALL OPEN  (*1170,EQEXIN,Z(BUF1),WRT)
      CALL WRITE (EQEXIN,Z,N1,1)
      CALL CLOSE (EQEXIN,CLSREW)
      MCB(1) = EQEXIN
      MCB(3) = 0
      CALL WRTTRL (MCB)
      CALL SSWTCH (36,K)
      IF (K    .EQ. 1) CALL DIAG36 (Z,BUF1,GPL,SIL,EQEXIN)
      IF (NOGO .NE. 0) CALL MESAGE (-61,0,0)
      RETURN
C
C     ABNORMAL EXIT FROM GP1
C
  980 CALL CLOSE (SCR1,CLSREW)
      CALL CLOSE (GEOM1,CLSREW)
      NOCSTM = -1
      RETURN
C
C     ===============================================================
C
C     INTERNAL SUBROUTINE TO TRANSFORM A RECTANGULAR GRID PT TO BASIC
C     I POINTS TO THE CSDT ENTRY WHERE THE TRANSFORMATION IS DEFINED
C     THE GRID PT TO BE TRANSFORMED IS STORED AT AX(1,2,3)
C     THE TRANSFORMED GRID PT WILL BE STORED AT AR(1,2,3)
C
  990 AR(1) = ZZ(I+ 7)*AX(1) + ZZ(I+ 8)*AX(2) + ZZ(I+ 9)*AX(3) + ZZ(I+4)
      AR(2) = ZZ(I+10)*AX(1) + ZZ(I+11)*AX(2) + ZZ(I+12)*AX(3) + ZZ(I+5)
      AR(3) = ZZ(I+13)*AX(1) + ZZ(I+14)*AX(2) + ZZ(I+15)*AX(3) + ZZ(I+6)
      GO TO NDX, (680,750,870)
C
C     INTERNAL SUBROUTINE TO TRANSFORM A CYLINDRICAL GRID PT TO BASIC
C     R,THETA,Z IS STORED AX(1,2,3)
C
 1000 R     = AX(1)
      AX(2) = DEGRA*AX(2)
      AX(1) = R*COS(AX(2))
      AX(2) = R*SIN(AX(2))
      GO TO 990
C
C
C     INTERNAL SUBROUTINE TO TRANSFORM A SPHERICAL GRID PT TO BASIC
C     RHO,THETA,PHI IS STORED AT AX(1,2,3)
C
 1010 AX(2) = DEGRA*AX(2)
      AX(3) = DEGRA*AX(3)
      RSTH  = AX(1)*SIN(AX(2))
      RCTH  = AX(1)*COS(AX(2))
      AX(1) = RSTH *COS(AX(3))
      AX(2) = RSTH *SIN(AX(3))
      AX(3) = RCTH
      GO TO 990
C
C
C     INTERNAL SUBROUTINE TO CALCULATE THE 3X3 TRANSFORMATION MATRIX
C     AND 3X1 TRANSLATION VECTOR GIVEN THREE POINTS IN THE BASIC SYSTEM
C     THE RESULTS ARE STORED BACK IN THE CSDT
C
C     STORE R0 = A IN THE CSDT
C
 1020 ZZ(II+4) = AA(1)
      ZZ(II+5) = AA(2)
      ZZ(II+6) = AA(3)
C
C     FORM B - A
C
      DO 1030 I = 1,3
 1030 AK(I) = AB(I) - AA(I)
C
C     FORM K = (B - A)/LENGTH(B - A)
C     FORM C - A
C
      LENGTH = SQRT(AK(1)**2 + AK(2)**2 + AK(3)**2)
      DO 1040 I = 1,3
      AK(I) = AK(I)/LENGTH
 1040 AC(I) = AC(I) - AA(I)
C
C     FORM K X (C - A)
C
      AJ(1) = AK(2)*AC(3) - AK(3)*AC(2)
      AJ(2) = AK(3)*AC(1) - AK(1)*AC(3)
      AJ(3) = AK(1)*AC(2) - AK(2)*AC(1)
C
C     FORM J = (K X (C-A))/LENGTH(K X (C-A))
C
      LENGTH =  SQRT(AJ(1)**2 + AJ(2)**2 + AJ(3)**2)
      DO 1050 I = 1,3
 1050 AJ(I) = AJ(I)/LENGTH
C
C     FORM I = J X K
C
      AI(1) = AJ(2)*AK(3) - AJ(3)*AK(2)
      AI(2) = AJ(3)*AK(1) - AJ(1)*AK(3)
      AI(3) = AJ(1)*AK(2) - AJ(2)*AK(1)
C
C     STORE 3X3 ROTATION MATRIX = ((IX,JX,KX),(IY,JY,KY),(IZ,JZ,KZ))
C     IN THE CSDT
C
      ZZ(II+ 7) = AI(1)
      ZZ(II+ 8) = AJ(1)
      ZZ(II+ 9) = AK(1)
      ZZ(II+10) = AI(2)
      ZZ(II+11) = AJ(2)
      ZZ(II+12) = AK(2)
      ZZ(II+13) = AI(3)
      ZZ(II+14) = AJ(3)
      ZZ(II+15) = AK(3)
C
C     SET WD 3 OF CSDT = 3 AND WD 4 = 0 TO INDICATE  SOLVED SYSTEM
C     INCREMENT SOLVED SYSTEM COUNT
C
      Z(II+2) = 3
      Z(II+3) = 0
      SOLV = SOLV + 1
      GO TO 700
C
C
C     INTERNAL SUBROUTINE TO PERFORM BINARY SEARCH ON FIRST ENTRY
C     OF A DOUBLE ENTRIED TABLE STORED AT Z(1) THRU Z(N+1)
C
 1060 KLO = 1
      KHI = KN
 1070 K = (KLO+KHI+1)/2
 1080 IF (A(1)-Z(2*K-1)) 1090,1150,1100
 1090 KHI = K
      GO TO 1110
 1100 KLO = K
 1110 IF (KHI-KLO-1) 1160,1120,1070
 1120 IF (K.EQ.KLO) GO TO 1130
      K = KLO
      GO TO 1140
 1130 K = KHI
 1140 KLO = KHI
      GO TO 1080
 1150 A(1) = Z(2*K)
      GO TO NDX,  (230,340,480,490,500)
 1160 CALL MESAGE (30,IERR,A(1))
      NOGO = 1
      GO TO NERR, (220,330,485,495,505)
 1165 CALL MESAGE (-61,0,0)
C
C
C     FATAL ERROR MESAGES
C
 1170 NDX = -1
      GO TO 1210
 1180 NDX = -2
      GO TO 1210
 1190 CALL MESAGE (-30,IERR,SPOINT)
 1200 NDX = -3
      GO TO 1210
 1210 CALL MESAGE (NDX,FILE,GP1AH)
 1220 SPOINT(1) = Z(K+3)
      SPOINT(2) = Z(II )
      IERR = 3
      GO TO 1190
 1230 SPOINT(1) = Z(II  )
      SPOINT(2) = Z(II+3)
      IERR = 4
      GO TO 1190
 1240 SPOINT(1) = 0
      SPOINT(2) = 0
      IERR = 5
      GO TO 1190
 1250 SPOINT(1) = A(1)
      SPOINT(2) = 0
      IERR = 12
      GO TO 1190
      END


text:
!/ ------------------------------------------------------------------- /
      MODULE W3PRO3MD
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         29-May-2009 |
!/                  +-----------------------------------+
!/
!/    27-Feb-2000 : Origination.                        ( version 2.08 )
!/    17-Sep-2000 : Clean-up.                           ( version 2.13 )
!/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
!/    16-Oct-2002 : Change INTENT for ATRN in W3XYP3.   ( version 3.00 )
!/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
!/    01-Aug-2003 : Moving grid GSE correction.         ( version 3.03 )
!/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/    07-Sep-2005 : Upgrade XY boundary conditions.     ( version 3.08 )
!/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
!/    05-Mar-2008 : Added NEC sxf90 compiler directives.
!/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
!/    01-Apr-2008 : Bug fix W3MAP3 MAPSTA range check.  ( version 3.13 )
!/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
!/
!/    Copyright 2009 National Weather Service (NWS),
!/       National Oceanic and Atmospheric Administration.  All rights
!/       reserved.  WAVEWATCH III is a trademark of the NWS.
!/       No unauthorized use without permission.
!/
!  1. Purpose :
!
!     Bundles routines for third order porpagation scheme in single
!     module.
!
!  2. Variables and types :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      TRNMIN    R.P.  Private   Minimum transparancy for local
!                                switching off of averaging.
!     ----------------------------------------------------------------
!
!     Also work arrays for W3KTP3 (private).
!
!  3. Subroutines and functions :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      W3MAP3    Subr. Public   Set up auxiliary maps.
!      W3MAPT    Subr. Public   Set up transparency map for GSE.
!      W3XYP3    Subr. Public   Third order spatial propagation.
!      W3KTP3    Subr. Public   Third order spectral propagation.
!     ----------------------------------------------------------------
!
!  4. Subroutines and functions used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      STRACE    Subr. W3SERVMD Subroutine tracing.
!      W3QCK1    Subr. W3UQCKMD Regular grid UQ scheme.
!      W3QCK2    Subr. W3UQCKMD Irregular grid UQ scheme.
!      W3QCK3    Subr. W3UQCKMD Regular grid UQ scheme + obstructions.
!     ----------------------------------------------------------------
!
!  5. Remarks :
!
!     - The averaging is not performed around semi-transparent grid
!       points to avoid that leaking through barriers occurs.
!
!  6. Switches :
!
!       !/C90   Cray FORTRAN 90 compiler directives.
!       !/NEC   NEC SXF90 compiler directives.
!
!       !/MGP   Moving grid corrections.
!       !/MGG   Moving grid corrections.
!
!       !/S     Enable subroutine tracing.
!       !/Tn    Enable test output.
!
!  7. Source code :
!
!/ ------------------------------------------------------------------- /
!/
!/ Public variables
!/
      PUBLIC
!/
!/ Private data
!/
      REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3MAP3
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         01-Apr-2008 |
!/                  +-----------------------------------+
!/
!/    27-Feb-2000 : Origination.                        ( version 2.08 )
!/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
!/                  (array allocation only.)
!/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
!/    01-Apr-2008 : Bug fix sec. 4 MAPSTA range check.  ( version 3.13 )
!/
!  1. Purpose :
!
!     Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
!
!  2. Method :
!
!     MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
!     grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
!     devided in three parts. For MAPX2, these ranges are :
!
!         1    - NMX0  Counters IXY for which grid point (IX,IY) and
!                      (IX+1,IY) both are active grid points.
!       NMX0+1 - NMX1  Id. only (IX,IY) active.
!       NMX1+1 - NMX2  Id. only (IX+1,IY) active.
!
!     The array MAPY2 has a similar layout varying IY instead of IX.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      W3WAVE    Subr. W3WAVEMD Wave model routine.
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!  7. Remarks :
!
!  8. Structure :
!
!     ------------------------------------------------------
!      1.   Map MAPX2
!        a  Range 1 to NMX0
!        b  Range NMX0+1 to NMX1
!        c  Range NMX1+1 to NMX2
!      2.   Map MAPY2
!        a  Range 1 to NMY0
!        b  Range NMY0+1 to NMY1
!        c  Range NMY1+1 to NMY2
!      3.   Map MAPAXY
!      4.   Map MAPCXY
!      5.   Maps for intra-spectral propagation
!        a  MAPTH2, MAPATK
!        b  MAPWN2
!     ------------------------------------------------------
!
!  9. Switches :
!
!     !/S   Enable subroutine tracing.
!     !/T   Enable test output.
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
      USE W3GDATMD, ONLY: NK, NTH, NSPEC, NX, NY, NSEA, MAPSTA, MAPSF
      USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT,   &
                          NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY,        &
                          MAPTH2, MAPWN2
      USE W3ODATMD, ONLY: NDST
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IX, IY, IXY0, IX2, IY2, IX0, IY0,    &
                                 ISEA, IK, ITH, ISP, ISP0, ISP2, NCENTC
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Map MAPX2 ------------------------------------------------------ *
! 1.a Range 1 to NMX0
!
      NMX0   = 0
      DO IX=1, NX
        IXY0   = (IX-1)*NY
        IX2    = 1 + MOD(IX,NX)
        DO IY=2, NY-1
          IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY,IX2).EQ.1 ) THEN
              NMX0   = NMX0 + 1
              MAPX2(NMX0) = IXY0 + IY
            END IF
          END DO
        END DO
!
! 1.b Range NMX0+1 to NMX1
!
      NMX1   = NMX0
      DO IX=1, NX
        IXY0   = (IX-1)*NY
        IX2    = 1 + MOD(IX,NX)
        DO IY=2, NY-1
          IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY,IX2).NE.1 ) THEN
              NMX1   = NMX1 + 1
              MAPX2(NMX1) = IXY0 + IY
            END IF
          END DO
        END DO
!
! 1.c Range NMX1+1 to NMX2
!
      NMX2   = NMX1
      DO IX=1, NX
        IXY0   = (IX-1)*NY
        IX2    = 1 + MOD(IX,NX)
        DO IY=2, NY-1
          IF ( MAPSTA(IY,IX).NE.1 .AND. MAPSTA(IY,IX2).EQ.1 ) THEN
              NMX2   = NMX2 + 1
              MAPX2(NMX2) = IXY0 + IY
            END IF
          END DO
        END DO
!
! 2.  Map MAPY2 ------------------------------------------------------ *
! 2.a Range 1 to NMY0
!
      NMY0   = 0
      DO IX=1, NX
        IXY0   = (IX-1)*NY
        DO IY=1, NY-1
          IY2    = IY + 1
          IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY2,IX).EQ.1 ) THEN
              NMY0   = NMY0 + 1
              MAPY2(NMY0) = IXY0 + IY
            END IF
          END DO
        END DO
!
! 2.b Range NMY0+1 to NMY1
!
      NMY1   = NMY0
      DO IX=1, NX
        IXY0   = (IX-1)*NY
        DO IY=1, NY-1
          IY2    = IY + 1
          IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY2,IX).NE.1 ) THEN
              NMY1   = NMY1 + 1
              MAPY2(NMY1) = IXY0 + IY
            END IF
          END DO
        END DO
!
! 2.c Range NMY1+1 to NMY2
!
      NMY2   = NMY1
      DO IX=1, NX
        IXY0   = (IX-1)*NY
        DO IY=1, NY-1
          IY2    = IY + 1
          IF ( MAPSTA(IY,IX).NE.1 .AND. MAPSTA(IY2,IX).EQ.1 ) THEN
              NMY2   = NMY2 + 1
              MAPY2(NMY2) = IXY0 + IY
            END IF
           END DO
         END DO
!
! 3.  Map MAPAXY ----------------------------------------------------- *
!
      NACT   = 0
      DO IX=1, NX
        IY0    = (IX-1)*NY
        DO IY=2, NY-1
          IF ( MAPSTA(IY,IX).EQ.1 ) THEN
              NACT         = NACT + 1
              MAPAXY(NACT) = IY0 + IY
            END IF
          END DO
        END DO
!
! 4.  Map MAPCXY ----------------------------------------------------- *
!
      NCENT  = 0
      NCENTC = NSEA
      MAPCXY = 0
!
      DO ISEA=1,  NSEA
        IX      = MAPSF(ISEA,1)
        IX0    = IX-1
        IX2    = IX+1
        IY      = MAPSF(ISEA,2)
        IY0    = IY-1
        IY2    = IY+1
        IF ( IX .EQ. NX ) IX2 = 1
        IF ( IX .EQ. 1 ) IX0 = NX
        IF ( MAPSTA(IY,IX).EQ.2 .OR. MAPSTA(IY,IX).LT.0 ) THEN
            MAPCXY(NCENTC) = ISEA
            NCENTC = NCENTC - 1
          ELSE IF ( MAPSTA(IY0,IX0).GE.1 .AND.                     &
                    MAPSTA(IY0,IX ).GE.1 .AND.                     &
                    MAPSTA(IY0,IX2).GE.1 .AND.                     &
                    MAPSTA(IY ,IX0).GE.1 .AND.                     &
                    MAPSTA(IY ,IX2).GE.1 .AND.                     &
                    MAPSTA(IY2,IX0).GE.1 .AND.                     &
                    MAPSTA(IY2,IX ).GE.1 .AND.                     &
                    MAPSTA(IY2,IX2).GE.1 ) THEN
            NCENT  = NCENT + 1
            MAPCXY(NCENT) = ISEA
          ELSE
            MAPCXY(NCENTC) = ISEA
            NCENTC = NCENTC - 1
          END IF
        END DO
!
! 5.  Maps for intra-spectral propagation ---------------------------- *
!
      IF ( MAPTH2(1) .NE. 0 ) RETURN
!
! 5.a MAPTH2 and MAPBTK
!
      DO IK=1, NK
        DO ITH=1, NTH
          ISP    = ITH + (IK-1)*NTH
          ISP2   = (IK+1) + (ITH-1)*(NK+2)
          MAPTH2(ISP) = ISP2
          END DO
        END DO
!
! 5.b MAPWN2
!
      ISP0   = 0
      DO IK=1, NK-1
        DO ITH=1, NTH
          ISP0   = ISP0 + 1
          ISP2   = (IK+1) + (ITH-1)*(NK+2)
          MAPWN2(ISP0) = ISP2
          END DO
        END DO
!
      DO ITH=1, NTH
        ISP0   = ISP0 + 1
        ISP2   = NK+1 + (ITH-1)*(NK+2)
        MAPWN2(ISP0) = ISP2
        END DO
!
      DO ITH=1, NTH
        ISP0   = ISP0 + 1
        ISP2   = 1 + (ITH-1)*(NK+2)
        MAPWN2(ISP0) = ISP2
        END DO
!
      RETURN
!
! Formats
!
!/
!/ End of W3MAP3 ----------------------------------------------------- /
!/
      END SUBROUTINE W3MAP3
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3MAPT
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         17-Dec-2004 |
!/                  +-----------------------------------+
!/
!/    10-Dec-2001 : Origination.                        ( version 2.14 )
!/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/
!  1. Purpose :
!
!     Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine
!     GSE alleviation with obstructions.
!
!  2. Method :
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      W3WAVE    Subr. W3WAVEMD Wave model routine.
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!     !/S   Enable subroutine tracing.
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
      USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF
      USE W3ADATMD, ONLY: ATRNX, ATRNY, MAPTRN
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: ISEA, IXY
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Map MAPTRN ----------------------------------------------------- *
!
      DO ISEA=1, NSEA
        IXY         = MAPSF(ISEA,3)
        MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNY(IXY,-1) ,              &
                           ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
        END DO
!
      RETURN
!
! Formats
!/
!/ End of W3MAPT ----------------------------------------------------- /
!/
      END SUBROUTINE W3MAPT
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3XYP3 ( ISP, FACX, FACY, DTG, MAPSTA, MAPFS, VQ,    &
                          VGX, VGY )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         05-Mar-2008 |
!/                  +-----------------------------------+
!/
!/    27-Feb-2000 : Origination.                        ( version 2.08 )
!/    17-Sep-2000 : Clean-up.                           ( version 2.13 )
!/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
!/    16-Oct-2002 : Change INTENT for ATRNRX/Y.         ( version 3.00 )
!/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
!/    01-Aug-2003 : Moving grid GSE correction.         ( version 3.03 )
!/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/    07-Sep-2005 : Upgrade XY boundary conditions.     ( version 3.08 )
!/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
!/    05-Mar-2008 : Added NEC sxf90 compiler directives.
!/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
!/
!  1. Purpose :
!
!     Propagation in phyiscal space for a given spectral component.
!
!  2. Method :
!
!     Third-order ULTIMATE QUICKEST scheme with averaging.
!
!  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.
!       MAPSTA  I.A.   I   Grid point status map.
!       MAPFS   I.A.   I   Storage map.
!       VQ      R.A.  I/O  Field to propagate.
!       VGX/Y   Real   I   Speed of grid.
!     ----------------------------------------------------------------
!
!     Local variables.
!     ----------------------------------------------------------------
!       NTLOC   Int   Number of local time steps.
!       DTLOC   Real  Local propagation time step.
!       VCFL0X  R.A.  Local courant numbers for absolute group vel.
!                     using local X-grid step.
!       VCFL0Y  R.A.  Id. in Y.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!       W3QCK3   Actual propagation algorithm
!
!       STRACE   Service routine.
!
!  5. Called by :
!
!       W3WAVE   Wave model routine.
!
!  6. Error messages :
!
!       None.
!
!  7. Remarks :
!
!     - Note that the ULTIMATE limiter does not guarantee non-zero
!       energies.
!     - The present scheme shows a strong distorsion when propaga-
!       ting a field under an angle with the grid in a truly 2-D
!       fashion. Propagation is therefore split along the two
!       axes.
!     - Two boundary treatments are available. The first uses real
!       boundaries in each space. In this case waves will not
!       penetrate in narrow straights under an angle with the grid.
!       This behavior is improved by using a 'soft' option, in
!       which the 'X' or 'Y' sweep allows for energy to go onto
!       the land. This improves the above behavior, but implies
!       that X-Y connenctions are required in barriers for them
!       to become inpenetrable.
!     - Note that unlike in W3XYP2, isotropic diffusion is never
!       used for growth.
!
!  8. Structure :
!
!     ---------------------------------------------
!       1.  Preparations
!         a Set constants
!         b Initialize arrays
!       2.  Prepare arrays
!         a Velocities and 'Q'
!       3.  Loop over sub-steps
!       ----------------------------------------
!         a Average
!         b Propagate
!         c Update boundaries.
!       ----------------------------------------
!       4.  Store Q field in spectra
!     ---------------------------------------------
!
!  9. Switches :
!
!       !/NEC   Enable NEC SXF90 compiler directives.
!
!       !/S     Enable subroutine tracing.
!
!       !/MGP   Moving grid corrections.
!       !/MGG   Moving grid corrections.
!
!       !/T     Enable general test output.
!       !/T0    Dump of precalcaulted interpolation data.
!       !/T1    Dump of input field and fluxes.
!       !/T2    Dump of output field.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3CONSTANTS
!
      USE W3TIMEMD, ONLY: DSEC21
!
      USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF, DTCFL, CLAT, CLATS,    &
                          GLOBAL, FLCX, FLCY, NK, NTH, DTH, XFR,      &
                          ECOS, ESIN, SIG, WDCG, WDTH, PFMOVE
      USE W3WDATMD, ONLY: TIME
      USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT,   &
                          NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY,        &
                          MAPTRN, CG, CX, CY, ATRNX, ATRNY, ITIME
      USE W3IDATMD, ONLY: FLCUR
      USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPI0, TBPIN,       &
                          ISBPI, BBPI0, BBPIN
      USE W3UQCKMD
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
      REAL, INTENT(IN)        :: FACX, FACY, DTG, VGX, VGY
      REAL, INTENT(INOUT)     :: VQ(1-NY:NY*(NX+2))
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: ITH, IK, NTLOC, ITLOC, ISEA, IXY, IP
      INTEGER                 :: IX, IY, IXC, IYC, IIXY1(NY),         &
                                 IIXY2(NY), IIXY3(NY), IIXY4(NY), IBI
      REAL                    :: CG0, CGA, CGN, CGX, CGY, CXC, CYC,   &
                                 CXMIN, CXMAX, CYMIN, CYMAX
      REAL                    :: CGC, FGSE = 1., FTH, FCG
      REAL                    :: DTLOC, CCOS, CSIN, CCURX, CCURY
      REAL                    :: DXCGN, DYCGN, DXCGS, DYCGS, DXCGC,   &
                                 DYCGC, RDI1(NY), RDI2(NY), RDI3(NY), &
                                 RDI4(NY), RD1, RD2, RD3, RD4
      LOGICAL                 :: YFIRST
!/
!/ Automatic work arrays
!/
      INTEGER                 :: MAPSTX(1-2*NY:NY*(NX+2))
      REAL                    :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
                                 AQ(1-NY:NY*(NX+2))
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Preparations --------------------------------------------------- *
! 1.a Set constants
!
      ITH    = 1 + MOD(ISP-1,NTH)
      IK     = 1 + (ISP-1)/NTH
!
      CG0    = 0.575 * GRAV / SIG(1)
      CGA    = 0.575 * GRAV / SIG(IK)
      CGX    = CGA * ECOS(ITH)
      CGY    = CGA * ESIN(ITH)
      CGC    = SQRT ( CGX**2 + CGY**2 )
!
      IF ( FLCUR ) THEN
          CXMIN  = MINVAL ( CX(1:NSEA) )
          CXMAX  = MAXVAL ( CX(1:NSEA) )
          CYMIN  = MINVAL ( CY(1:NSEA) )
          CYMAX  = MAXVAL ( CY(1:NSEA) )
          IF ( ABS(CGX+CXMIN) .GT. ABS(CGX+CXMAX) ) THEN
              CGX    = CGX + CXMIN
            ELSE
              CGX    = CGX + CXMAX
            END IF
          IF ( ABS(CGY+CYMIN) .GT. ABS(CGY+CYMAX) ) THEN
              CGY    = CGY + CYMIN
            ELSE
              CGY    = CGY + CYMAX
            END IF
          CXC    = MAX ( ABS(CXMIN) , ABS(CXMAX) )
          CYC    = MAX ( ABS(CYMIN) , ABS(CYMAX) )
        ELSE
          CXC    = 0.
          CYC    = 0.
        END IF
!
      CGN    = MAX ( ABS(CGX) , ABS(CGY) , CXC, CYC, 0.001*CG0 )
!
      NTLOC  = 1 + INT(DTG/(DTCFL*CG0/CGN))
      DTLOC  = DTG / REAL(NTLOC)
!
      CCOS   = FACX * ECOS(ITH) / REAL(NTLOC)
      CSIN   = FACY * ESIN(ITH) / REAL(NTLOC)
!
      CCURX  = FACX / REAL(NTLOC)
      CCURY  = FACY / REAL(NTLOC)
!
      YFIRST = MOD(ITIME,2) .EQ. 0
!
! 1.b Initialize arrays
!
      VLCFLX = 0.
      VLCFLY = 0.
!
      MAPSTX(1:NX*NY) = MAPSTA(1:NX*NY)
!
      IF ( GLOBAL ) THEN
          MAPSTX(1-2*NY:0)            = MAPSTA((NX-2)*NY+1:NX*NY)
          MAPSTX(NX*NY+1:NX*NY+2*NY)  = MAPSTA(1:2*NY)
        ELSE
          MAPSTX(1-2*NY:0)            = 0
          MAPSTX(NX*NY+1:NX*NY+2*NY)  = 0
        END IF
!
! 1.c Pre-calculate interpolation info
!
      FTH    =   FGSE * WDTH * DTH / REAL(NTLOC)
      FCG    =   FGSE * WDCG * 0.5 * (XFR-1./XFR) / REAL(NTLOC)
!
      DO IY=1, NY
!
! 1.c.1 Normal and parallel width ...
!
        DXCGN  = - FTH * FACX * ESIN(ITH) / CLAT(IY)
        DYCGN  =   FTH * FACY * ECOS(ITH)
        DXCGS  =   FCG * CCOS / CLAT(IY)
        DYCGS  =   FCG * CSIN
!
! 1.c.2 "Sum" corner (and mirror image) ...
!
        DXCGC  = DXCGN + DXCGS
        DYCGC  = DYCGN + DYCGS
!
        IXC    = NY
        IF ( DXCGC .LT. 0. ) IXC = - IXC
        IYC    = 1
        IF ( DYCGC .LT. 0. ) IYC = - IYC
!
        IIXY1(IY) = IXC + IYC
        IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
            IIXY2(IY) = IXC
            RDI1 (IY) = ABS(DYCGC/DXCGC)
            RDI2 (IY) = ABS(DXCGC)
          ELSE
            IIXY2(IY) = IYC
            IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
                RDI1(IY) = ABS(DXCGC/DYCGC)
              ELSE
                RDI1(IY) = 1.
              END IF
            RDI2(IY) = ABS(DYCGC)
          END IF
!
! 1.c.2 "Difference" corner (and mirror image) ...
!
        DXCGC  = DXCGN - DXCGS
        DYCGC  = DYCGN - DYCGS
!
        IXC    = NY
        IF ( DXCGC .LT. 0. ) IXC = - IXC
        IYC    = 1
        IF ( DYCGC .LT. 0. ) IYC = - IYC
!
        IIXY3(IY) = IXC + IYC
        IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
            IIXY4(IY) = IXC
            RDI3 (IY) = ABS(DYCGC/DXCGC)
            RDI4 (IY) = ABS(DXCGC)
          ELSE
            IIXY4(IY) = IYC
            IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
                RDI3(IY) = ABS(DXCGC/DYCGC)
              ELSE
                RDI3(IY) = 1.
              END IF
            RDI4(IY) = ABS(DYCGC)
          END IF
!
        END DO
!
! 2.  Calculate velocities and diffusion coefficients ---------------- *
! 2.a Velocities
!
!     Q     = ( A / CG * CLATS )
!     LCFLX = ( COS*CG / CLATS ) * DT / DX
!     LCFLY = (     SIN*CG )     * DT / DY
!
      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)
        END DO
!
      IF ( FLCUR ) THEN
          DO ISEA=1, NSEA
            IXY         = MAPSF(ISEA,3)
            VLCFLX(IXY) = VLCFLX(IXY) + CCURX*CX(ISEA)/CLATS(ISEA)
            VLCFLY(IXY) = VLCFLY(IXY) + CCURY*CY(ISEA)
           END DO
        END IF
!
! 3.  Loop over sub-steps -------------------------------------------- *
!
      DO ITLOC=1, NTLOC
!
! 3.a Average
!
        AQ     = VQ
        VQ     = 0.
!
! 3.a.1 Central points
!
        DO IP=1, NCENT
          ISEA    = MAPCXY(IP)
          IXY     = MAPSF(ISEA,3)
          IY      = MAPSF(ISEA,2)
          IF ( MAPTRN(IXY) ) THEN
              VQ(IXY) = AQ(IXY)
            ELSE
              RD1     = RDI1(IY)
              RD2     = MIN ( 1. , RDI2(IY) * CG(IK,ISEA) )
              RD3     = RDI3(IY)
              RD4     = MIN ( 1. , RDI4(IY) * CG(IK,ISEA) )
              VQ(IXY          ) = VQ(IXY          )                   &
                                   + AQ(IXY) * (3.-RD2-RD4)/3.
              VQ(IXY+IIXY1(IY)) = VQ(IXY+IIXY1(IY))                   &
                                   + AQ(IXY) * RD2*RD1/6.
              VQ(IXY+IIXY2(IY)) = VQ(IXY+IIXY2(IY))                   &
                                   + AQ(IXY) * (1.-RD1)*RD2/6.
              VQ(IXY+IIXY3(IY)) = VQ(IXY+IIXY3(IY))                   &
                                   + AQ(IXY) * RD4*RD3/6.
              VQ(IXY+IIXY4(IY)) = VQ(IXY+IIXY4(IY))                   &
                                   + AQ(IXY) * (1.-RD3)*RD4/6.
              VQ(IXY-IIXY1(IY)) = VQ(IXY-IIXY1(IY))                   &
                                   + AQ(IXY) * RD2*RD1/6.
              VQ(IXY-IIXY2(IY)) = VQ(IXY-IIXY2(IY))                   &
                                   + AQ(IXY) * (1.-RD1)*RD2/6.
              VQ(IXY-IIXY3(IY)) = VQ(IXY-IIXY3(IY))                   &
                                   + AQ(IXY) * RD4*RD3/6.
              VQ(IXY-IIXY4(IY)) = VQ(IXY-IIXY4(IY))                   &
                                   + AQ(IXY) * (1.-RD3)*RD4/6.
            END IF
          END DO
!
! 3.a.2 Near-coast points
!
        DO IP=NCENT+1, NSEA
          ISEA    = MAPCXY(IP)
          IX      = MAPSF(ISEA,1)
          IY      = MAPSF(ISEA,2)
          IXY     = MAPSF(ISEA,3)
          IF ( MAPSTA(IXY) .LE. 0 ) CYCLE
          IF ( MAPTRN(IXY) ) THEN
              VQ(IXY) = AQ(IXY)
            ELSE
              RD1     = RDI1(IY)
              RD3     = RDI3(IY)
              RD2     = MIN ( 1. , RDI2(IY) * CG(IK,ISEA) )
              RD4     = MIN ( 1. , RDI4(IY) * CG(IK,ISEA) )
              VQ(IXY          ) = VQ(IXY          )                   &
                                   + AQ(IXY) * (3.-RD2-RD4)/3.
!
              IXC    = SIGN(NY,IIXY1(IY))
              IYC    = IIXY1(IY) - IXC
              IF ( MAPSTX(IXY+IIXY1(IY)) .GE. 1 .AND.                 &
                   .NOT. ( MAPSTX(IXY+IXC).LE.0 .AND.                 &
                           MAPSTX(IXY+IYC).LE.0 ) ) THEN
                       VQ(IXY+IIXY1(IY)) = VQ(IXY+IIXY1(IY))          &
                                            + AQ(IXY) * RD2*RD1/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * RD2*RD1/6.
                END IF
              IF ( MAPSTX(IXY-IIXY1(IY)) .GE. 1 .AND.                 &
                   .NOT. ( MAPSTX(IXY-IXC).LE.0 .AND.                 &
                           MAPSTX(IXY-IYC).LE.0 ) ) THEN
                       VQ(IXY-IIXY1(IY)) = VQ(IXY-IIXY1(IY))          &
                                            + AQ(IXY) * RD2*RD1/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * RD2*RD1/6.
                END IF

              IF ( MAPSTX(IXY+IIXY2(IY)) .GE. 1 ) THEN
                       VQ(IXY+IIXY2(IY)) = VQ(IXY+IIXY2(IY))          &
                                            + AQ(IXY) * (1.-RD1)*RD2/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * (1.-RD1)*RD2/6.
                END IF
              IF ( MAPSTX(IXY-IIXY2(IY)) .GE. 1 ) THEN
                       VQ(IXY-IIXY2(IY)) = VQ(IXY-IIXY2(IY))          &
                                            + AQ(IXY) * (1.-RD1)*RD2/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * (1.-RD1)*RD2/6.
                END IF
!
              IXC    = SIGN(NY,IIXY3(IY))
              IYC    = IIXY3(IY) - IXC
              IF ( MAPSTX(IXY+IIXY3(IY)) .GE. 1 .AND.                 &
                   .NOT. ( MAPSTX(IXY+IXC).LE.0 .AND.                 &
                           MAPSTX(IXY+IYC).LE.0 ) ) THEN
                       VQ(IXY+IIXY3(IY)) = VQ(IXY+IIXY3(IY))          &
                                            + AQ(IXY) * RD4*RD3/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * RD4*RD3/6.
                END IF
              IF ( MAPSTX(IXY-IIXY3(IY)) .GE. 1 .AND.                 &
                   .NOT. ( MAPSTX(IXY-IXC).LE.0 .AND.                 &
                           MAPSTX(IXY-IYC).LE.0 ) ) THEN
                       VQ(IXY-IIXY3(IY)) = VQ(IXY-IIXY3(IY))          &
                                            + AQ(IXY) * RD4*RD3/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * RD4*RD3/6.
                END IF
!
              IF ( MAPSTX(IXY+IIXY4(IY)) .GE. 1 ) THEN
                       VQ(IXY+IIXY4(IY)) = VQ(IXY+IIXY4(IY))          &
                                            + AQ(IXY) * (1.-RD3)*RD4/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * (1.-RD3)*RD4/6.
                END IF
              IF ( MAPSTX(IXY-IIXY4(IY)) .GE. 1 ) THEN
                       VQ(IXY-IIXY4(IY)) = VQ(IXY-IIXY4(IY))          &
                                            + AQ(IXY) * (1.-RD3)*RD4/6.
                ELSE
                       VQ(IXY          ) = VQ(IXY          )          &
                                            + AQ(IXY) * (1.-RD3)*RD4/6.
                END IF
!
            END IF
          END DO
!
! 3.a.3 Restore boundary data
!
        DO IXY=1, NX*NY
          IF ( MAPSTA(IXY).EQ.2 ) VQ(IXY) = AQ(IXY)
          END DO
!
! 3.a.4 Global closure
!
        IF ( GLOBAL ) THEN
            DO IY=1, NY
              VQ(IY          ) = VQ(IY          ) + VQ(NX*NY+IY)
              VQ((NX-1)*NY+IY) = VQ((NX-1)*NY+IY) + VQ(IY-NY)
              END DO
          END IF
!
! 3.b Propagate fields
!
        IF ( YFIRST ) THEN
            IF ( FLCY ) CALL W3QCK3                                   &
                        (NX, NY, NX, NY, VLCFLY, ATRNY, VQ,           &
                         .FALSE., 1, MAPAXY, NACT, MAPY2, NMY0,       &
                         NMY1, NMY2, NDSE, NDST )
            IF ( FLCX ) CALL W3QCK3                                   &
                        (NX, NY, NX, NY, VLCFLX, ATRNX, VQ,           &
                         GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0,       &
                         NMX1, NMX2, NDSE, NDST )
          ELSE
            IF ( FLCX ) CALL W3QCK3                                   &
                        (NX, NY, NX, NY, VLCFLX, ATRNX, VQ,           &
                         GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0,       &
                         NMX1, NMX2, NDSE, NDST )
            IF ( FLCY ) CALL W3QCK3                                   &
                        (NX, NY, NX, NY, VLCFLY, ATRNY, VQ,           &
                         .FALSE., 1, MAPAXY, NACT, MAPY2, NMY0,       &
                         NMY1, NMY2, NDSE, NDST )
          END IF
!
! 3.c Update boundaries
!
        IF ( FLBPI ) THEN
            RD1    = DSEC21 ( TBPI0, TIME ) - DTG *                   &
                                      REAL(NTLOC-ITLOC)/REAL(NTLOC)
            RD2    = DSEC21 ( TBPI0, TBPIN )
            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
              IXY     = MAPSF(ISBPI(IBI),3)
              VQ(IXY) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) )   &
                          / CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
              END DO
          END IF
!
        YFIRST = .NOT. YFIRST
        END DO
!
! 4.  Store results in VQ in proper format --------------------------- *
!
      DO ISEA=1, NSEA
        IXY    = MAPSF(ISEA,3)
        IF ( MAPSTA(IXY) .GT. 0 ) THEN
            VQ(IXY) =  MAX ( 0. , CG(IK,ISEA)/CLATS(ISEA)*VQ(IXY) )
          END IF
        END DO
!
      RETURN
!
! Formats
!
!/
!/ End of W3XYP3 ----------------------------------------------------- /
!/
      END SUBROUTINE W3XYP3
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3KTP3 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH,    &
                          DDDX, DDDY, CX, CY, DCXDX, DCXDY,           &
                          DCYDX, DCYDY, VA )
!/
!/    *** THIS ROUTINE SHOULD BE IDENTICAL TO W3KTP2 ***
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         17-Dec-2004 |
!/                  +-----------------------------------+
!/
!/    14-Feb-2000 : Origination.                        ( version 2.08 )
!/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/
!  1. Purpose :
!
!     Propagation in spectral space.
!
!  2. Method :
!
!     Third order QUICKEST scheme with ULTIMATE limiter.
!
!     As with the spatial propagation, the two spaces are considered
!     independently, but the propagation is performed in a 2-D space.
!     Compared to the propagation in physical space, the directions
!     rerpesent a closed space and are therefore comparable to the
!     longitudinal or 'X' propagation. The wavenumber space has to be
!     extended to allow for boundary treatment. Using a simple first
!     order boundary treatment at both sided, two points need to
!     be added. This implies that the spectrum needs to be extended,
!     shifted and rotated, as is performed using MAPTH2 as set
!     in W3MAP3.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       ISEA    Int.   I   Number of sea point.
!       FACTH/K Real   I   Factor in propagation velocity.
!       CTHG0   Real   I   Factor in great circle refracftion term.
!       MAPxx2  I.A.   I   Propagation and storage maps.
!       CG      R.A.   I   Local group velocities.
!       WN      R.A.   I   Local wavenumbers.
!       DEPTH   R.A.   I   Depth.
!       DDDx    Real   I   Depth gradients.
!       CX/Y    Real   I   Current components.
!       DCxDx   Real   I   Current gradients.
!       VA      R.A.  I/O  Spectrum.
!     ----------------------------------------------------------------
!
!     Local variables.
!     ----------------------------------------------------------------
!       DSDD    R.A.  Partial derivative of sigma for depth.
!       FDD, FDU, FDG, FCD, FCU
!               R.A.  Directionally varying part of depth, current and
!                     great-circle refraction terms and of consit.
!                     of Ck term.
!       CFLT-K  R.A.  Propagation velocities of local fluxes.
!       DB      R.A.  Wavenumber band widths at cell centers.
!       DM      R.A.  Wavenumber band widths between cell centers and
!                     next cell center.
!       Q       R.A.  Extracted spectrum
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!       W3QCK1   Actual propagation routine.
!       W3QCK2   Actual propagation routine.
!       STRACE   Service routine.
!
!  5. Called by :
!
!       W3WAVE   Wave model routine.
!
!  6. Error messages :
!
!       None.
!
!  8. Structure :
!
!     -----------------------------------------------------------------
!       1.  Preparations
!         a Initialize arrays
!         b Set constants and counters
!       2.  Point  preparations
!         a Calculate DSDD
!         b Extract spectrum
!       3.  Refraction velocities
!         a Filter level depth reffraction.
!         b Depth refratcion velocity.
!         c Current refraction velocity.
!       4.  Wavenumber shift velocities
!         a Prepare directional arrays
!         b Calcuate velocity.
!       5.  Propagate.
!       6.  Store results.
!     -----------------------------------------------------------------
!
!  9. Switches :
!
!       !/S     Enable subroutine tracing.
!       !/T     Enable general test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3CONSTANTS
      USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, SIG, DSIP, ECOS, ESIN, &
                          EC2, ESC, ES2, FACHFA, MAPWN, FLCTH, FLCK,  &
                          CTMAX
      USE W3ADATMD, ONLY: MAPTH2, MAPWN2, ITIME
      USE W3IDATMD, ONLY: FLCUR
      USE W3ODATMD, ONLY: NDSE, NDST
      USE W3UQCKMD
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: ISEA
      REAL, INTENT(IN)        :: FACTH, FACK, CTHG0, CG(0:NK+1),      &
                                 WN(0:NK+1), DEPTH, DDDX, DDDY,       &
                                 CX, CY, DCXDX, DCXDY, DCYDX, DCYDY
      REAL, INTENT(INOUT)     :: VA(NSPEC)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: ITH, IK, ISP
      REAL                    :: FDDMAX, FDG, FKD, FKD0, DCYX,        &
                                 DCXXYY, DCXY, DCXX, DCXYYX, DCYY
      REAL                    :: DSDD(0:NK+1), FRK(NK), FRG(NK),      &
                                 FKC(NTH), VQ(-NK-1:NK2*(NTH+2)),     &
                                 DB(NK2,NTH+1), DM(NK2,0:NTH+1),      &
                                 VCFLT(NK2*(NTH+1)), CFLK(NK2,NTH)
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Preparations --------------------------------------------------- *
! 1.a Initialize arrays
!
      IF ( FLCK ) VQ = 0.
!
! 2.  Preparation for point ------------------------------------------ *
! 2.a Array with partial derivative of sigma versus depth
!
      DO IK=0, NK+1
        IF ( DEPTH*WN(IK) .LT. 5. ) THEN
            DSDD(IK) = MAX ( 0. ,                                     &
                CG(IK)*WN(IK)-0.5*SIG(IK) ) / DEPTH
          ELSE
            DSDD(IK) = 0.
          END IF
        END DO
!
! 2.b Extract spectrum
!
      DO ISP=1, NSPEC
        VQ(MAPTH2(ISP)) = VA(ISP)
        END DO
!
! 3.  Refraction velocities ------------------------------------------ *
!
      IF ( FLCTH ) THEN
!
! 3.a Set slope filter for depth refraction
!
          FDDMAX = 0.
          FDG    = FACTH * CTHG0
!
          DO ITH=1, NTH/2
            FDDMAX = MAX(FDDMAX,ABS(ESIN(ITH)*DDDX-ECOS(ITH)*DDDY))
            END DO
!
          DO IK=1, NK
            FRK(IK) = FACTH * DSDD(IK) / WN(IK)
            FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX )
            FRG(IK) = FDG * CG(IK)
            END DO
!
! 3.b Depth refraction and great-circle propagation
!
          DO ISP=1, NSPEC
            VCFLT(MAPTH2(ISP)) = FRG(MAPWN(ISP)) * ECOS(ISP)          &
             + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DDDX - ECOS(ISP)*DDDY )
            END DO
!
! 3.c Current refraction
!
          IF ( FLCUR ) THEN
!
              DCYX   = FACTH *   DCYDX
              DCXXYY = FACTH * ( DCXDX - DCYDY )
              DCXY   = FACTH *   DCXDY
!
              DO ISP=1, NSPEC
                VCFLT(MAPTH2(ISP)) = VCFLT(MAPTH2(ISP)) +             &
                  ES2(ISP)*DCYX  + ESC(ISP)*DCXXYY - EC2(ISP)*DCXY
                  END DO
!
            END IF
!
        END IF
!
! 4.  Wavenumber shift velocities ------------------------------------ *
!     FACK is just the time step, which is accounted for in W3QCK2
!
      IF ( FLCK ) THEN
!
! 4.a Directionally dependent part
!
          DCXX   =  -   DCXDX
          DCXYYX =  - ( DCXDY + DCYDX )
          DCYY   =  -   DCYDY
          FKD    =    ( CX*DDDX + CY*DDDY )
!
          DO ITH=1, NTH
            FKC(ITH) = EC2(ITH)*DCXX +                                &
                       ESC(ITH)*DCXYYX + ES2(ITH)*DCYY
            END DO
!
! 4.b Velocities
!
          DO IK=0, NK+1
            FKD0   = FKD / CG(IK) * DSDD(IK)
            DO ITH=1, NTH
              CFLK(IK+1,ITH) = FKD0 + WN(IK)*FKC(ITH)
              END DO
            END DO
!
! 4.c Band widths
!
          DO IK=0, NK
            DB(IK+1,1) = DSIP(IK) / CG(IK)
            DM(IK+1,1) = WN(IK+1) - WN(IK)
            END DO
          DB(NK+2,1) = DSIP(NK+1) / CG(NK+1)
          DM(NK+2,1) = 0.
!
          DO ITH=2, NTH
            DO IK=1, NK+2
              DB(IK,ITH) = DB(IK,1)
              DM(IK,ITH) = DM(IK,1)
              END DO
            END DO
!
        END IF
!
! 5.  Propagate ------------------------------------------------------ *
!
      IF ( MOD(ITIME,2) .EQ. 0 ) THEN
          IF ( FLCK ) THEN
              DO ITH=1, NTH
                VQ(NK+2+(ITH-1)*NK2) = FACHFA * VQ(NK+1+(ITH-1)*NK2)
                END DO
              CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM,   &
                            VQ, .FALSE., 1, MAPTH2, NSPEC,            &
                            MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH,      &
                            NDSE, NDST )
            END IF
          IF ( FLCTH ) THEN
              CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE.,    &
                            NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, &
                            NSPEC, NDSE, NDST )
            END IF
        ELSE
          IF ( FLCTH ) THEN
              CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE.,    &
                            NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, &
                            NSPEC, NDSE, NDST )
            END IF
          IF ( FLCK )  THEN
              DO ITH=1, NTH
                VQ(NK+2+(ITH-1)*NK2) = FACHFA * VQ(NK+1+(ITH-1)*NK2)
                END DO
              CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM,   &
                            VQ, .FALSE., 1, MAPTH2, NSPEC,            &
                            MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH,      &
                            NDSE, NDST )
          END IF
        END IF
!
! 6.  Store reults --------------------------------------------------- *
!
      DO ISP=1, NSPEC
        VA(ISP) = VQ(MAPTH2(ISP))
        END DO
!
      RETURN
!
! Formats
!
!/
!/ End of W3KTP3 ----------------------------------------------------- /
!/
      END SUBROUTINE W3KTP3
!/
!/ End of module W3PRO3MD -------------------------------------------- /
!/
      END MODULE W3PRO3MD

max_stars_repo_name: cbegeman/OceanNGD

max_stars_repo_licenses: zlib-acknowledgement RSA-MD FTL

max_issues_repo_name: cbegeman/OceanNGD

max_issues_repo_licenses: zlib-acknowledgement RSA-MD FTL

max_forks_repo_name: cbegeman/OceanNGD

max_forks_repo_licenses: zlib-acknowledgement RSA-MD FTL