text:
subroutine z_hydpres_nhfull(mmax      ,nmax      ,j         ,nmmaxj    ,nmmax     , &
                          & kmax      ,nst       ,icx       ,icy       ,nsrc      , &
                          & norow     ,nocol     ,irocol    ,s1        ,s0        , &
                          & u1        ,v1        ,w1        ,guu       ,gvv       , &
                          & u0        ,v0        ,w0        ,circ2d    ,            &
                          & guv       ,gvu       ,gsqs      ,guz       ,gvz       , &
                          & kfu       ,kfv       ,kfs       ,kcs       ,kfuz0     , &
                          & kfvz0     ,kfsz0     ,kfsmin    ,kfsmx0    ,kcshyd    , &
                          & mnksrc    ,umean     ,vmean     ,dps       , &
                          & dzu0      ,dzv0      ,dzs0      ,disch     , &
                          & p1        ,p0        ,pnhcor    ,aak       ,zk        , &
                          & bbk       ,cck       ,ddk       ,bbka      ,bbkc      , &
                          & aak2      ,cck2      ,dinv      ,pbbk      ,pbbkc     , &
                          & vj        ,rj        ,pj        ,rjshadow  ,sj        , &
                          & sjprec    ,tj        ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: z_hydpres_nhfull.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/non_hydro/z_hydpres_nhfull.f90 $
!!--description-----------------------------------------------------------------
!
! The difference equations for the non-hydrostatic
! pressure are derived from the discretized
! incompressibility condition.
! The non-hydrostatic pressure is computed for
! a restricted horizontal area.
! At the open boundaries the discretization
! is such that it leads to a symmetric matrix.
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
   !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp)               , pointer :: hdt
    real(fp)               , pointer :: rhow
    real(fp)               , pointer :: ag
    real(fp)               , pointer :: tetaq
    real(fp)               , pointer :: tetaz
    integer                , pointer :: iro
    integer                , pointer :: m1_nhy
    integer                , pointer :: m2_nhy
    integer                , pointer :: n1_nhy
    integer                , pointer :: n2_nhy
    logical                , pointer :: flag_pp
    character(8)           , pointer :: krylov
    integer                , pointer :: lundia
!
! Global variables
!
    integer                                                        :: icx
    integer                                                        :: icy
    integer                                                        :: j
    integer                                                        :: kmax      !  Description and declaration in esm_alloc_int.f90
    integer                                                        :: mmax      !  Description and declaration in esm_alloc_int.f90
    integer                                                        :: nmax      !  Description and declaration in esm_alloc_int.f90
    integer                                                        :: nmmax     !  Description and declaration in dimens.igs
    integer                                                        :: nmmaxj    !  Description and declaration in dimens.igs
    integer                                           , intent(in) :: nocol     !  Description and declaration in esm_alloc_int.f90
    integer                                           , intent(in) :: norow     !  Description and declaration in esm_alloc_int.f90
    integer                                           , intent(in) :: nsrc      !  Description and declaration in esm_alloc_int.f90
    integer                                                        :: nst
    integer, dimension(5, norow + nocol)              , intent(in) :: irocol    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(7, nsrc)                       , intent(in) :: mnksrc    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)         , intent(in) :: kcshyd    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)         , intent(in) :: kcs       !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)                      :: kfs       !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)         , intent(in) :: kfsmin    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)         , intent(in) :: kfsmx0    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)         , intent(in) :: kfu       !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)         , intent(in) :: kfv       !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                :: kfsz0     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)   , intent(in) :: kfuz0     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)   , intent(in) :: kfvz0     !  Description and declaration in esm_alloc_int.f90
    real(prec), dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in) :: dps       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: gsqs      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: guu       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: guv       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: gvu       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: gvv       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: guz       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: gvz       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)        , intent(in) :: s0        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                     :: s1        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                     :: umean     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                     :: vmean     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax), intent(in) :: w1        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: aak
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: aak2
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: vj        ! Renamed apj into v for use in BiCGSTAB, Ullmann 21/02/2008
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: bbk
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: bbka
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: bbkc
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: cck
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: cck2
    real(fp), dimension(4, norow+nocol)               , intent(in) :: circ2d    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: ddk
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: dinv
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: dzs0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: dzu0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: dzv0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: p0        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: p1        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: pbbk
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: pbbkc
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: pj
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: rjshadow  ! Extra array for BiCGSTAB, Ullmann 21/02/2008
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: sj        ! Extra array for BiCGSTAB, Ullmann 21/02/2008
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: sjprec    ! Extra array for BiCGSTAB, Ullmann 21/02/2008
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: tj        ! Extra array for BiCGSTAB, Ullmann 21/02/2008
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: pnhcor    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)               :: rj
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: u1        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: v1        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: u0        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)  , intent(in) :: v0        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax), intent(in) :: w0        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nsrc)                         , intent(in) :: disch     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax)                       , intent(in) :: zk        !  Description and declaration in esm_alloc_real.f90
!
! Local variables
!
    integer       :: ddb
    integer       :: i
    integer       :: ibf
    integer       :: ibl
    integer       :: ic
    integer       :: icxy
    integer       :: k
    integer       :: kdo
    integer       :: kenm
    integer       :: kk
    integer       :: kup
    integer       :: m
    integer       :: mf
    integer       :: mfu
    integer       :: ml
    integer       :: n
    integer       :: ndelta
    integer       :: ndm
    integer       :: nf
    integer       :: nfm
    integer       :: nfu
    integer       :: nfum
    integer       :: nfdm
    integer       :: nl
    integer       :: nlm
    integer       :: nlum
    integer       :: nm
    integer       :: nmd
    integer       :: nmf
    integer       :: nmfu
    integer       :: nml
    integer       :: nmlu
    integer       :: nmrow
    integer       :: nmst
    integer       :: nmstart
    integer       :: nmu
    integer       :: num
    real(fp)      :: corr
    real(fp)      :: dt
    real(fp)      :: dxd
    real(fp)      :: dxu
    real(fp)      :: dyd
    real(fp)      :: dyu
    real(fp)      :: dzdo
    real(fp)      :: dzup
    real(fp)      :: dz
    real(fp)      :: alfa
    real(fp)      :: alfai
    real(fp)      :: alfae
    real(fp)      :: coef
    real(fp)      :: alfak
    real(fp)      :: alfakm1
    real(fp)      :: z_k
    real(fp)      :: z_km1
    character(6)  :: method
    character(20) :: errtxt
!
!! executable statements -------------------------------------------------------
!
    lundia    => gdp%gdinout%lundia
    m1_nhy   => gdp%gdnonhyd%m1_nhy
    m2_nhy   => gdp%gdnonhyd%m2_nhy
    n1_nhy   => gdp%gdnonhyd%n1_nhy
    n2_nhy   => gdp%gdnonhyd%n2_nhy
    flag_pp  => gdp%gdnonhyd%flag_pp
    tetaq    => gdp%gdnonhyd%tetaq
    tetaz    => gdp%gdnonhyd%tetaz
    krylov   => gdp%gdnonhyd%krylov
    rhow     => gdp%gdphysco%rhow
    ag       => gdp%gdphysco%ag
    iro      => gdp%gdphysco%iro
    hdt      => gdp%gdnumeco%hdt
    !
    ddb     = gdp%d%ddbound
    icxy    = max(icx, icy)
    dt      = 2.0_fp * hdt
    ndelta  = n2_nhy - n1_nhy
    nmstart = (n1_nhy+ddb) + (m1_nhy-1+ddb)*icxy
    !
    !
    !
    ! initialise all coefficients
    !
    aak  = 0.0_fp
    aak2 = 0.0_fp
    bbk  = 1.0_fp
    cck  = 0.0_fp
    cck2 = 0.0_fp
    ddk  = 0.0_fp
    bbka = 0.0_fp
    bbkc = 0.0_fp
    p1   = 0.0_fp
    !
    ! inner area, closed boundaries: homogeneous Neumann horizontal gradients of velocity
    ! coefficients:
    !
    !   horizontal:             cck2(  m,n+1)                  vertical:  bbkc(k+1)
    !                            |                                         |
    !            aak(m-1,  n) - bbk (  m,  n) - cck (m+1,n)               bbk (  k)
    !                                   |                   |
    !                           aak2(  m,n-1)                             bbka(k-1)
    !

    do nm = 1,nmmax
       if (kcs(nm)*kfs(nm) == 1) then
          nmd = nm - icx
          nmu = nm + icx
          ndm = nm - icy
          num = nm + icy
          do k = kfsmin(nm), kfsmx0(nm)
             !
             ! Construction of full system of equations for pressure
             ! loop over all cells where continuity equation is applied, i.e., over
             ! all (partially) filled cells, which includes bottom cells and surface cells
             !
             !    Construction includes boundary conditions at vertical faces (boundaries left/right,
             !    vertical bottom walls, etc.) by means of using switches kf...
             !
             kdo  = max( k-1, kfsmin(nm) )
             kup  = min( k+1, kfsmx0(nm) )
             dxu  = gvu(nm)
             dxd  = gvu(nmd)
             dyu  = guv(nm)
             dyd  = guv(ndm)
             dzup        = 0.5_fp * (zk(k)-zk(k-1)+zk(kup)-zk(kup-1))
             dzdo        = 0.5_fp * (zk(k)-zk(k-1)+zk(kdo)-zk(kdo-1))
             aak2(nm, k) = -  dt * dzv0(ndm, k)*gvv(ndm)/ (dyd*rhow)       &
                         &       * kfvz0(ndm, k) * kfsz0(ndm, k)
             cck2(nm, k) = -  dt * dzv0(nm , k)*gvv(nm) / (dyu*rhow)       &
                         &       * kfvz0(nm , k) * kfsz0(num, k)
             bbka(nm, k) = -  dt *gsqs(nm)/ (dzdo*rhow)
             bbkc(nm, k) = -  dt *gsqs(nm)/ (dzup*rhow)
             aak (nm, k) = -  dt * dzu0(nmd, k)*guu(nmd) / (dxd*rhow)      &
                         &       * kfuz0(nmd, k) * kfsz0(nmd, k)
             cck (nm, k) = -  dt * dzu0( nm, k)*guu(nm ) / (dxu*rhow)      &
                         &       * kfuz0(nm, k) * kfsz0(nmu, k)
             bbk (nm, k) = - aak (nm, k) - cck (nm, k) - aak2(nm, k)       &
                         & - cck2(nm, k) - bbka(nm, k) - bbkc(nm, k)
             ddk (nm, k) = - ( u1(nm , k)*dzu0(nm , k)*guu(nm )            &
                         &   - u1(nmd, k)*dzu0(nmd, k)*guu(nmd) )          &
                         & - ( v1(nm , k)*dzv0(nm , k)*gvv(nm )            &
                         &   - v1(ndm, k)*dzv0(ndm, k)*gvv(ndm) )          &
                         & - w1(nm, k)*gsqs(nm) + w1(nm, k - 1)*gsqs(nm)
             !
             ! If outflow at a vertical free-surface wall then hor.mom.eq. applied => horizontal
             ! pressure gradient required, but no pressure outside flow domain => increase only
             ! diagonal of system of equations
             !
             !coef = 0.0_fp
             !if ( kfvz0(ndm,k)==1 .and. kfsz0(ndm,k)==0 ) then
             !   coef = coef + dt * dzv0(ndm, k)*gvv(ndm) / (dyd*rhow)
             !endif
             !if ( kfvz0(nm ,k)==1 .and. kfsz0(num,k)==0 ) then
             !   coef = coef + dt * dzv0(nm , k)*gvv(nm ) / (dyu*rhow)
             !endif
             !if ( kfuz0(nmd,k)==1 .and. kfsz0(nmd,k)==0 ) then
             !   coef = coef + dt * dzu0(nmd, k)*guu(nmd) / (dxd*rhow)
             !endif
             !if ( kfuz0(nm, k)==1 .and. kfsz0(nmu,k)==0 ) then
             !   coef = coef + dt * dzu0(nm , k)*guu(nm ) / (dxu*rhow)
             !endif
             !bbk(nm,k) = bbk(nm,k) + coef
          enddo
       endif
    enddo
    !
    ! Boundary conditions at free surface and at bottom
    ! First the free surface, then the bottom, because the free-surface procedure uses pressure
    ! in 2 grid layers, which in case of flow in only 1 layer are put together in the bottom
    ! procedure.
    !
    do nm = 1,nmmax
       nmu = nm + icx
       num = nm + icy
       nmd = nm - icx
       ndm = nm - icy
       if (kcs(nm)*kfs(nm) == 1) then
          !
          k = kfsmx0(nm)
          !
          ! surface: insert dp/dz=0 at surface and add dzeta/dt effect
          !
          bbk (nm,k) = bbk(nm,k) + bbkc(nm,k)
          bbkc(nm,k) = 0.0_fp
          !
          if ( flag_pp ) then
             !
             ! Non-hydrostatic pressure at free surface equal to zero
             !
             z_k = 0.5_fp * ( zk(k) + zk(k-1) )
             if ( k > 2 ) then
                !
                ! interpolation possible
                !
                z_km1 = 0.5_fp * ( zk(k-1) + zk(k-2) )
             else
                !
                ! extrapolation
                !
                z_km1 = 1.5_fp * zk(k-1) - 0.5_fp * zk(k)
             endif
             coef = gsqs(nm) / dt / tetaq / tetaz                                       &
                  / ( (z_k - z_km1) *ag*rhow ) ! - p0(nm,k) + p0(nm,k-1) ) ! extra linearization
             !
             ! Correction for neglecting dq/dz at the free surface
             !
             corr    = 1.0_fp
             !corr   = 1.0_fp / (1.0_fp-((p0(nm,k)-p0(nm,k-1))/(ag*rhow*(z_k-z_km1))))
             !
             alfak   = coef * (s0(nm) - z_km1 ) * corr
             alfakm1 = coef * (z_k    - s0(nm)) * corr
          else
             !
             ! Top non-hydrostatic pressure equal to zero => hydrostatic pressure in top layer
             !
             alfak   = gsqs(nm) / dt / tetaq / tetaz / ( ag*rhow )
             alfakm1 = 0.0_fp
          endif
          !
          bbk (nm, k) = bbk (nm,k) + alfak
          bbka(nm, k) = bbka(nm,k) + alfakm1
          !
          ddk (nm, k) = ddk(nm,k)     + w1(nm,k)*gsqs(nm)     & ! Should not be necessary
                      - (1.0_fp-tetaz)/tetaz *                                        &
                                      (  (  u0(nm , k)*dzu0(nm , k)*guu(nm )          &
                                          - u0(nmd, k)*dzu0(nmd, k)*guu(nmd) )        &
                                       + (  v0(nm , k)*dzv0(nm , k)*gvv(nm )          &
                                          - v0(ndm, k)*dzv0(ndm, k)*gvv(ndm) )        &
                                       - w0(nm, k - 1)*gsqs(nm)                )      &
                      - alfak * p0(nm,k) !!! + ag*rhow*(s0(nm)-z_k  ) ) )
          if (k > kfsmin(nm)) then
             ddk (nm, k) = ddk(nm,k) - alfakm1 * p0(nm,k-1) !!! + ag*rhow*(s0(nm)-z_km1) )
          endif
          !
          k = kfsmin(nm)
          !
          ! bottom: insert dp/dz = 0 at bottom
          !
          ! dp/dn = 0 at vertical bottom walls as well. Already set above because those sides
          ! are closed => switches kfuz1 etc. set equal to 0 => coefficients cck, etc. set to 0
          !
          bbk (nm,k) = bbk(nm,k) + bbka(nm,k)
          bbka(nm,k) = 0.0_fp
       endif
    enddo
    !
    ! addition of discharges to continuity equation only when it lies inside NH domain
    !
    do i = 1, nsrc
       nm   = (mnksrc(2,i)+ddb)+((mnksrc(1,i)-1)+ddb) * icxy
       if (kcs(nm) == 1) then
          k    = mnksrc(3, i)
          kenm = min(1, kfu(nm)+kfu(nm-icx)+kfv(nm)+kfv(nm-icy))
          if (kenm/=0 .or. disch(i)>=0.0_fp) then
             if (k/=0 .or. kfsmx0(nm)<=kfsmin(nm)) then
                !
                ! The order is important at dry points
                !
                if (k > kfsmx0(nm)) then
                   k = kfsmx0(nm)
                endif
                if (k < kfsmin(nm)) then
                   k = kfsmin(nm)
                endif
                ddk(nm,k) = ddk(nm,k) + disch(i)/gsqs(nm)
             else
                do kk = kfsmin(nm), kfsmx0(nm)
                   ddk(nm,kk) = ddk(nm,kk) + disch(i)*dzs0(nm,kk)                &
                              & /((real(dps(nm),fp)+s0(nm))*gsqs(nm))
                enddo
             endif
          endif
       endif
       !
       ! in case of an intake for an intake/outfall combination:
       !
       if (mnksrc(7,i) >= 2) then
          nm   = (mnksrc(2,i)+ddb)+((mnksrc(1,i)-1)+ddb) * icxy
          k    = mnksrc(3,i)
          kenm = min(1, kfu(nm)+kfu(nm-icx)+kfv(nm)+kfv(nm-icy))
          if (kenm/=0 .or. -disch(i)>=0.0_fp) then
             if (k/=0 .or. kfsmx0(nm)<=kfsmin(nm)) then
                !
                ! The order is important at dry points
                !
                if (k > kfsmx0(nm)) then
                   k = kfsmx0(nm)
                endif
                if (k < kfsmin(nm)) then
                   k = kfsmin(nm)
                endif
                ddk(nm,k) = ddk(nm,k) - disch(i)/gsqs(nm)
             else
                do kk = kfsmin(nm), kfsmx0(nm)
                      ddk(nm,kk) = ddk(nm,kk) - disch(i)*dzs0(nm,kk)            &
                                 & /((real(dps(nm),fp)+s0(nm))*gsqs(nm))
                enddo
             endif
          elseif (mnksrc(7,i) /= 3) then
             !
             ! in case of a culvert no warning generated
             !
             write (errtxt, '(i0,i3)') nst, i
             call prterr(lundia, 'S208', trim(errtxt))
          else
          endif
       endif
    enddo
    !
    !
    ! Open boundary conditions
    ! At water level boundaries a Dirichlet bc. and for other other boundary types
    ! a Neumann bc are prescribed.
    !
    do ic = 1, norow
       n    = irocol(1, ic)
       mf   = irocol(2, ic) - 1
       ml   = irocol(3, ic)
       ibf  = irocol(4, ic)
       ibl  = irocol(5, ic)
       nmf  = (n+ddb)*icy + (mf+ddb)*icx - icxy
       nmfu = nmf + icx
       nml  = (n+ddb)*icy + (ml+ddb)*icx - icxy
       nmlu =  nml + icx
       if (ibf == 2) then
          do k = kfsmin(nmfu), kfsmx0(nmfu)
             ddk(nmfu,k) = ddk(nmfu,k) - aak(nmfu,k)*tetaq*ag*rhow*(circ2d(1,ic)-s0(nmf))
             aak(nmfu,k) = 0.0_fp
          enddo
       endif
       if (ibf > 2) then
          do k = kfsmin(nmfu), kfsmx0(nmfu)
             bbk(nmfu,k) = bbk(nmfu,k) + aak(nmfu,k)
             aak(nmfu,k) = 0.0_fp
          enddo
       endif
       if (ibl == 2) then
          do k = kfsmin(nml), kfsmx0(nml)
             ddk(nml,k) = ddk(nml,k) - cck(nml,k)*tetaq*ag*rhow*(circ2d(2,ic)-s0(nmlu))
             cck(nml,k) = 0.0_fp
          enddo
       endif
       if (ibl > 2) then
          do k = kfsmin(nml), kfsmx0(nml)
             bbk(nml,k) = bbk(nml,k) + cck(nml,k)
             cck(nml,k) = 0.0_fp
          enddo
       endif
    enddo
    do ic = norow + 1, norow + nocol
       m    = irocol(1, ic)
       nf   = irocol(2, ic) - 1
       nl   = irocol(3, ic)
       ibf  = irocol(4, ic)
       ibl  = irocol(5, ic)
       nfm  = (nf+ddb)*icy + (m+ddb)*icx - icxy
       nfum = nfm + icy
       nlm  = (nl+ddb)*icy + (m+ddb)*icx - icxy
       nlum = nlm + icy
       if (ibf == 2) then
          do k = kfsmin(nfum), kfsmx0(nfum)
             ddk (nfum,k) = ddk(nfum, k)- aak2(nfum,k)*tetaq*ag*rhow*(circ2d(1,ic)-s0(nfm))
             aak2(nfum,k) = 0.0_fp
          enddo
       endif
       if (ibf > 2) then
          do k = kfsmin(nfum), kfsmx0(nfum)
             bbk (nfum,k) = bbk(nfum, k) + aak2(nfum, k)
             aak2(nfum,k) = 0.0_fp
          enddo
       endif
       if (ibl == 2) then
          do k = kfsmin(nlm), kfsmx0(nlm)
             ddk (nlm,k) = ddk(nlm, k)-cck2(nlm, k)*tetaq*ag*rhow*(circ2d(2,ic)-s0(nlum))
             cck2(nlm,k) = 0.0_fp
          enddo
       endif
       if (ibl > 2) then
          do k = kfsmin(nlm), kfsmx0(nlm)
             bbk (nlm,k)  = bbk(nlm,k) + cck2(nlm,k)
             cck2(nlm,k) = 0.0_fp
          enddo
       endif
    enddo
    !
    ! solve system of equations for pressure correction
    ! Now p0 as a better initial guess
    !
    call z_initcg_nhfull(aak       ,bbk       ,cck       ,aak2      ,cck2      , &
                       & bbka      ,bbkc      ,ddk       ,kmax      ,icx       , &
                       & icy       ,nmmax     ,nst       ,kfsz0     ,dinv      , &
                       & pbbk      ,pbbkc     ,pnhcor    ,rj        ,kfs       , &
                       & kfsmin    ,kfsmx0    ,gdp       )
    if ( krylov == 'bicgstab') then
        call z_solbicgstab(aak       ,bbk       ,cck      ,aak2      ,cck2      , &
                         & bbka      ,bbkc      ,ddk      ,kmax      ,icx       , &
                         & icy       ,nmmax     ,nst      ,kfsz0     ,pnhcor    , &
                         & pj        ,rj        ,vj       ,dinv      ,pbbk      , &
                         & pbbkc     ,p1        ,rjshadow ,sj        ,sjprec    , &
                         & tj        ,kfs       ,kfsmin   ,kfsmx0    ,gdp       )

    elseif ( krylov == 'cg') then
       !
       ! old implementation (see version nh-15) and rename "apj" by "vj"
       !
       call z_solcg_nhfull(aak       ,bbk       ,cck       ,aak2      ,cck2      , &
                         & bbka      ,bbkc      ,ddk       ,kmax      ,icx       , &
                         & icy       ,nmmax     ,nst       ,kfsz0     ,pnhcor    , &
                         & pj        ,rj        ,vj        ,dinv      ,pbbk      , &
                         & pbbkc     ,p1        ,kfs       ,kfsmin    ,kfsmx0    , &
                         & gdp       )
    else
        write (lundia,*) '*** ERROR:z_hydpres_nhfull: THIS IS NOT ALLOWED'
    endif
    !
    !   pressure at open boundaries
    !
    do ic = 1, norow
       n    = irocol(1, ic)
       mf   = irocol(2, ic) - 1
       ml   = irocol(3, ic)
       ibf  = irocol(4, ic)
       ibl  = irocol(5, ic)
       nmf  = (n+ddb)*icy + (mf+ddb)*icx - icxy
       nmfu = nmf + icx
       nml  = (n+ddb)*icy + (ml+ddb)*icx - icxy
       nmlu = nml + icx
       if (ibf == 2) then
          s1(nmf) = circ2d(1,ic)
          do k = kfsmin(nmf), kfsmx0(nmf)
             p1(nmf,k) = tetaq * ag * rhow * (s1(nmf)-s0(nmf))
             p0(nmf,k) = 0.0_fp
          enddo
       endif
       if (ibl == 2) then
          s1(nmlu) = circ2d(2,ic)
          do k = kfsmin(nmlu), kfsmx0(nmlu)
             p1(nmlu,k) = tetaq * ag * rhow * (s1(nmlu)-s0(nmlu))
             p0(nmlu,k) = 0.0_fp
          enddo
       endif
    enddo
    do ic = norow + 1, norow + nocol
       m    = irocol(1, ic)
       nf   = irocol(2, ic) - 1
       nl   = irocol(3, ic)
       ibf  = irocol(4, ic)
       ibl  = irocol(5, ic)
       nfm  = (nf+ddb)*icy + (m+ddb)*icx - icxy
       nlm  = (nl+ddb)*icy + (m+ddb)*icx - icxy
       nfum = nfm + icy
       nlum = nlm + icy
       if (ibf == 2) then
          s1(nfm) = circ2d(1,ic)
          do k = kfsmin(nfm), kfsmx0(nfm)
             p1(nfm,k) = tetaq * ag * rhow * (s1(nfm)-s0(nfm))
             p0(nfm,k) = 0.0_fp
          enddo
       endif
       if (ibl == 2) then
          s1(nlum) = circ2d(2,ic)
          do k = kfsmin(nlum), kfsmx0(nlum)
             p1(nlum,k) = tetaq * ag * rhow * (s1(nlum)-s0(nlum))
             p0(nlum,k) = 0.0_fp
          enddo
       endif
    enddo
end subroutine z_hydpres_nhfull


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006, 2010, 2015, 2017
!        Salvatore Filippone    Cranfield University
!        Alfredo Buttari        CNRS-IRIT, Toulouse
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine psi_iaxpby(m,n,alpha, x, beta, y, info)

  use psb_const_mod
  use psb_error_mod
  implicit none
  integer(psb_ipk_), intent(in)      :: m, n
  integer(psb_ipk_), intent (in)       ::  x(:,:)
  integer(psb_ipk_), intent (inout)    ::  y(:,:)
  integer(psb_ipk_), intent (in)       ::  alpha, beta
  integer(psb_ipk_), intent(out)     :: info
  integer(psb_ipk_) :: err_act
  integer(psb_ipk_) :: lx, ly
  integer(psb_ipk_) :: ierr(5)
  character(len=20)        :: name, ch_err

  name='psb_geaxpby'
  if(psb_get_errstatus() /= 0) return
  info=psb_success_
  call psb_erractionsave(err_act)

  if (m < 0) then
    info = psb_err_iarg_neg_
    ierr(1) = 1; ierr(2) = m
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if
  if (n < 0) then
    info = psb_err_iarg_neg_
    ierr(1) = 2; ierr(2) = n
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if
  lx = size(x,1)
  ly = size(y,1)
  if (lx < m) then
    info = psb_err_input_asize_small_i_
    ierr(1) = 4; ierr(2) = m
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if
  if (ly < m) then
    info = psb_err_input_asize_small_i_
    ierr(1) = 6; ierr(2) = m
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if

  if ((m>0).and.(n>0)) call iaxpby(m,n,alpha,x,lx,beta,y,ly,info)

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return
end subroutine psi_iaxpby

subroutine psi_iaxpbyv(m,alpha, x, beta, y, info)

  use psb_const_mod
  use psb_error_mod
  implicit none
  integer(psb_ipk_), intent(in)      :: m
  integer(psb_ipk_), intent (in)       ::  x(:)
  integer(psb_ipk_), intent (inout)    ::  y(:)
  integer(psb_ipk_), intent (in)       :: alpha, beta
  integer(psb_ipk_), intent(out)     :: info
  integer(psb_ipk_) :: err_act
  integer(psb_ipk_) :: lx, ly
  integer(psb_ipk_) :: ierr(5)
  character(len=20)        :: name, ch_err

  name='psb_geaxpby'
  if(psb_get_errstatus() /= 0) return
  info=psb_success_
  call psb_erractionsave(err_act)

  if (m < 0) then
    info = psb_err_iarg_neg_
    ierr(1) = 1; ierr(2) = m
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if
  lx = size(x,1)
  ly = size(y,1)
  if (lx < m) then
    info = psb_err_input_asize_small_i_
    ierr(1) = 3; ierr(2) = m
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if
  if (ly < m) then
    info = psb_err_input_asize_small_i_
    ierr(1) = 5; ierr(2) = m
    call psb_errpush(info,name,i_err=ierr)
    goto 9999
  end if

  if (m>0) call iaxpby(m,ione,alpha,x,lx,beta,y,ly,info)

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psi_iaxpbyv

subroutine psi_igthmv(n,k,idx,alpha,x,beta,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, k, idx(:)
  integer(psb_ipk_) :: x(:,:), y(:),alpha,beta

  ! Locals
  integer(psb_ipk_) :: i, j, pt

  if (beta == izero) then
    if (alpha == izero) then
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = izero
        end do
      end do
    else if (alpha == ione) then
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = x(idx(i),j)
        end do
      end do
    else if (alpha == -ione) then
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = -x(idx(i),j)
        end do
      end do
    else
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = alpha*x(idx(i),j)
        end do
      end do
    end if
  else
    if (beta == ione) then
      ! Do nothing
    else if (beta == -ione) then
      y(1:n*k) = -y(1:n*k)
    else
      y(1:n*k) = beta*y(1:n*k)
    end if

    if (alpha == izero) then
      ! do nothing
    else if (alpha == ione) then
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = y(pt) + x(idx(i),j)
        end do
      end do
    else if (alpha == -ione) then
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = y(pt) - x(idx(i),j)
        end do
      end do
    else
      pt=0
      do j=1,k
        do i=1,n
          pt=pt+1
          y(pt) = y(pt) + alpha*x(idx(i),j)
        end do
      end do
    end if
  end if

end subroutine psi_igthmv

subroutine psi_igthv(n,idx,alpha,x,beta,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, idx(:)
  integer(psb_ipk_) :: x(:), y(:),alpha,beta

  ! Locals
  integer(psb_ipk_) :: i
  if (beta == izero) then
    if (alpha == izero) then
      do i=1,n
        y(i) = izero
      end do
    else if (alpha == ione) then
      do i=1,n
        y(i) = x(idx(i))
      end do
    else if (alpha == -ione) then
      do i=1,n
        y(i) = -x(idx(i))
      end do
    else
      do i=1,n
        y(i) = alpha*x(idx(i))
      end do
    end if
  else
    if (beta == ione) then
      ! Do nothing
    else if (beta == -ione) then
      y(1:n) = -y(1:n)
    else
      y(1:n) = beta*y(1:n)
    end if

    if (alpha == izero) then
      ! do nothing
    else if (alpha == ione) then
      do i=1,n
        y(i) = y(i) + x(idx(i))
      end do
    else if (alpha == -ione) then
      do i=1,n
        y(i) = y(i) - x(idx(i))
      end do
    else
      do i=1,n
        y(i) = y(i) + alpha*x(idx(i))
      end do
    end if
  end if

end subroutine psi_igthv

subroutine psi_igthzmm(n,k,idx,x,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, k, idx(:)
  integer(psb_ipk_) :: x(:,:), y(:,:)

  ! Locals
  integer(psb_ipk_) :: i

  do i=1,n
    y(i,1:k)=x(idx(i),1:k)
  end do

end subroutine psi_igthzmm

subroutine psi_igthzmv(n,k,idx,x,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, k, idx(:)
  integer(psb_ipk_) :: x(:,:), y(:)

  ! Locals
  integer(psb_ipk_) :: i, j, pt

  pt=0
  do j=1,k
    do i=1,n
      pt=pt+1
      y(pt)=x(idx(i),j)
    end do
  end do

end subroutine psi_igthzmv

subroutine psi_igthzv(n,idx,x,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, idx(:)
  integer(psb_ipk_) :: x(:), y(:)

  ! Locals
  integer(psb_ipk_) :: i

  do i=1,n
    y(i)=x(idx(i))
  end do

end subroutine psi_igthzv

subroutine psi_isctmm(n,k,idx,x,beta,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, k, idx(:)
  integer(psb_ipk_) :: beta, x(:,:), y(:,:)

  ! Locals
  integer(psb_ipk_) :: i, j

  if (beta == izero) then
    do i=1,n
      y(idx(i),1:k) = x(i,1:k)
    end do
  else if (beta == ione) then
    do i=1,n
      y(idx(i),1:k) = y(idx(i),1:k)+x(i,1:k)
    end do
  else
    do i=1,n
      y(idx(i),1:k) = beta*y(idx(i),1:k)+x(i,1:k)
    end do
  end if
end subroutine psi_isctmm

subroutine psi_isctmv(n,k,idx,x,beta,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, k, idx(:)
  integer(psb_ipk_) :: beta, x(:), y(:,:)

  ! Locals
  integer(psb_ipk_) :: i, j, pt

  if (beta == izero) then
    pt=0
    do j=1,k
      do i=1,n
        pt=pt+1
        y(idx(i),j) = x(pt)
      end do
    end do
  else if (beta == ione) then
    pt=0
    do j=1,k
      do i=1,n
        pt=pt+1
        y(idx(i),j) = y(idx(i),j)+x(pt)
      end do
    end do
  else
    pt=0
    do j=1,k
      do i=1,n
        pt=pt+1
        y(idx(i),j) = beta*y(idx(i),j)+x(pt)
      end do
    end do
  end if
end subroutine psi_isctmv

subroutine psi_isctv(n,idx,x,beta,y)

  use psb_const_mod
  implicit none

  integer(psb_ipk_) :: n, idx(:)
  integer(psb_ipk_) :: beta, x(:), y(:)

  ! Locals
  integer(psb_ipk_) :: i

  if (beta == izero) then
    do i=1,n
      y(idx(i)) = x(i)
    end do
  else if (beta == ione) then
    do i=1,n
      y(idx(i)) = y(idx(i))+x(i)
    end do
  else
    do i=1,n
      y(idx(i)) = beta*y(idx(i))+x(i)
    end do
  end if
end subroutine psi_isctv
subroutine  iaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
  use psb_const_mod
  use psb_error_mod
  implicit none
  integer n, m, lldx, lldy, info
  integer(psb_ipk_) X(lldx,*), Y(lldy,*)
  integer(psb_ipk_) alpha, beta
  integer(psb_ipk_) :: i, j
  integer(psb_ipk_) :: int_err(5)
  character  name*20
  name='iaxpby'

  !
  !     Error handling
  !
  info = psb_success_
  if (m.lt.0) then
    info=psb_err_iarg_neg_
    int_err(1)=1
    int_err(2)=m
    call fcpsb_errpush(info,name,int_err)
    goto 9999
  else if (n.lt.0) then
    info=psb_err_iarg_neg_
    int_err(1)=1
    int_err(2)=n
    call fcpsb_errpush(info,name,int_err)
    goto 9999
  else if (lldx.lt.max(1,m)) then
    info=psb_err_iarg_not_gtia_ii_
    int_err(1)=5
    int_err(2)=1
    int_err(3)=lldx
    int_err(4)=m
    call fcpsb_errpush(info,name,int_err)
    goto 9999
  else if (lldy.lt.max(1,m)) then
    info=psb_err_iarg_not_gtia_ii_
    int_err(1)=8
    int_err(2)=1
    int_err(3)=lldy
    int_err(4)=m
    call fcpsb_errpush(info,name,int_err)
    goto 9999
  endif

  if (alpha.eq.izero) then
    if (beta.eq.izero) then
      do j=1, n
        do i=1,m
          y(i,j) = izero
        enddo
      enddo
    else if (beta.eq.ione) then
      !
      !     Do nothing!
      !

    else if (beta.eq.-ione) then
      do j=1,n
        do i=1,m
          y(i,j) = - y(i,j)
        enddo
      enddo
    else
      do j=1,n
        do i=1,m
          y(i,j) =  beta*y(i,j)
        enddo
      enddo
    endif

  else if (alpha.eq.ione) then

    if (beta.eq.izero) then
      do j=1,n
        do i=1,m
          y(i,j) = x(i,j)
        enddo
      enddo
    else if (beta.eq.ione) then
      do j=1,n
        do i=1,m
          y(i,j) = x(i,j) + y(i,j)
        enddo
      enddo

    else if (beta.eq.-ione) then
      do j=1,n
        do i=1,m
          y(i,j) = x(i,j) - y(i,j)
        enddo
      enddo
    else
      do j=1,n
        do i=1,m
          y(i,j) = x(i,j) + beta*y(i,j)
        enddo
      enddo
    endif

  else if (alpha.eq.-ione) then

    if (beta.eq.izero) then
      do j=1,n
        do i=1,m
          y(i,j) = -x(i,j)
        enddo
      enddo
    else if (beta.eq.ione) then
      do j=1,n
        do i=1,m
          y(i,j) = -x(i,j) + y(i,j)
        enddo
      enddo

    else if (beta.eq.-ione) then
      do j=1,n
        do i=1,m
          y(i,j) = -x(i,j) - y(i,j)
        enddo
      enddo
    else
      do j=1,n
        do i=1,m
          y(i,j) = -x(i,j) + beta*y(i,j)
        enddo
      enddo
    endif

  else

    if (beta.eq.izero) then
      do j=1,n
        do i=1,m
          y(i,j) = alpha*x(i,j)
        enddo
      enddo
    else if (beta.eq.ione) then
      do j=1,n
        do i=1,m
          y(i,j) = alpha*x(i,j) + y(i,j)
        enddo
      enddo

    else if (beta.eq.-ione) then
      do j=1,n
        do i=1,m
          y(i,j) = alpha*x(i,j) - y(i,j)
        enddo
      enddo
    else
      do j=1,n
        do i=1,m
          y(i,j) = alpha*x(i,j) + beta*y(i,j)
        enddo
      enddo
    endif

  endif

  return

9999 continue
  call fcpsb_serror()
  return

end subroutine iaxpby


text:
varying vec4 diffuse,ambient;
varying vec3 normal,lightDir,halfVector;

#define hasSpecMap %hasSpecMap%
#define hasBumpMap %hasBumpMap%
#define hasAlphaMap %hasAlphaMap%
#define hasColorMap %hasColorMap%
#define hasLumiMap %hasLumiMap%
#define hasShadowMap0 %hasShadowMap0%
#define hasShadowMap1 %hasShadowMap1%
#define hasShadowMap2 %hasShadowMap2%
#define hasShadowMap3 %hasShadowMap3%
#define hasShadowMap4 %hasShadowMap4%
#define hasShadowMap5 %hasShadowMap5%
#define hasShadowMap6 %hasShadowMap6%
#define hasShadowMap7 %hasShadowMap7%
#define hasEnvSphereMap %hasEnvSphereMap%
#define hasReflectionMap %hasReflectionMap%
#define hasLightMap %hasLightMap%
#define hasParallaxMap %hasParallaxMap%

#if hasSpecMap
uniform sampler2D specMap;
#endif

#if hasBumpMap
uniform sampler2D bumpMap;
#endif

//#if hasAlphaMap
//uniform sampler2D alphaMap;
//#endif

#if hasColorMap
uniform sampler2D colorMap;
#endif

#if hasLumiMap
uniform sampler2D lumiMap;
#endif

#if hasLightMap
uniform sampler2D lightMap;
#endif

#if hasEnvSphereMap
	uniform sampler2D envSphereMap;
	uniform vec4 projEnv;
	uniform float envAmount;
#if hasBumpMap
	varying vec3 u;
#endif
#endif

#if hasReflectionMap
	uniform sampler2D reflectionMap;
#endif

#if hasShadowMap0|hasShadowMap1|hasShadowMap2|hasShadowMap3|hasShadowMap4|hasShadowMap5|hasShadowMap6|hasShadowMap7
	uniform sampler2D shadowMap0,shadowMap1,shadowMap2;
	varying vec4 ProjShadow0;
	varying vec4 ProjShadow1;
	varying vec4 ProjShadow2;
#endif

#if hasParallaxMap
	uniform sampler2D parallaxMap;
	uniform float parallaxScale;
	uniform float parallaxBias;
#endif
#if hasParallaxMap||hasEnvSphereMap
	varying vec3 eyeVec;
	varying mat3 TBNMatrix;
#endif

float shadowVal(vec4 shadowProj, sampler2D shadowMap)
{
		vec4 shadowCoordinateWdivide;
		float distanceFromLight;

		shadowCoordinateWdivide = shadowProj / shadowProj.w;
		shadowCoordinateWdivide.w += 0.0005;

		distanceFromLight = texture2D(shadowMap,shadowCoordinateWdivide.st).z;

		return (distanceFromLight < shadowCoordinateWdivide.z && shadowProj.w > 0.00) ? 0.0 : 1.0;
}

void main()
{
#if hasShadowMap0|hasShadowMap1|hasShadowMap2|hasShadowMap3|hasShadowMap4|hasShadowMap5|hasShadowMap6|hasShadowMap7
		float shadowv = shadowVal(ProjShadow0, shadowMap0)*shadowVal(ProjShadow1, shadowMap1)*shadowVal(ProjShadow2, shadowMap2);
#endif

	vec3 n,halfV;
	float NdotL,NdotHV;

	vec4 color;
	vec4 specVal;

#if hasParallaxMap
	  float height = texture2D(parallaxMap, gl_TexCoord[0].xy).r;
	  float v = (height) * parallaxScale - parallaxBias;
	  vec3 eye = normalize(eyeVec);
	  vec2 texCoord = gl_TexCoord[0].xy + (eye.xy * v);
#else
	  vec2 texCoord = gl_TexCoord[0].xy;
#endif

#if hasBumpMap
			vec3 bumpNorm = vec3(texture2D(bumpMap, texCoord));
//			n = normalize(normal);
			n = normalize(normal+((bumpNorm-0.5)*2.0));
#else
			n = normalize(normal);
#endif

		/* compute the dot product between normal and ldir */

		NdotL = max(dot(n,lightDir),0.0);

		if (NdotL > 0.0)
		{
#if hasSpecMap
					specVal = texture2D(specMap, texCoord);
#else
					specVal = gl_FrontMaterial.specular;
#endif

		}

		color = diffuse;

#if hasColorMap
		color *= texture2D(colorMap, texCoord);
#endif

		vec4 litColor = vec4(0.0,0.0,0.0,1.0);

		litColor += NdotL;
//		halfV = normalize(halfVector);
//		NdotHV = max(dot(n,halfV),0.0);
//		color += specVal * gl_LightSource[%LIGHTID%].specular * pow(NdotHV, gl_FrontMaterial.shininess);

		NdotHV = max(dot(n, halfVector.xyz),0.001);

		litColor += specVal * gl_LightSource[%LIGHTID%].specular * pow(NdotHV,gl_FrontMaterial.shininess);

#if hasShadowMap0|hasShadowMap1|hasShadowMap2|hasShadowMap3|hasShadowMap4|hasShadowMap5|hasShadowMap6|hasShadowMap7
		litColor *= shadowv;
#endif

#if hasLightMap
		litColor.rgb += (texture2D(lightMap, gl_TexCoord[2].st).rgb*2.0)*texture2D(colorMap, texCoord).rgb;
#endif

#if hasLumiMap
		litColor.rgb += (vec3(1.0,1.0,1.0)-color.rgb) * 2.0 * texture2D(lumiMap, texCoord).rgb;
#endif

	litColor += ambient;

#if hasEnvSphereMap

#if hasBumpMap
	vec3 r = reflect( u, n );
	float m = 2.0 * sqrt( r.x*r.x + r.y*r.y + (r.z+1.0)*(r.z+1.0) );
//
	vec3 coord;
	coord.s = r.x/m + 0.5;
	coord.t = r.y/m + 0.5;

	#if hasReflectionMap
		color += litColor * texture2D( envSphereMap, coord.st) * texture2D( reflectionMap, texCoord);
	#else
		color += litColor * texture2D( envSphereMap, coord.st) * envAmount;
	#endif

#else
	#if hasReflectionMap
		color += litColor * texture2D( envSphereMap, gl_TexCoord[1].st) * texture2D( reflectionMap, texCoord);
	#else
		color += litColor * texture2D( envSphereMap, gl_TexCoord[1].st) * envAmount;
	#endif
#endif

#endif

		color *= litColor;

#if hasAlphaMap
		color.a = texture2D(colorMap, texCoord).a;
#else
		color.a = gl_FrontLightModelProduct.sceneColor.a;
#endif

	gl_FragColor = color;
}

ext: f

max_stars_repo_name: cjcliffe/CubicVR

max_issues_repo_name: mcanthony/CubicVR

max_forks_repo_name: mcanthony/CubicVR


text:
C$Procedure ZZDDHGTU ( Private --- DDH Get Unit )

      SUBROUTINE ZZDDHGTU ( UTCST, UTHAN, UTLCK, UTLUN, NUT, UINDEX )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines.  Users should not call this routine directly due
C     to the volatile nature of this routine.
C
C     Get or prepare an entry in the unit table to receive a new
C     file.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     PRIVATE
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE              'zzddhman.inc'

      INTEGER               UTCST  ( * )
      INTEGER               UTHAN  ( * )
      LOGICAL               UTLCK  ( * )
      INTEGER               UTLUN  ( * )
      INTEGER               NUT

      INTEGER               UINDEX

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     UTCST,
C     UTHAN,
C     UTLCK,
C     UTLUN,    I/O   Unit table.
C     NUT       I/O   Number of entries in the unit table.
C     UINDEX     O    Row in the unit table that can be replaced.
C
C$ Detailed_Input
C
C     UTCST,
C     UTHAN,
C     UTLCK,
C     UTLUN,     are the arrays respectively containing the cost,
C                handle, locked, and logical unit columns of the
C                unit table.
C
C     NUT        is the number of entries in the unit table.
C
C$ Detailed_Output
C
C     UTCST,
C     UTHAN,
C     UTLCK,
C     UTLUN,     are the arrays respectively containing the cost,
C                handle, locked, and logical unit columns of the
C                unit table.  This may change as a new unit is
C                added or old ones are removed.
C
C     NUT        is the number of entries in the unit table.  This may
C                change as new entries are added.
C
C     UINDEX     is the index of the row where the new unit should
C                be attached.
C
C$ Parameters
C
C     None.
C
C$ Files
C
C     This routine may disconnect a file from its logical unit, to
C     successfully process the caller's request for a unit.
C
C$ Exceptions
C
C     1) If GETLUN fails to assign a logical unit for any reason to
C        the row of interest, this routine sets the logical unit to -1,
C        since negative logical units in Fortran are not permitted.
C
C$ Particulars
C
C      This routine only manipulates the contents of the unit table.
C      Any "zero" cost rows in the table indicate rows where the
C      listed logical unit has been reserved, but no file is currently
C      attached.
C
C      Callers of this routine should check FAILED since this
C      routine may invoke GETLUN.
C
C$ Examples
C
C     See ZZDDHHLU for sample usage.
C
C$ Restrictions
C
C     1) This routine must not be used to retrieve a unit for a
C        file that is already connected to a unit listed in the
C        unit table.
C
C$ Author_and_Institution
C
C     F.S. Turner     (JPL)
C
C$ Literature_References
C
C     None.
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 29-MAY-2001 (FST)
C
C
C-&

C
C     SPICELIB Functions
C
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local Variables
C
      INTEGER               I
      INTEGER               ORDERV ( UTSIZE )

      LOGICAL               DONE

C
C     Standard SPICE discovery error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      END IF

C
C     First check the case when the unit table is completely empty.
C
      IF ( NUT .EQ. 0 ) THEN

         NUT    = 1
         UINDEX = 1

         UTCST(UINDEX) = 0
         UTHAN(UINDEX) = 0
         UTLCK(UINDEX) = .FALSE.

         CALL GETLUN ( UTLUN(UINDEX) )

C
C        Check FAILED to see if GETLUN signaled an error.  If so, then
C        return an invalid unit to the caller.
C
         IF ( FAILED() ) THEN
            UTLUN(UINDEX) = -1
            RETURN
         END IF

C
C        If we end up here, then GETLUN succeeded and we have the new
C        unit.  Now return.
C
         RETURN

      END IF

C
C     If we reach here, then the table contains at least one entry.
C     Order the table rows by cost.
C
      CALL ORDERI ( UTCST, NUT, ORDERV )

C
C     Now check to for '0' cost rows as this indicates rows whose
C     logical units are reserved for this suite of routines usage,
C     but are not currently assigned a file.
C
      IF ( UTCST(ORDERV(1)) .LE. 0 ) THEN

         UINDEX = ORDERV(1)

C
C        '0' cost rows end up in the unit table as the result of a
C        row deletion, occurring when excess files are present.
C        When this process occurs, the logical unit listed in this
C        row is reserved for this module's usage only with RESLUN.
C        Free it, since we're about to reassign it.
C
         CALL FRELUN ( UTLUN(UINDEX) )

         RETURN

      END IF

C
C     Now if no '0' cost rows exist, check to see if we can
C     expand the table.
C
      IF ( NUT .LT. UTSIZE ) THEN

C
C        Now increment NUT and set UINDEX.
C
         NUT    = NUT + 1
         UINDEX = NUT

C
C        Prepare the default values for the new row.
C
         UTCST(UINDEX) = 0
         UTHAN(UINDEX) = 0
         UTLCK(UINDEX) = .FALSE.

         CALL GETLUN ( UTLUN(UINDEX) )

C
C        Check FAILED to see if GETLUN signaled an error.  If so, then
C        return an invalid unit to the caller.
C
         IF ( FAILED() ) THEN
            UTLUN(UINDEX) = -1
            RETURN
         END IF

C
C        If we end up here, then GETLUN worked properly.  Now return.
C
         RETURN

      END IF

C
C     If we reach here, then we have no zero-cost rows and a full unit
C     table.  Now it's time to determine which entry in the table to
C     bump.  We do this by stepping through the order vector until
C     we find the first 'non-locked' row.
C
      I    = 0
      DONE = .FALSE.

      DO WHILE ( (.NOT. DONE) .AND. (I .NE. NUT) )

         I    = I + 1
         DONE = .NOT. UTLCK( ORDERV(I) )

      END DO

C
C     Before going any further, signal an error if we discover
C     we have not found a row.
C
      IF ( .NOT. DONE ) THEN

         UINDEX = 0
         CALL CHKIN  ( 'ZZDDHGTU' )
         CALL SETMSG ( 'The unit table is full and all entries are '
     .   //            'locked.  This should never happen. Contact '
     .   //            'NAIF.'                                       )
         CALL SIGERR ( 'SPICE(BUG)'                                  )
         CALL CHKOUT ( 'ZZDDHGTU'                                    )
         RETURN

      END IF

C
C     Clear UTCST and UTHAN since we intend to disconnect
C     the unit upon return.
C
      UTCST( ORDERV(I) ) = 0
      UTHAN( ORDERV(I) ) = 0

C
C     Set UINDEX and CLSLUN, then return.
C
      UINDEX = ORDERV(I)

C
C     At this point we need to close the unit from the row of interest.
C
      CLOSE ( UNIT = UTLUN(UINDEX) )

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
c     Sorting subroutine
c     Feb12/2018: arr must be DOUBLE PRECISION
c     23Feb2020: taken for usaf update files

      SUBROUTINE sort_doublelist_pvwave2(n,arr,brr)

      IMPLICIT NONE
c****
      INTEGER          n
      DOUBLE PRECISION arr(n)
      DOUBLE PRECISION crr(n)
      INTEGER          brr(n)
      INTEGER          i,j
      DOUBLE PRECISION a
      INTEGER          b
c****
c     Copy original array
      DO i=1,n
       brr(i)=i         !simple ascending index
       crr(i)=arr(i)
      ENDDO
c****
      DO j=2,n
       a=crr(j)
       b=brr(j)
       DO i=j-1,1,-1
        IF (crr(i).LE.a) GOTO 10
        crr(i+1)=crr(i)
        brr(i+1)=brr(i)
       ENDDO
       i=0
10     crr(i+1)=a
       brr(i+1)=b
      ENDDO

      RETURN
      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text:
C------------------------------------------------------------------------------
C
C                          NEKTON 2.6  2/8/90
C
C			Copyright (C) 1990, by the
C
C		Massachusetts Institute of Technology  and Nektonics, Inc.
C
C All Rights Reserved
C
C This program is a licenced product of MIT and Nektonics, Inc.,  and it is
C not to be disclosed to others, copied, distributed, or displayed
C without prior authorization.
C
C------------------------------------------------------------------------------
c-----------------------------------------------------------------------
      SUBROUTINE SORT(A,IND,N)
C
C     Use Heap Sort (p 233 Num. Rec.), 5/26/93 pff.
C
      DIMENSION A(1),IND(1)
C
      if (n.le.1) return
      DO 10 J=1,N
         IND(j)=j
   10 continue
C
      L=n/2+1
      ir=n
  100 CONTINUE
         IF (l.gt.1) THEN
            l=l-1
            indx=ind(l)
            q=a(indx)
         ELSE
            indx=ind(ir)
            q=a(indx)
            ind(ir)=ind(1)
            ir=ir-1
            if (ir.eq.1) then
               ind(1)=indx
               return
            endif
         ENDIF
         i=l
         j=l+l
  200    CONTINUE
         IF (J.le.IR) THEN
            IF (J.lt.IR) THEN
               IF ( A(IND(j)).lt.A(IND(j+1)) ) j=j+1
            ENDIF
            IF (q.lt.A(IND(j))) THEN
               IND(I)=IND(J)
               I=J
               J=J+J
            ELSE
               J=IR+1
            ENDIF
         GOTO 200
         ENDIF
         IND(I)=INDX
      GOTO 100
      END
c-----------------------------------------------------------------------
      SUBROUTINE SWAP(A,W,IND,N)
C
C     Use IND to sort array A   (p 233 Num. Rec.), 5/26/93 pff.
C
      DIMENSION A(1),W(1),IND(1)
C
      if (n.le.1) return
      DO 10 J=1,N
         W(j)=A(j)
   10 continue
C
      DO 20 J=1,N
         A(j)=W(ind(j))
   20 continue
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE ISWAP(A,W,IND,N)
C
C     Use IND to sort array A
C
      INTEGER A(1),W(1),IND(1)
C
      if (n.le.1) return
      DO 10 J=1,N
         W(j)=A(j)
   10 continue
C
      DO 20 J=1,N
         A(j)=W(ind(j))
   20 continue
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CHSWAP(A,W,L,IND,N)
C
C     Use IND to sort array A
C
      CHARACTER*1 A(L,1),W(L,1)
      INTEGER IND(1)
C
      if (n.le.1) return
      DO 10 J=1,L*N
         W(j,1)=A(j,1)
   10 continue
C
      DO 20 J=1,N
         ij=IND(J)
         DO 20 k=1,l
            A(l,j)=W(l,ij)
   20    continue
   30 continue
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE LJUST(STRING)
C     left justify string
      CHARACTER*1 STRING(80)
C
      IF (STRING(1).NE.' ') RETURN
C
      DO 100 I=2,80
C
         IF (STRING(I).NE.' ') THEN
            DO 20 J=1,81-I
               IJ=I+J-1
               STRING(J)=STRING(IJ)
   20       CONTINUE
            DO 30 J=82-I,80
               STRING(J)=' '
   30       CONTINUE
            RETURN
         ENDIF
C
  100 CONTINUE
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CAPIT(LETTRS,N)
C     Capitalizes string of length n
      CHARACTER LETTRS(N)
C
      DO 5 I=1,N
         INT=ICHAR(LETTRS(I))
         IF(INT.GE.97 .AND. INT.LE.122) THEN
            INT=INT-32
            LETTRS(I)=CHAR(INT)
         ENDIF
5     CONTINUE
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE PARSE(LINES,ARGS,IA)
C       Capitalizes LINE and splits it into 5 ARGuments of 10 Characters.
C
      CHARACTER ARGS(10,5),COMAND*10,LINES(70)
      LOGICAL GAP
      INTEGER IAS(80),IES(80)
C
      IA=0
      IE=0
        GAP=.TRUE.
      DO 3 I=1,10
         DO 3 J=1,5
            ARGS(I,J)= ' '
3     CONTINUE
C       Capitalize LINE
      CALL CAPIT(LINES,70)
C
      DO 10 I=1,70
         IF(LINES(I).NE.' ' .AND. LINES(I).NE.',') THEN
            IF(GAP) THEN
C! Found the start of new argument
               IA=IA+1
               IAS(IA)=I
            ELSE
            ENDIF
            GAP = .FALSE.
         ELSE
            IF(.NOT.(GAP)) THEN
C!Found End of argument
               IE=IE+1
               IES(IE)=I-1
            ELSE
            ENDIF
            GAP = .TRUE.
         ENDIF
10    CONTINUE
C       COMMAND
C     Max 5 Arguments
      IF(IA.GT.5)IA=5
      DO 50 I=1,IA
         DO 50 IC=IAS(I),IES(I)
            ARGS(IC-IAS(I)+1,I)= LINES(IC)
50    CONTINUE
      DO 60 I=1,10
60            COMAND(I:I)=ARGS(I,1)
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE BLANK(STRING,N)
      CHARACTER*1 STRING(N)
      CHARACTER*1   BLNK
      DATA BLNK/' '/
C
      DO 100 I=1,N
         STRING(I)=BLNK
  100 CONTINUE
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CADD2 (A,B,C,N)
      DIMENSION A(1),B(1)
      DO 100 I = 1, N
 100     A(I) = B(I) + C
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CADD (A,B,N)
      DIMENSION A(1)
      DO 100 I = 1, N
 100     A(I) = A(I) + B
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CMULT(A,B,N)
      DIMENSION A(1)
      DO 100 I = 1, N
 100     A(I) = B*A(I)
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CHCOPY(A,B,N)
      CHARACTER*1 A(1), B(1)
      DO 100 I = 1, N
 100     A(I) = B(I)
      RETURN
      END
C
c-----------------------------------------------------------------------
      FUNCTION LTRUNC(STRING,L)
      CHARACTER*1 STRING(L)
      CHARACTER*1   BLNK
      DATA BLNK/' '/
      DO 100 I=L,1,-1
         L1=I
         IF (STRING(I).NE.BLNK) GOTO 200
  100 CONTINUE
      L1=0
  200 CONTINUE
      LTRUNC=L1
      if (l1.eq.l) then
         write(6,*) 'LTRUNC: string:',l
         write(6,80) (string(k),k=1,l)
   80    format(80a1)
      endif
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE ROOTS(XVAL,NROOTS,FTARGT,F,N)
C
C     Find all values of I which yield F(I) close to FTARGT
C     and give the additional increment DI such that F(I+DI)+FTARGT,
C     where F(I+DI) is interpreted as the linear interpolation of F.
C
      DIMENSION XVAL(2,N),F(0:N)
C
      NROOTS=0
      DO 100 I=1,N
         F1=F(I-1)
         F2=F(I)
         DF1=FTARGT-F1
         DF2=F2-FTARGT
         IF (DF1*DF2.GE.0.0) THEN
            NROOTS=NROOTS+1
            DF=F2-F1
            IF (DF.EQ.0) THEN
               DI=0.5
            ELSE
               DI=(FTARGT-F1)/DF
            ENDIF
            XVAL(1,NROOTS)=FLOAT(I-1)
            XVAL(2,NROOTS)=DI
         ENDIF
  100 CONTINUE
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE VSQRT(A,N)
      DIMENSION  A(1)
      do i=1,n
         if (a(i).gt.0) a(i) = sqrt(a(i))
      enddo
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE VCNVERT(A,N)
      REAL*4 A(1)
      DO 100 I=1,N
         CALL CONVERT(A(I))
  100 CONTINUE
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE CONVERT(T)
C
      CHARACTER*12 VC
      CHARACTER*1  V1(12)
      EQUIVALENCE (V1,VC)
      CHARACTER*1  CSIGN
C
      REAL*4 W
      CHARACTER*4  WC
      CHARACTER*1  W1(4)
      EQUIVALENCE (W1,W)
      EQUIVALENCE (W1,WC)
C
      CHARACTER*1 ALPH64(0:63)
      SAVE        ALPH64
      DATA        ALPH64
     $  /'1','2','3','4','5','6','7','8','9','0'
     $   ,'a','b','c','d','e','f','g','h','i','j'
     $   ,'k','l','m','n','o','p','q','r','s','t'
     $   ,'u','v','w','x','y','z'
     $   ,'A','B','C','D','E','F','G','H','I','J'
     $   ,'K','L','M','N','O','P','Q','R','S','T'
     $   ,'U','V','W','X','Y','Z','+','-'/
C
C     Find out the usual decimal format for T
C
      WRITE(VC,10) T
   10 FORMAT(E12.5)
C
C     Begin converting the mantissa to base 64
C
      READ(VC,11) MANTIS
   11 FORMAT(3X,I5)
C
C     Sign?
C
      READ(VC,12) CSIGN
   12 FORMAT(A1)
      IF (CSIGN.EQ.'-') MANTIS=-MANTIS
      MANTIS=MANTIS+131072
C
C     ONES,TENS, HUNDREDS
C
      IONE=MOD(MANTIS,64)
      ITMP=MANTIS/64
      ITEN=MOD(ITMP,64)
      ITMP=ITMP/64
      IHUN=MOD(ITMP,64)
C
C     Exponent
C
      READ(VC,21) IEXP
   21 FORMAT(9X,I3)
C     We assume that the exponent is bounded by 31.
      IEXP=IEXP+31
C
C     Compute alpha equivalent
C
      W1(1)=ALPH64(IHUN)
      W1(2)=ALPH64(ITEN)
      W1(3)=ALPH64(IONE)
      W1(4)=ALPH64(IEXP)
C
C     Convert the input value
      T=W
      RETURN
      END
c-----------------------------------------------------------------------
      SUBROUTINE DELTMP
#     include "basics.inc"
      character*80 command
C     Remove all tmp.* files
C
      CALL BLANK(COMMAND,80)
      IF(IFVMS)THEN
C        VMS
         IERR=0
      ELSE
C        Ultrix
         write(command,10)
   10    format('rm tmp.*')
         CALL SYSTEM(command)
      ENDIF
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION ILMIN(A,N)
      integer A(N),temp
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = min(TEMP,A(I))
   10 CONTINUE
      ILmin=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION ILMAX(A,N)
      integer A(N),temp
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MAX(TEMP,A(I))
   10 CONTINUE
      ILMAX=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION VLMAX(A,N)
      DIMENSION A(N)
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MAX(TEMP,A(I))
   10 CONTINUE
      VLMAX=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION VLAMAX(A,N)
      DIMENSION A(N)
      TEMP = abs(A(1))
      DO 10 I=1,N
         TEMP = MAX(TEMP,abs(A(I)))
   10 CONTINUE
      VLAMAX=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION IGLSUM(A,N)
      integer A(N)
      ITEMP = 0
      DO 10 I=1,N
         ITEMP = ITEMP+A(I)
   10 CONTINUE
      IGLSUM=ITEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION IGLMIN(A,N)
      INTEGER A(N),TEMP
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MIN(TEMP,A(I))
   10 CONTINUE
      IGLMIN=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION IGLMAX(A,N)
      INTEGER A(N),TEMP
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MAX(TEMP,A(I))
   10 CONTINUE
      IGLMAX=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION IVLMAX(A,N)
      INTEGER A(N),TEMP
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MAX(TEMP,A(I))
   10 CONTINUE
      IVLMAX=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION GLMAX(A,N)
      DIMENSION A(N)
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MAX(TEMP,A(I))
   10 CONTINUE
      GLMAX=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION VLMIN(A,N)
      DIMENSION A(N)
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MIN(TEMP,A(I))
   10 CONTINUE
      VLMIN=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION GLMIN(A,N)
      DIMENSION A(N)
      TEMP = A(1)
      DO 10 I=1,N
         TEMP = MIN(TEMP,A(I))
   10 CONTINUE
      GLMIN=TEMP
      RETURN
      END
c-----------------------------------------------------------------------
      FUNCTION MOD1(I,N)
C
C     Yields MOD(I,N) with the exception that if I=K*N, result is N.
C
      MOD1 = 1
      IF (N.EQ.0) THEN
         CALL PRS(
     $  'WARNING:  Attempt to take MOD(I,0) in FUNCTION MOD1.$')
         RETURN
      ENDIF
      II = I+N-1
      MOD1 = MOD(II,N)+1
      RETURN
      END
c-----------------------------------------------------------------------
      INTEGER FUNCTION INDX1(S1,S2,L2)
      CHARACTER*80 S1,S2
C
      N1=80-L2+1
      INDX1=0
      IF (N1.LT.1) RETURN
C
      DO 100 I=1,N1
         I2=I+L2-1
         IF (S1(I:I2).EQ.S2(1:L2)) THEN
            INDX1=I
            RETURN
         ENDIF
  100 CONTINUE
C
      RETURN
      END
c-----------------------------------------------------------------------
      INTEGER FUNCTION NINDX1(S1,S2,L2)
C
C     Return index of first character Not equal to S2
C
      CHARACTER*80 S1,S2
C
      N1=80-L2+1
      NINDX1=0
      IF (N1.LT.1) RETURN
C
      DO 100 I=1,N1
         I2=I+L2-1
         IF (S1(I:I2).NE.S2(1:L2)) THEN
            NINDX1=I
            RETURN
         ENDIF
  100 CONTINUE
C
      RETURN
      END
c-----------------------------------------------------------------------
      subroutine isort(a,ind,n)
C
C     Use Heap Sort (p 231 Num. Rec., 1st Ed.)
C
      integer ind(1)
      integer    a(1)
      integer    aa
C
      dO 10 j=1,n
         ind(j)=j
   10 continue
C
      if (n.le.1) return
      L=n/2+1
      ir=n
  100 continue
         if (l.gt.1) then
            l=l-1
            aa  = a  (l)
            ii  = ind(l)
         else
                 aa =   a(ir)
                 ii = ind(ir)
              a(ir) =   a( 1)
            ind(ir) = ind( 1)
            ir=ir-1
            if (ir.eq.1) then
                 a(1) = aa
               ind(1) = ii
               return
            endif
         endif
         i=l
         j=l+l
  200    continue
         if (j.le.ir) then
            if (j.lt.ir) then
               if ( a(j).lt.a(j+1) ) j=j+1
            endif
            if (aa.lt.a(j)) then
                 a(i) = a(j)
               ind(i) = ind(j)
               i=j
               j=j+j
            else
               j=ir+1
            endif
         GOTO 200
         endif
           a(i) = aa
         ind(i) = ii
      GOTO 100
      end
c-----------------------------------------------------------------------
      subroutine sortit(a,ind,n)
C
C     Use Heap Sort (p 231 Num. Rec., 1st Ed.)
C
      integer ind(1)
      real    a(1)
      real    aa
C
      dO 10 j=1,n
         ind(j)=j
   10 continue
C
      if (n.le.1) return
      L=n/2+1
      ir=n
  100 continue
         if (l.gt.1) then
            l=l-1
            aa  = a  (l)
            ii  = ind(l)
         else
                 aa =   a(ir)
                 ii = ind(ir)
              a(ir) =   a( 1)
            ind(ir) = ind( 1)
            ir=ir-1
            if (ir.eq.1) then
                 a(1) = aa
               ind(1) = ii
               return
            endif
         endif
         i=l
         j=l+l
  200    continue
         if (j.le.ir) then
            if (j.lt.ir) then
               if ( a(j).lt.a(j+1) ) j=j+1
            endif
            if (aa.lt.a(j)) then
                 a(i) = a(j)
               ind(i) = ind(j)
               i=j
               j=j+j
            else
               j=ir+1
            endif
         GOTO 200
         endif
           a(i) = aa
         ind(i) = ii
      GOTO 100
      end
c-----------------------------------------------------------------------
      subroutine add3s12(x,y,cy,z,cz,n)
      real x(1),y(1),z(1),cy,cz
c
      do i=1,n
         x(i) = cy*y(i) + cz*z(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine add3s2(x,y,z,c,n)
      real x(1),y(1),z(1),c
c
      do i=1,n
         x(i) = y(i) + c*z(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine cmult2(x,y,c,n)
      real x(1),y(1),c
c
      do i=1,n
         x(i) = c*y(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine add2s1(x,y,c,n)
      real x(1),y(1),c
c
      do i=1,n
         x(i) = c*x(i) + y(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine add2s2(x,y,c,n)
      real x(1),y(1),c
c
      do i=1,n
         x(i) = x(i) + c*y(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      function glsc2(x,y,n)
C
C     Perform inner-product
C
      real x(1), y(1)
c
      sum = 0.0
      do i=1,n
         sum=sum+x(i)*y(i)
      enddo
c
      glsc2 = sum
c
      return
      end
c-----------------------------------------------------------------------
      function vlsum(x,n)
      real x(1)
      sum = 0.0
      do i=1,n
         sum=sum+x(i)
      enddo
      vlsum = sum
      return
      end
c-----------------------------------------------------------------------
      function glsc1(x,n)
C     Perform inner-product
      real x(1)
c
      sum = 0.0
      do i=1,n
         sum=sum+x(i)*x(i)
      enddo
      glsc1 = sum
      return
      end
c-----------------------------------------------------------------------
      SUBROUTINE IRANK(A,IND,N)
C
C     Use Heap Sort (p 233 Num. Rec.)
C
      INTEGER A(1),IND(1)
      INTEGER Q
c
      if (n.le.1) return
      DO 10 J=1,N
         IND(j)=j
   10 continue
C
      L=n/2+1
      ir=n
  100 CONTINUE
         IF (l.gt.1) THEN
            l=l-1
            indx=ind(l)
            q=a(indx)
         ELSE
            indx=ind(ir)
            q=a(indx)
            ind(ir)=ind(1)
            ir=ir-1
            if (ir.eq.1) then
               ind(1)=indx
               return
            endif
         ENDIF
         i=l
         j=l+l
  200    CONTINUE
         IF (J.le.IR) THEN
            IF (J.lt.IR) THEN
               IF ( A(IND(j)).lt.A(IND(j+1)) ) j=j+1
            ENDIF
            IF (q.lt.A(IND(j))) THEN
               IND(I)=IND(J)
               I=J
               J=J+J
            ELSE
               J=IR+1
            ENDIF
         GOTO 200
         ENDIF
         IND(I)=INDX
      GOTO 100
      END
c-----------------------------------------------------------------------
      SUBROUTINE RANK(A,IND,N)
C
C     Use Heap Sort (p 233 Num. Rec.), 5/26/93 pff.
C
      DIMENSION A(1),IND(1)
C
      if (n.le.1) return
      DO 10 J=1,N
         IND(j)=j
   10 continue
C
      L=n/2+1
      ir=n
  100 CONTINUE
         IF (l.gt.1) THEN
            l=l-1
            indx=ind(l)
            q=a(indx)
         ELSE
            indx=ind(ir)
            q=a(indx)
            ind(ir)=ind(1)
            ir=ir-1
            if (ir.eq.1) then
               ind(1)=indx
               return
            endif
         ENDIF
         i=l
         j=l+l
  200    CONTINUE
         IF (J.le.IR) THEN
            IF (J.lt.IR) THEN
               IF ( A(IND(j)).lt.A(IND(j+1)) ) j=j+1
            ENDIF
            IF (q.lt.A(IND(j))) THEN
               IND(I)=IND(J)
               I=J
               J=J+J
            ELSE
               J=IR+1
            ENDIF
         GOTO 200
         ENDIF
         IND(I)=INDX
      GOTO 100
      END
c-----------------------------------------------------------------------
      SUBROUTINE COL3(A,B,C,N)
      REAL A(N),B(N),C(N)
      DO 100 I=1,N
         A(I)=B(I)*C(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE SUBCOL3(A,B,C,N)
      REAL A(N),B(N),C(N)
      DO 100 I=1,N
         A(I)=A(I)-B(I)*C(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE ADDCOL3(A,B,C,N)
      REAL A(N),B(N),C(N)
      DO 100 I=1,N
         A(I)=A(I)+B(I)*C(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE ADDCOL4(A,B,C,D,N)
      REAL A(N),B(N),C(N),D(N)
      DO 100 I=1,N
         A(I)=A(I)+B(I)*C(I)*D(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE SUBCOL4(A,B,C,D,N)
      REAL A(N),B(N),C(N),D(N)
      DO 100 I=1,N
         A(I)=A(I)-B(I)*C(I)*D(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE IZERO (IA,N)
      DIMENSION  IA(1)
      DO 100 I = 1, N
 100     IA(I) = 0
      RETURN
      END
      SUBROUTINE RZERO (A,N)
      DIMENSION  A(1)
      DO 100 I = 1, N
 100     A(I ) = 0.0
      RETURN
      END
      SUBROUTINE VDOT(DOT,U1,U2,U3,V1,V2,V3,N)
C
C     Compute an elemental based Cartesian vector dot product.
C
      DIMENSION DOT(N)
      DIMENSION U1(N),U2(N),U3(N)
      DIMENSION V1(N),V2(N),V3(N)
C
      DO 100 I=1,N
         DOT(I) = U1(I)*V1(I) + U2(I)*V2(I) + U3(I)*V3(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE INVCOL2 (A,B,N)
      REAL A(1),B(1)
      DO 100 I=1,N
         A(I)=A(I)/B(I)
 100  CONTINUE
      RETURN
      END
      SUBROUTINE COPY(A,B,N)
      REAL A(N),B(N)
      DO 10 I=1,N
         A(I)=B(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE ICOPY(A,B,N)
      INTEGER A(1),B(1)
      DO 10 I=1,N
         A(I)=B(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE SUB2(A,B,N)
      REAL A(1),B(1)
      DO 10 I=1,N
         A(I)=A(I)-B(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE SUB3(C,A,B,N)
      REAL C(1),A(1),B(1)
      DO 10 I=1,N
         C(I)=A(I)-B(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE ADD3(C,A,B,N)
      REAL C(1),A(1),B(1)
      DO 10 I=1,N
         C(I)=A(I)+B(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE ADD2(A,B,N)
      REAL A(N),B(N)
      DO 1 I=1,N
         A(I)=A(I)+B(I)
1     CONTINUE
      RETURN
      END
      FUNCTION ALMIN(A,N)
C
C     Compute the local minimum of quantity A
C
      REAL A(N)
C
C  Compute the local sum
C
      TMIN = A(1)
      DO 100 I=1,N
         TMIN=AMIN1(TMIN,A(I))
  100 CONTINUE
      ALMIN = TMIN
      RETURN
      END
      FUNCTION ALMAX(A,N)
C
C     Compute the local maximum of quantity A
C
      REAL A(N)
C
C  Compute the local sum
C
      TMAX = A(1)
      DO 100 I=1,N
         TMAX=AMAX1(TMAX,A(I))
  100 CONTINUE
      ALMAX = TMAX
      RETURN
      END
c-----------------------------------------------------------------------
      subroutine iint(in,n)
      integer in(1)
      do i=1,n
         in(i) = i
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine irank_vec(ind,nn,a,m,n,key,nkey,aa)
c
c     Compute rank of each unique entry a(1,i)
c
c     Output:   ind(i)  i=1,...,n    (global) rank of entry a(*,i)
c               nn  = max(rank)
c               a(j,i) is destroyed
c
c     Input:    a(j,i) j=1,...,m;  i=1,...,n
c               m      :   leading dim. of v  (ldv must be .ge. m)
c               key    :   sort key
c               nkey   :
c
c     Although not mandatory, this ranking procedure is probably
c     most effectively employed when the keys are pre-sorted. Thus,
c     the option is provided to sort vi() prior to the ranking.
c
c
      integer ind(n),a(m,n)
      integer key(nkey),aa(m)
      logical iftuple_ianeb,a_ne_b
c
      if (m.eq.1) then
c
         write(6,*)
     $        'WARNING: For single key, not clear that rank is unique!'
         call irank(a,ind,n)
         return
      endif
c
c
      nk = min(nkey,m)
      call ituple_sort(a,m,n,key,nk,ind,aa)
c
c     Find unique a's
c
      nn=1
c
      call icopy(aa,a,m)
      a(1,1) = nn
      a(2,1)=ind(1)
c
      do i=2,n
         a_ne_b = iftuple_ianeb(aa,a(1,i),key,nk)
c        write(6,1) nn,i,a_ne_b,aa,(a(k,i),k=1,nk)
c  1     format(2i7,1x,l4,1x,12i6)
         if (a_ne_b) then
            call icopy(aa,a(1,i),m)
            nn = nn+1
         endif
         a(1,i) = nn
         a(2,i) = ind(i)
      enddo
c
c     Set ind() to rank
c
      do i=1,n
         iold=a(2,i)
         ind(iold) = a(1,i)
      enddo
c
      return
      end
c-----------------------------------------------------------------------
      subroutine cfill(a,b,n)
      real a(1),b
      DO 100 I = 1, N
 100     A(I) = B
      RETURN
      END
c-----------------------------------------------------------------------
      subroutine ifill(a,b,n)
      integer a(1),b
      DO 100 I = 1, N
 100     A(I) = B
      RETURN
      END
c-----------------------------------------------------------------------
      subroutine semhat(a,b,c,d,dp,jp,z,n,w)
c
c     Generate matrices for single element, 1D operators:
c
c        a  = Laplacian
c        b  = diagonal mass matrix
c        c  = convection operator b*d
c        d  = derivative matrix
c        dp = derivative matrix, mapping from pressure nodes to velocity
c        jp = interpolation matrix, mapping from pressure nodes to velocity
c        z  = GLL points
c        n  = polynomial degree (velocity space)
c        w  = work array of size 2*n+2
c
c     Currently, this is set up for pressure nodes on the interior GLL pts.
c
c
      real a(0:n,0:n),b(0:n),c(0:n,0:n),d(0:n,0:n),z(0:n)
      real dp(0:n,1:n-1),jp(0:n,1:n-1)
      real w(0:1)
c
      np = n+1
      nm = n-1
      n2 = n-2
c
      call zwgll (z,b,np)
c
      do i=0,n
         call fd_weights_full(z(i),z,n,1,w)
         do j=0,n
            d(i,j) = w(j+np)                   !  Derivative matrix
         enddo
      enddo
c
      do i=0,n
         call fd_weights_full(z(i),z(1),n2,1,w(1))
         do j=1,nm
            jp(i,j) = w(j   )                  !  Interpolation matrix
            dp(i,j) = w(j+nm)                  !  Derivative    matrix
         enddo
      enddo
c
      call rzero(a,np*np)
      do j=0,n
      do i=0,n
         do k=0,n
            a(i,j) = a(i,j) + d(k,i)*b(k)*d(k,j)
         enddo
         c(i,j) = b(i)*d(i,j)
      enddo
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine generalev(a,b,lam,n,w)
c
c     Solve the generalized eigenvalue problem  A x = lam B x
c
c     A -- symm.
c     B -- symm., pos. definite
c
      integer wdsize
      save    wdsize
      data    wdsize/0/
c
      real a(n,n),b(n,n),lam(n),w(n,n)
c
c
c     First call: determine word size for dgemm/sgemm discrimination, below.
      if (wdsize.eq.0) then
         one = 1.0
         eps = 1.e-12
         wdsize = 8
         if (one+eps.eq.1.0) wdsize = 4
      endif
c     write(6,*) 'in generalev, =',n,' wdsize',wdsize
c
c
c     call outmat(a,n,n,'Agn')
c     call outmat(b,n,n,'Bgn')
c
      lwork = max(n*n,3*n-1)
      call ssygv(1,'V','U',n,a,n,b,n,lam,w,lwork,info)
c
c
      if (info.ne.0) then
         ninf = n-info
         write(6,*) 'Error in generalev, info=',info,n,ninf
c
         call outmat(a,n,n,'Agn')
         call outmat(b,n,n,'Bgn')
c
         call exitt
      endif
c
c     call outmat(a,n,n,'Agn')
c
      return
      end
c-----------------------------------------------------------------------

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
!==============================================================================!
  subroutine Surf_Mod_Count_Elements_Neighbours(surf)
!------------------------------------------------------------------------------!
!   Counts elements neighbours                                                 !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Surf_Type), target :: surf
!-----------------------------------[Locals]-----------------------------------!
  type(Elem_Type), pointer :: elem(:)
  integer,         pointer :: ne
  integer                  :: e
!==============================================================================!

  ! Take aliases
  ne   => surf % n_elems
  elem => surf % elem

  ! Do the actual count
  elem(1:ne) % nne = 0
  do e = 1, ne
    if(elem(e) % ei .ne. 0) elem(e) % nne = elem(e) % nne + 1
    if(elem(e) % ej .ne. 0) elem(e) % nne = elem(e) % nne + 1
    if(elem(e) % ek .ne. 0) elem(e) % nne = elem(e) % nne + 1
  end do

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows