text:
!----------------------------------------------------------------------------
! This file is part of UCLALES.
!
! UCLALES 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; either version 3 of the License, or
! (at your option) any later version.
!
! UCLALES 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/>.
!
! Copyright 1999-2007, Bjorn B. Stevens, Dep't Atmos and Ocean Sci, UCLA
!----------------------------------------------------------------------------
!
module advl

  implicit none

contains
  !
  ! ----------------------------------------------------------------------
  ! LADVECT:  This does the leap-frog advection using a fourth
  ! order centered in space algorithm.  Boundary point fluxes are not
  ! computed... hence the routine is hardwired to a cyclic domain.
  ! Vertical advection is density weighted consistent with the anelastic
  ! approximation
!
  subroutine ladvect

    use grid, only : a_uc, a_vc, a_wc, a_ut, a_vt, a_wt, a_scr1, a_scr2,      &
         nxp, nyp, nzp, dzt, dzm, dxi, dyi, dn0
    use stat, only : sflg, updtst, acc_tend
    use util, only : get_avg3

    real ::  v1da(nzp), v1db(nzp), v1dc(nzp), v1dd(nzp), v1de(nzp)
    !
    ! prepare density weights for use vertical advection
    !
    if (sflg) call acc_tend(nzp,nxp,nyp,a_uc,a_vc,a_wc,a_ut,a_vt,a_wt,        &
         v1dc,v1dd,v1de,1,'adv')

    call advl_prep(nzp,nxp,nyp,a_wc,a_scr1,dn0,dzt,dzm,v1da,v1db)
    !
    ! advection of u by (u,v,w) all at current timelevel.  also when flag
    ! is set updated statistical array with uw flux derived from ladvzu
    !
    call ladvxu(nzp,nxp,nyp,a_uc,a_ut,a_scr2,dxi)
    call ladvyu(nzp,nxp,nyp,a_uc,a_ut,a_vc,a_scr2,dyi)
    call ladvzu(nzp,nxp,nyp,a_uc,a_ut,a_scr1,a_scr2,v1da)
    if (sflg) then
       call get_avg3(nzp,nxp,nyp,a_scr2,v1da)
       call updtst(nzp,'adv',-1,v1da,1)
    end if
    !
    ! advection of v by (u,v,w) all at current timelevel.  also when flag
    ! is set updated statistical array with uw flux derived from ladvzu
    !
    call ladvxv(nzp,nxp,nyp,a_uc,a_vc,a_vt,a_scr2,dxi)
    call ladvyv(nzp,nxp,nyp,a_vc,a_vt,a_scr2,dyi)
    call ladvzv(nzp,nxp,nyp,a_vc,a_vt,a_scr1,a_scr2,v1da)
    if (sflg) then
       call get_avg3(nzp,nxp,nyp,a_scr2,v1da)
       call updtst(nzp,'adv',-2,v1da,1)
    end if
    !
    ! advection of w by (u,v,w) all at current timelevel.  also when flag
    ! is set updated statistical array with uw flux derived from ladvzu
    !
    call ladvxw(nzp,nxp,nyp,a_uc,a_wc,a_wt,a_scr2,dxi)
    call ladvyw(nzp,nxp,nyp,a_vc,a_wc,a_wt,a_scr2,dyi)
    call ladvzw(nzp,nxp,nyp,a_wc,a_wt,a_scr1,a_scr2,v1db)

    if (sflg) then
       call get_avg3(nzp,nxp,nyp,a_scr2,v1da)
       call updtst(nzp,'adv',-3,v1da,1)
    end if

    if (sflg) call acc_tend(nzp,nxp,nyp,a_uc,a_vc,a_wc,a_ut,a_vt,a_wt,        &
         v1dc,v1dd,v1de,2,'adv')

  end subroutine ladvect
  !
  ! ----------------------------------------------------------------------
  ! LADVXU: Advection of U, by U, div is a scratch array,
  ! tendencies are accumulated if ut variable, dxi is the inverse of
  ! delta-x the grid spacing in the horizontal direction.
  !
  subroutine ladvxu(n1,n2,n3,uc,ut,flx,dxi)

    integer, intent (in) ::  n1,n2,n3
    real, intent(in)     :: uc(n1,n2,n3),dxi
    real, intent(inout)  :: ut(n1,n2,n3)
    real, intent(out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=1,n3
       do i=3,n2-1
          do k=2,n1-1
             flx(k,i,j)=(0.583333*(uc(k,i,j)+uc(k,i-1,j)) - 0.083333          &
                  *(uc(k,i+1,j)+uc(k,i-2,j)))*.5*(uc(k,i,j)+uc(k,i-1,j))
          end do
       end do

       do i=3,n2-2
          do k=2,n1-1
             ut(k,i,j)=ut(k,i,j)-(flx(k,i+1,j)-flx(k,i,j))*dxi
          end do
       end do
    end do

  end subroutine ladvxu
  !
  ! ----------------------------------------------------------------------
  ! LADVYU: Advection of U, by V, div is a scratch array,
  ! tendencies are accumulated if fu variable, dyi is the inverse of
  ! delta-y.
  !
  subroutine ladvyu(n1,n2,n3,uc,ut,vc,flx,dyi)

    integer, intent (in)  ::  n1,n2,n3
    real, intent (in)     :: uc(n1,n2,n3),vc(n1,n2,n3),dyi
    real, intent (in out) :: ut(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=2,n3-2
       do i=1,n2-1
          do k=2,n1-1
             flx(k,i,j)=(0.583333*(uc(k,i,j)+uc(k,i,j+1)) - 0.083333          &
                  *(uc(k,i,j-1)+uc(k,i,j+2)))*.5*(vc(k,i,j)+vc(k,i+1,j))
          end do
       end do
    end do

    do j=3,n3-2
       do i=3,n2-2
          do k=2,n1-1
             ut(k,i,j)=ut(k,i,j)-(flx(k,i,j)-flx(k,i,j-1))*dyi
          end do
       end do
    end do

  end subroutine ladvyu
  !
  ! ----------------------------------------------------------------------
  ! Subroutine ladvzu: Advection of U, by W*rho, div is a scratch array,
  ! tendencies are accumulated if ut variable, v1 is the inverse density
  ! array valid at thermo points
  !
  subroutine ladvzu(n1,n2,n3,uc,ut,wm,flx,v1)

    integer, intent (in)  ::  n1,n2,n3
    real, intent (in)     :: uc(n1,n2,n3),wm(n1,n2,n3),v1(n1)
    real, intent (in out) :: ut(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=1,n3
       do i=1,n2-1
          flx(1,i,j)=0.
          do k=2,n1-2
             flx(k,i,j)=(0.583333*(uc(k,i,j)+uc(k+1,i,j)) - 0.083333          &
                  *(uc(k-1,i,j)+uc(k+2,i,j)))*.5*(wm(k,i,j)+wm(k,i+1,j))
          end do
          flx(n1-1,i,j)=0.
          flx(n1,i,j)=0.
       end do

       do i=1,n2-1
          do k=2,n1-1
             ut(k,i,j)=ut(k,i,j)-(flx(k,i,j)-flx(k-1,i,j))*v1(k)
          end do
       end do
    end do

  end subroutine ladvzu
  !
  ! ----------------------------------------------------------------------
  ! LADVXV: Advection of V, by U, div is a scratch array,
  ! tendencies are accumulated if vt variable, dxi is the inverse of
  ! delta-x the grid spacing in the horizontal direction.
  !
  subroutine ladvxv(n1,n2,n3,uc,vc,vt,flx,dxi)

    integer, intent (in)  ::  n1,n2,n3
    real, intent (in)     :: uc(n1,n2,n3),vc(n1,n2,n3),dxi
    real, intent (in out) :: vt(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=1,n3-1
       do i=2,n2-2
          do k=2,n1-1
             flx(k,i,j)=(0.583333*(vc(k,i,j)+vc(k,i+1,j)) - 0.083333          &
                  *(vc(k,i-1,j)+vc(k,i+2,j)))*.5*(uc(k,i,j)+uc(k,i,j+1))
          end do
       end do

       do i=3,n2-2
          do k=2,n1-1
             vt(k,i,j)=vt(k,i,j)-(flx(k,i,j)-flx(k,i-1,j))*dxi
          end do
       end do
    end do

  end subroutine ladvxv
  !
  ! ----------------------------------------------------------------------
  ! LADVYV: Advection of V, by V, div is a scratch array,
  ! tendencies are accumulated if vt variable, dyi is the inverse of
  ! delta-y.
  !
  subroutine ladvyv(n1,n2,n3,vc,vt,flx,dyi)

    integer, intent (in)  ::  n1,n2,n3
    real, intent (in)     :: vc(n1,n2,n3),dyi
    real, intent (in out) :: vt(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=3,n3-1
       do i=1,n2
          do k=2,n1-1
             flx(k,i,j)=(0.583333*(vc(k,i,j)+vc(k,i,j-1))-0.083333           &
                  *(vc(k,i,j+1)+vc(k,i,j-2)))*.5*(vc(k,i,j)+vc(k,i,j-1))
          end do
       end do
    end do

    do j=3,n3-2
       do i=1,n2
          do k=2,n1-1
             vt(k,i,j)=vt(k,i,j)-(flx(k,i,j+1)-flx(k,i,j))*dyi
          end do
       end do
    end do

  end subroutine ladvyv
  !
  ! ----------------------------------------------------------------------
  ! LADVZV: Advection of V, by W*rho, div is a scratch array,
  ! tendencies are accumulated if vt variable, v1 is the inverse density
  ! array valid at thermo points
  !
  subroutine ladvzv(n1,n2,n3,vc,vt,wm,flx,v1)

    integer, intent (in)  ::  n1,n2,n3
    real, intent (in)     :: vc(n1,n2,n3),wm(n1,n2,n3),v1(n1)
    real, intent (in out) :: vt(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=1,n3-1
       do i=1,n2
          flx(1,i,j)=0.
          do k=2,n1-2
             flx(k,i,j)=(0.583333*(vc(k,i,j)+vc(k+1,i,j))-0.083333            &
                  *(vc(k-1,i,j)+vc(k+2,i,j)))*.5*(wm(k,i,j)+wm(k,i,j+1))
          end do
          flx(n1-1,i,j)=0.
          flx(n1,i,j)=0.
       end do

       do i=1,n2
          do k=2,n1-1
             vt(k,i,j)=vt(k,i,j)-(flx(k,i,j)-flx(k-1,i,j))*v1(k)
          end do
       end do
    end do

  end subroutine ladvzv
  !
  ! ----------------------------------------------------------------------
  ! LADVXW: Advection of W, by U, div is a scratch array,
  ! tendencies are accumulated if wt variable, dxi is the inverse of
  ! delta-x the grid spacing in the horizontal direction.
  !
  subroutine ladvxw(n1,n2,n3,uc,wc,wt,flx,dxi)

    integer, intent (in)  :: n1,n2,n3
    real, intent (in)     :: uc(n1,n2,n3),wc(n1,n2,n3),dxi
    real, intent (in out) :: wt(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=1,n3
       do i=2,n2-2
          do k=2,n1-2
             flx(k,i,j)=(0.583333*(wc(k,i,j)+wc(k,i+1,j))-0.083333*          &
                  (wc(k,i-1,j)+wc(k,i+2,j)))*.5*(uc(k,i,j)+uc(k+1,i,j))
          end do
       end do

       do i=3,n2-2
          do k=2,n1-2
             wt(k,i,j)=wt(k,i,j)-(flx(k,i,j)-flx(k,i-1,j))*dxi
          end do
       end do
    end do

  end subroutine ladvxw
  !
  ! ----------------------------------------------------------------------
  ! LADVYW: Advection of W, by V, div is a scratch array,
  ! tendencies are accumulated if wt variable, dyi is the inverse of
  ! delta-y.
  !
  subroutine ladvyw(n1,n2,n3,vm,wc,wt,flx,dyi)

    integer, intent (in)  ::  n1,n2,n3
    real, intent (in)     :: vm(n1,n2,n3),wc(n1,n2,n3),dyi
    real, intent (in out) :: wt(n1,n2,n3)
    real, intent (out)    :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=2,n3-2
       do i=1,n2
          do k=2,n1-2
             flx(k,i,j)=(0.583333*(wc(k,i,j)+wc(k,i,j+1))-0.083333            &
                  *(wc(k,i,j-1)+wc(k,i,j+2)))*.5*(vm(k,i,j)+vm(k+1,i,j))
          end do
       end do
    end do

    do j=3,n3-2
       do i=1,n2
          do k=2,n1-2
             wt(k,i,j)=wt(k,i,j)-(flx(k,i,j)-flx(k,i,j-1))*dyi
          end do
       end do
    end do

  end subroutine ladvyw
  !
  ! ----------------------------------------------------------------------
  ! LADVZW: Advection of W, by W*rho, div is a scratch array,
  ! tendencies are accumulated if wt variable, v2 is the inverse density
  ! array valid at w-points
  !
  subroutine ladvzw(n1,n2,n3,wc,wt,wm,flx,v2)

    integer, intent (in) ::  n1,n2,n3
    real, intent (in)    :: wm(n1,n2,n3),wc(n1,n2,n3),v2(n1)
    real, intent (inout) :: wt(n1,n2,n3)
    real, intent (out)   :: flx(n1,n2,n3)

    integer :: i,j,k

    do j=1,n3
       do i=1,n2
          flx(1,i,j)=0.
          flx(2,i,j)=wm(2,i,j)*.5*(0.583333*(wc(2,i,j))-0.083333*(wc(3,i,j)))
          do k=3,n1-1
             flx(k,i,j)=(0.583333*(wc(k,i,j)+wc(k-1,i,j))-0.083333            &
                  *(wc(k+1,i,j)+wc(k-2,i,j)))*.5*(wm(k,i,j)+wm(k-1,i,j))
          end do
          flx(n1,i,j)=0.

          do k=2,n1-2
             wt(k,i,j)=wt(k,i,j)-(flx(k+1,i,j)-flx(k,i,j))*v2(k)
          end do
       end do
    end do

  end subroutine ladvzw
  !
  ! ----------------------------------------------------------------------
  ! ADVL_PREP: prepares two scratch arrays with the inverse
  ! densities as they locate on thermo levels (v1) and w-levels (v2)
  !
  subroutine advl_prep(n1,n2,n3,w,wm,dn0,dzt,dzm,v1,v2)

    integer, intent (in) :: n1,n2,n3
    real, intent (in)    ::  w(n1,n2,n3),dn0(n1),dzt(n1),dzm(n1)
    real, intent (out)   ::  wm(n1,n2,n3),v1(n1),v2(n1)

    integer :: k

    do k=1,n1-1
       wm(k,:,:)=w(k,:,:)*(dn0(k)+dn0(k+1))*.5
       v1(k)=dzt(k)/dn0(k)
       v2(k)=2.*dzm(k)/(dn0(k)+dn0(k+1))
    end do

  end subroutine advl_prep

end module advl

max_stars_repo_name: waqarnabi/tybec

max_issues_repo_name: waqarnabi/tybec

max_forks_repo_name: waqarnabi/tybec


text: The Perimeter Greenbelt is the section of the The Greenbelt Davis Greenbelt that connects to Northstar Park and the Senda Nueva Greenbelt.  It is located in North Davis, along the Covell Drainage Channel.  The North Davis Fitness Trail has various exercise stations along the trail.  Near the southwest end of the perimeter, there is an Greenbelt Secret Path unmarked path that connects to the south end of the Senda Nueva Greenbelt Senda Nueva greenbelt.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
subroutine sud(dischy    ,nst       ,icreep    ,betac     ,mmax      , &
             & nmax      ,j         ,nmmaxj    ,nmmax     ,kmax      , &
             & lstsci    ,nsrc      ,lsecfl    ,norow     ,icx       , &
             & icy       ,dismmt    ,irocol               ,mnksrc    , &
             & kfu       ,kfv       ,kfs       ,kcs       ,kspu      , &
             & kadu      ,kadv      ,kcu       ,kfumin    ,kfumax    , &
             & porosu    ,s0        ,s1        ,u0        ,u1        , &
             & v1        ,w1        ,r0        ,qxk       ,qyk       , &
             & qzk       ,guu       ,gvv       ,gvu       ,gsqs      , &
             & gud       ,gvd       ,gvz       ,gsqiu     ,dteu      , &
             & circ2d    ,circ3d    ,disch     ,                       &
             & umdis     ,umean     ,hu        ,hv        ,dpu       ,dzu1      , &
             & dpdksi    ,thick     ,sig       ,dps       ,taubpu    , &
             & taubsu    ,rho       ,sumrho    ,wsu       ,fxw       , &
             & wsbodyu   ,idry      ,crbc      ,vicuv     ,hu0       , &
             & vnu2d     ,vicww     ,rxx       ,rxy       ,dfu       , &
             & deltau    ,tp        ,rlabda    ,cfurou    ,cfvrou    , &
             & rttfu     ,diapl     ,rnpl      , &
             & windsu    ,patm      ,fcorio    ,evap      ,ubrlsu    , &
             & uwtypu    ,hkru      ,pship     ,tgfsep    ,a         , &
             & b         ,c         ,d         ,aa        ,bb        , &
             & cc        ,dd        ,tetau     ,aak       ,bbk       , &
             & cck       ,ddk       ,d0        ,d0k       ,bbka      , &
             & bbkc      ,ua        ,ub        ,soumud    ,dis_nf    , &
             & precip    ,ustokes   ,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: sud.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/compute/sud.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: SUD evaluates/solves the implicitly coupled
!              momentum and continuity equation at each
!              half time step.
!              Special approximation for pressure term,
!              based on limiter to avoid artificial flow.
!              Switch which makes it possible to use
!              upwind-approach for wet cross section in shallow
!              areas or if the model area contains structures.
! Method used: A.D.I.-scheme is used.
!              Upwind-approach for wet cross section in shallow
!              areas or if the model area contains structures.
!              At 2D Weir points:
!              - depth value DPU is not corrected like in general
!                (3D) weir case (see CALDPU)
!              - crest height is explicitly taken into account
!                in drying check
!              - 2D Turbulence model at depth points
!
! UA and UB are workarrays WRKB15 and WRKB16; Used only in CUCNP
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use flow2d3d_timers
    use globaldata
    use dfparall
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    include 'flow_steps_f.inc'
    real(fp)               , pointer :: eps
    integer                , pointer :: maseva
    integer                , pointer :: lundia
    integer                , pointer :: ntstep
    real(fp)               , pointer :: dco
    real(fp)               , pointer :: drycrt
    real(fp)               , pointer :: dryflc
    real(fp)               , pointer :: hdt
    integer                , pointer :: iter1
    character(6)           , pointer :: momsol
    character(8)           , pointer :: dpuopt
    real(fp)               , pointer :: rhow
    real(fp)               , pointer :: ag
    integer                , pointer :: iro
    logical                , pointer :: wind
    logical                , pointer :: culvert
    logical                , pointer :: mudlay
    logical                , pointer :: nfl
    logical                , pointer :: zmodel
    logical                , pointer :: wavcmp
!
! Global variables
!
    integer                                                                     :: icreep  !  Description and declaration in tricom.igs
    integer                                                       , intent(in)  :: icx     !!  Increment in the X-dir., if ICX= NMAX then computation proceeds in the X-dir. If icx=1 then computation proceeds in the Y-dir.
    integer                                                       , intent(in)  :: icy     !!  Increment in the Y-dir. (see ICX)
    integer                                                                     :: idry
    integer                                                                     :: j       !!  Begin pointer for arrays which have been transformed into 1D arrays. Due to the shift in the 2nd (M-) index, J = -2*NMAX + 1
    integer                                                                     :: kmax    !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: lsecfl  !  Description and declaration in dimens.igs
    integer                                                                     :: lstsci  !  Description and declaration in esm_alloc_int.f90
    integer                                                       , intent(in)  :: mmax    !  Description and declaration in esm_alloc_int.f90
    integer                                                       , intent(in)  :: 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                                                                     :: norow   !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: nsrc    !  Description and declaration in esm_alloc_int.f90
    integer                                                       , intent(in)  :: nst     !!  Time step number
    integer      , dimension(5, norow)                                          :: irocol  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(7, nsrc)                                           :: mnksrc  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: kcs     !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: kcu     !  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)                             :: kfu     !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: kfumax  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: kfumin  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: kfv     !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)                     :: kspu    !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: kadu    !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: kadv    !  Description and declaration in esm_alloc_int.f90
    real(fp)                                                                    :: betac   !  Description and declaration in tricom.igs
    real(fp)     , dimension(12, norow)                                         :: crbc    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(4, norow)                                          :: circ2d  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: a       !!  Internal work array, tridiagonal matrix water levels lower diagonal
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: aa      !!  Internal work array, coupling mean velocity with water level point in (N,M,K) left (down)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: b       !!  Internal work array, tridiagonal matrix water levels main diagonal
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: bb      !!  Internal work array, coefficient mean velocity
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: c       !!  Internal work array, tridiagonal matrix water levels upper diagonal
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: cc      !!  Internal work array, coupling mean velocity with water level point right (upper)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: d       !!  Internal work array, Right Hand side of the Continuity equation
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: d0      !!  Internal work array, Explicit part of the Right Hand side
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: dd      !!  Internal work array, Right hand side of the momentum eq. at (N,M)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: deltau  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: dfu     !  Description and declaration in esm_alloc_real.f90
    real(prec)   , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: dps     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: dpu     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: dteu    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in)  :: evap    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: fcorio  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: fxw     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gsqiu   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gsqs    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gud     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: guu     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gvd     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gvu     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gvv     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: gvz     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: hkru    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: hu0
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: hu      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: hv      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: patm    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in)  :: precip  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: pship   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: rlabda  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: 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)               , intent(in)  :: soumud  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: taubpu  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: taubsu  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: tetau   !!  Factor for upwind approach S0 can be 0.0, 0.5 or 1.0 depending on value of HU, DCO, KSPU and UMEAN
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: tgfsep  !!  Water elev. induced by tide gen.force
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: tp      !  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)                             :: uwtypu  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: vnu2d   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: windsu  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: wsu     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             :: wsbodyu !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)       , intent(out) :: qzk     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)                     :: vicww   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)                     :: w1      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, 3)                          :: cfurou  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, 3)                          :: cfvrou  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: aak     !!  Internal work array (in CUCNP & UZD)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: bbk     !!  Internal work array (in CUCNP & UZD)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: bbka    !!  Internal work array
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: bbkc    !!  Internal work array
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: cck     !!  Internal work array (in CUCNP & UZD)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: d0k     !!  Internal work array
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: ddk     !!  Internal work array, diagonal space at (N,M,K)
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: diapl   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: dpdksi  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: dzu1    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: porosu  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: qxk     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in)  :: qyk     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: rho     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: rnpl    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: rttfu   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: rxx     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: rxy     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: sumrho  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: u0      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: u1      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: ua
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: ub
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: ubrlsu  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: ustokes !  Description and declaration in trisol.igs
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: v1      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax+2)                     :: vicuv   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                       :: dis_nf  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)               :: r0      !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(kmax)                                              :: sig     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(kmax)                                              :: thick   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(kmax, 2, norow)                                    :: circ3d  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrc)                                              :: disch   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrc)                                              :: umdis   !  Description and declaration in esm_alloc_real.f90
    character(1) , dimension(nsrc)                                              :: dismmt  !  Description and declaration in esm_alloc_char.f90
    character(8)                                                                :: dischy  !  Description and declaration in tricom.igs
!
! Local variables
!
    integer       :: ddb
    integer       :: i
    integer       :: icxy
    integer       :: icol
    integer       :: ierror
    integer       :: intdir
    integer       :: iter
    integer       :: itr
    integer       :: k
    integer       :: kenm
    integer       :: kk
    integer       :: m
    integer       :: mmaxddb
    integer       :: n
    integer       :: nhystp
    integer       :: nm
    integer       :: nmaxddb
    integer       :: nmd
    integer       :: nmf
    integer       :: nmlu
    integer       :: nmu
    logical       :: error   ! Flag for detection of closure error in mass-balance
    real(hp)      :: bi
    real(hp)      :: fac
    real(fp)      :: epsomb
    real(fp)      :: hdti
    real(fp)      :: hnm
    real(fp)      :: hucres
    real(fp)      :: humean  ! Mean value for H in U-points
    real(fp)      :: drytrsh
    real(fp)      :: pr
    real(fp)      :: dxiu
    real(fp)      :: dxid
    character(80) :: errtxt
    integer       :: nm_pos ! indicating the array to be exchanged has nm index at the 2nd place, e.g., dbodsd(lsedtot,nm)
!
!! executable statements -------------------------------------------------------
!
    eps         => gdp%gdconst%eps
    maseva      => gdp%gdheat%maseva
    lundia      => gdp%gdinout%lundia
    ntstep      => gdp%gdinttim%ntstep
    dco         => gdp%gdnumeco%dco
    drycrt      => gdp%gdnumeco%drycrt
    dryflc      => gdp%gdnumeco%dryflc
    hdt         => gdp%gdnumeco%hdt
    iter1       => gdp%gdnumeco%iter1
    momsol      => gdp%gdnumeco%momsol
    dpuopt      => gdp%gdnumeco%dpuopt
    rhow        => gdp%gdphysco%rhow
    ag          => gdp%gdphysco%ag
    iro         => gdp%gdphysco%iro
    wind        => gdp%gdprocs%wind
    culvert     => gdp%gdprocs%culvert
    mudlay      => gdp%gdprocs%mudlay
    nfl         => gdp%gdprocs%nfl
    zmodel      => gdp%gdprocs%zmodel
    wavcmp      => gdp%gdprocs%wavcmp
    !
    call timer_start(timer_sud_rest, gdp)
    !
    ddb     = gdp%d%ddbound
    nmaxddb = nmax + 2*gdp%d%ddbound
    mmaxddb = mmax + 2*gdp%d%ddbound
    hdti    = 1.0_fp / hdt
    icxy    = max(icx, icy)
    drytrsh = drycrt
    !
    nm_pos  = 1
    !
    if (idry == 1) then
       !
       ! This is necessary because SUD can be repeated in case of drying
       !    in DRYCHK (DRYFLP <> NO)
       !
       do nm = 1, nmmax
          umean(nm) = 0.0
          if (kfu(nm)==1) then
             do k = 1, kmax
                umean(nm) = umean(nm) + thick(k)*u0(nm, k)
             enddo
          endif
       enddo
       call upwhu   (j         ,nmmaxj    ,nmmax     ,kmax      ,icx       , &
                   & zmodel    ,kcs       ,kcu       ,kspu      ,dps       , &
                   & s0        ,dpu       ,umean     ,hu        ,gdp       )
    endif
    do nm = 1, nmmax
       hu0(nm) = hu(nm)
    enddo
    !
    ! AFTER CALCULATION OF HU FOR CUCNP THE VALUE SHOULD BE CORRECTED
    !     AND TETAU SHOULD BE CALCULATED FOR UPWIND APPROACH
    !
    nmu = +icx
    do nm = 1, nmmax
       nmu       = nmu + 1
       hu(nm)    = max(hu(nm), 0.01_fp)
       tetau(nm) = 0.5
       if (kfu(nm)==1) then
          humean = 0.5*(s0(nm) + s0(nmu)) + dpu(nm)
          if (humean<dco .or. kspu(nm, 0)>0 .or. dpuopt=='UPW' .or. momsol == 'flood ') then
             if (umean(nm)>=0.001) then
                tetau(nm) = 1.0
             elseif (umean(nm)<= - 0.001) then
                tetau(nm) = 0.0
             else
                tetau(nm) = 1.0
                if (s0(nmu)>s0(nm)) tetau(nm) = 0.0
             endif
          endif
       endif
    enddo
    call timer_stop(timer_sud_rest, gdp)
    !
    call timer_start(timer_sud_cucnp, gdp)
    call cucnp(dischy    ,icreep    ,dpdksi    ,s0        ,u0        , &
             & v1        ,w1        ,hu        ,hv        ,dps       ,dpu       , &
             & umean     ,guu       ,gvv       ,gvu       ,gsqs      , &
             & gvd       ,gud       ,gvz       ,gsqiu     ,qxk       , &
             & qyk       ,disch     ,umdis     ,mnksrc    ,dismmt    ,j         , &
             & nmmaxj    ,nmmax     ,kmax      ,icx       ,icy       , &
             & nsrc      ,lsecfl    ,lstsci    ,betac     ,aak       , &
             & bbk       ,cck       ,ddk       ,bbka      ,bbkc      , &
             & thick     ,sig       ,rho       ,sumrho    ,vicuv     , &
             & vnu2d     ,vicww     ,wsu       ,fxw       ,wsbodyu   , &
             & rxx       ,rxy       ,kcs       ,kcu       ,kfu       ,kfv       , &
             & kfs       ,kspu      ,kadu      ,kadv      ,dfu       ,deltau    , &
             & tp        ,rlabda    ,cfurou    ,cfvrou    ,rttfu     , &
             & r0        ,diapl     ,rnpl      ,taubpu    ,taubsu    , &
             & windsu    ,patm      ,fcorio    ,ubrlsu    ,uwtypu    , &
             & hkru      ,pship     ,tgfsep    ,dteu      ,ua        , &
             & ub        ,ustokes   ,.false.   ,u1        ,s1        ,gdp       )
    call timer_stop(timer_sud_cucnp, gdp)
    !
    ! INITIALISATION OF ITERATION OVER CONTINUITY EQUATION
    !
    call timer_start(timer_sud_rest, gdp)
    do k = 1, kmax
       do nm = 1, nmmax
          if (kcs(nm) > 0) then
             d0k(nm, k) = gsqs(nm)*thick(k)*s0(nm)*hdti - qyk(nm, k) + qyk(nm - icy, k)
          endif
       enddo
    enddo
    !
    ! IN LAYER 1 DUE TO PRECIPITATION/EVAPORATION
    !     FOR TIME DEPENDENT INPUT OR HEAT MODEL WITH SPECIAL REQUEST
    !
    if (maseva>0) then
       do nm = 1, nmmax
          if (kcs(nm)==1) then
             d0k(nm, 1) = d0k(nm, 1) + precip(nm)*gsqs(nm)
             if (kfs(nm)==1) then
                d0k(nm, 1) = d0k(nm, 1) - (evap(nm)/rhow)*gsqs(nm)
             endif
          endif
       enddo
    endif
    !
    ! ADDITION OF DISCHARGES (suction only permitted if the point isn't dry)
    !
    do i = 1, nsrc
       nm   = (mnksrc(5, i) + ddb) + ((mnksrc(4, i) - 1) + ddb)*icxy
       k    = mnksrc(6, i)
       if (k .eq. -1) cycle
       kenm = min(1, kfu(nm) + kfu(nm - icx) + kfv(nm) + kfv(nm - icy))
       if (kenm/=0 .or. disch(i)>=0.0) then
          if (k/=0) then
             d0k(nm, k) = d0k(nm, k) + disch(i)
          else
             do kk = 1, kmax
                d0k(nm, kk) = d0k(nm, kk) + disch(i)*thick(kk)
             enddo
          endif
       else
          write (errtxt, '(i0,i3)') nst, i
          call prterr(lundia    ,'S208'    ,trim(errtxt))
       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) then
             if (k/=0) then
                d0k(nm, k) = d0k(nm, k) - disch(i)
             else
                do kk = 1, kmax
                   d0k(nm, kk) = d0k(nm, kk) - disch(i)*thick(kk)
                enddo
             endif
          !
          ! in case of a culvert no warning generated
          !
          elseif (mnksrc(7, i)/=3) then
             write (errtxt, '(i0,i3)') nst, i
             call prterr(lundia    ,'S208'    ,trim(errtxt))
          else
          endif
       endif
    enddo
    !
    ! ADDITION OF DISCHARGES from near field model
    !
    if (nfl) then
       do nm = 1, nmmax
          do k = 1, kmax
             d0k(nm,k) = d0k(nm,k) + dis_nf(nm,k)
          enddo
       enddo
    endif
    !
    ! add sources/sinks mud layer if mudlay == .true.
    !
    if (mudlay) then
       do nm = 1, nmmax
          !
          ! kfs mask array removed, since then the layer cannot increase
          ! any more after it has gone dry
          !
          if (kcs(nm)==1) then
             d0k(nm, 1) = d0k(nm, 1) + gsqs(nm)*soumud(nm)
          endif
       enddo
    endif
    !
    ! Initialise arrays a - dd for all (nm)
    !
    a  = 0.0
    b  = 1.0
    c  = 0.0
    d0 = 0.0
    aa = 0.0
    bb = 1.0
    cc = 0.0
    dd = 0.0
    !
    do nm = 1, nmmax
       d(nm)  = s0(nm)
    enddo
    !
    do k = 1, kmax
       do nm = 1, nmmax
          if (kfu(nm)==1) then
             fac    = real(porosu(nm,k),hp) * real(thick(k),hp) / real(bbk(nm,k),hp)
             aa(nm) = real(aa(nm),hp) + real(aak(nm,k),hp)*fac
             cc(nm) = real(cc(nm),hp) + real(cck(nm,k),hp)*fac
             dd(nm) = real(dd(nm),hp) + real(ddk(nm,k),hp)*fac
          endif
       enddo
    enddo
    do k = 1, kmax
       do nm = 1, nmmax
          if (kcs(nm)==1) d0(nm) = d0(nm) + d0k(nm, k)
       enddo
    enddo
    call timer_stop(timer_sud_rest, gdp)
    !
    ! ITERATIVE LOOP OVER CURRENT ROW. CALCULATION OF CONTINUITY EQ.
    !
 9999 continue
    itr = 0
    !
    do iter = 1, iter1
       !
       ! BOUNDARY CONDITIONS
       !
       call timer_start(timer_sud_cucbp, gdp)
       call cucbp(kmax      ,norow     ,icx       , &
                & icy       ,zmodel    ,irocol    ,kcs       ,kfu       , &
                & kfumin    ,kfumax    ,s0        ,u0        ,dpu       , &
                & hu        ,umean     ,tetau     ,guu       ,gvu       , &
                & dzu1      ,thick     ,circ2d    ,circ3d    ,a         , &
                & b         ,c         ,d         ,aa        ,bb        , &
                & cc        ,dd        ,aak       ,bbk       ,cck       , &
                & ddk       ,crbc      ,wavcmp    ,gdp       )
       call timer_stop(timer_sud_cucbp, gdp)
       !
       ! SET UP SYSTEM OF EQUATIONS FOR INTERIOR POINTS
       !
       call timer_start(timer_sud_rest, gdp)
       if (momsol=='flood ') then
          nmd  = -icx
          do nm = 1, nmmax
             nmd = nmd + 1
             if (kcs(nm)==1) then
                dxid  = hu0(nmd)*guu(nmd)
                dxiu  = hu0(nm) *guu(nm)
                a(nm) = dxid*aa(nmd)
                b(nm) = hdti*gsqs(nm) + dxid*cc(nmd) - dxiu*aa(nm)
                c(nm) = -dxiu*cc(nm)
                d(nm) = d0(nm) - dxiu*dd(nm) + dxid*dd(nmd)
            endif
          enddo
       else
          nmd  = -icx
          do nm = 1, nmmax
             nmd  = nmd  + 1
             if (kcs(nm)==1) then
                a(nm) = guu(nmd)*(hu(nmd)*aa(nmd) - tetau(nmd)*dd(nmd))
                b(nm) = hdti*gsqs(nm)                                             &
                      & + guu(nmd)*(hu(nmd)*cc(nmd) - (1.0 - tetau(nmd))*dd(nmd)) &
                      & - guu(nm) *(hu(nm) *aa(nm)  - tetau(nm)         *dd(nm) )
                c(nm) = -guu(nm)*(hu(nm)*cc(nm) - (1.0 - tetau(nm))*dd(nm))
                d(nm) = d0(nm) - (guu(nm)*dpu(nm)*dd(nm) - guu(nmd)*dpu(nmd)*dd(nmd))
             endif
          enddo
       endif
       call timer_stop(timer_sud_rest, gdp)
       !
       ! Domain decomposition:
       !        Give Mapper chance to build the coupling equations
       !        Note that this is a two-stage process. First, coupling equations
       !        are built for the coupling points start+1;end-1
       !        Secondly, coupling equations are built for the `end' coupling points,
       !        start and end.
       !
       !
       nhystp = nxtstp(d3dflow_build_adi_zeta, gdp)
       !
       ! End of Domain decomposition addition
       !
       !
       !***SCALE ROWS OF MATRIX/RIGHT HAND SIDE VECTOR
       !
       call timer_start(timer_sud_rowsc, gdp)
       do nm = 1, nmmax
          bi    = 1.0_hp / real(b(nm),hp)
          a(nm) = real(a(nm),hp) * bi
          b(nm) = 1.0_fp
          c(nm) = real(c(nm),hp) * bi
          d(nm) = real(d(nm),hp) * bi
       enddo
       call timer_stop(timer_sud_rowsc, gdp)
       !
       ! SOLUTION TRIDIAGONAL SYSTEM FOR THE WATERLEVELS
       !
       if (nhystp==noneighbors) then
          !
          ! Single domain case without domain decomposition
          ! The next piece of code in this IF-statement works for both serial and parallel runs
          ! In case of parallel runs twisted factorization technique is employed which is
          ! perfectly parallizable for two processors only. In case of more than 2 processors,
          ! this technique is combined with the block Jacobi approach at coupling points between
          ! pairs of "twisted" processors. Improvement in convergence is achieved by means of
          ! alternating the pairs of twisted processors at each iteration.
          !
          if ( nproc > 2 ) then
             icol = mod(iter,2)
          else
             icol = 1
          endif
          !
          call timer_start(timer_sud_solve, gdp)
          if ( mod(inode,2) == icol ) then
             !
             ! FORWARD SWEEP (elimination)
             !
             ! Division by the pivot for nmf is not needed anymore
             ! because of row scaling
             !
             do m = 2, mmaxddb
                nm = m*icx - icxy
                do n = 1, nmaxddb
                   nm = nm + icy
                   if (kcs(nm) > 0) then
                      bi    = 1.0_hp / (real(b(nm),hp) - real(a(nm),hp)*real(c(nm-icx),hp))
                      c(nm) = real(c(nm),hp) * bi
                      d(nm) = (real(d(nm),hp) - real(a(nm),hp)*real(d(nm-icx),hp)) * bi
                   endif
                enddo
             enddo
          else
             !
             ! BACKWARD SWEEP (elimination)
             !
             ! Division by the pivot for nmlu is not needed anymore
             ! because of row scaling
             !
             do m = mmaxddb-1, 1, -1
                nm = m*icx - icxy
                do n = 1, nmaxddb
                   nm = nm + icy
                   if (kcs(nm) > 0) then
                      bi    = 1.0_hp / (real(b(nm),hp) - real(c(nm),hp)*real(a(nm+icx),hp))
                      a(nm) = real(a(nm),hp) * bi
                      d(nm) = (real(d(nm),hp) - real(c(nm),hp)*real(d(nm+icx),hp)) * bi
                   endif
                enddo
             enddo
          endif
          !
          ! exchange coefficients a, b, c and d with neighbours for parallel runs
          !
          call dfexchg ( a, 1, 1, dfloat, nm_pos, gdp )
          call dfexchg ( b, 1, 1, dfloat, nm_pos, gdp )
          call dfexchg ( c, 1, 1, dfloat, nm_pos, gdp )
          call dfexchg ( d, 1, 1, dfloat, nm_pos, gdp )
          call dfsync(gdp)
          !
          if ( mod(inode,2) == icol ) then
             !
             ! FORWARD SWEEP in coupling points (elimination)
             !
             do m = 1, mmaxddb
                nm = m*icx - icxy
                do n = 1, nmaxddb
                   nm = nm + icy
                   if (kcs(nm) == -1) then
                      bi     = 1.0_fp / (real(b(nm),hp) - real(a(nm),hp)*real(c(nm-icx),hp))
                      c (nm) = real(c(nm),hp) * bi
                      d (nm) = (real(d(nm),hp) - real(a(nm),hp)*real(d(nm-icx),hp)) * bi
                      s1(nm) = d(nm)
                   endif
                enddo
             enddo
             !
             ! BACKWARD SWEEP (substitution)
             !
             nmlu = mmaxddb*icx - icxy
             do n = 1, nmaxddb
                nmlu = nmlu + icy
                if (kcs(nmlu) > 0) s1(nmlu) = d(nmlu)
             enddo
             do m = mmaxddb - 1, 1, -1
                nm = m*icx - icxy
                do n = 1, nmaxddb
                   nm = nm + icy
                   if (kcs(nm) > 0) then
                      d(nm)  = d(nm) - c(nm)*d(nm + icx)
                      s1(nm) = d(nm)
                  endif
                enddo
             enddo
          else
             !
             ! BACKWARD SWEEP in coupling points (elimination)
             !
             do m = mmaxddb, 1, -1
                nm = m*icx - icxy
                do n = 1, nmaxddb
                   nm = nm + icy
                   if (kcs(nm) == -1) then
                      bi     = 1.0_hp / (real(b(nm),hp) - real(c(nm),hp)*real(a(nm+icx),hp))
                      a (nm) = real(a(nm),hp) * bi
                      d (nm) = (real(d(nm),hp) - real(c(nm),hp)*real(d(nm+icx),hp)) * bi
                      s1(nm) = d(nm)
                   endif
                enddo
             enddo
             !
             ! FORWARD SWEEP (substitution)
             !
             nmf = icx - icxy
             do n = 1, nmaxddb
                nmf = nmf + icy
                if (kcs(nmf) > 0) s1(nmf) = d(nmf)
             enddo
             do m = 2, mmaxddb
                nm = m*icx - icxy
                do n = 1, nmaxddb
                   nm = nm + icy
                   if (kcs(nm) > 0) then
                      d (nm) = d(nm) - a(nm)*d(nm - icx)
                      s1(nm) = d(nm)
                   endif
                enddo
             enddo
          endif
          !
          ! exchange s1 with neighbours for parallel runs
          !
          call dfexchg ( s1, 1, 1, dfloat, nm_pos, gdp )
          !
          ! insert block Jacobi equation in coupling points
          !
          do nm = 1, nmmax
             if ( kcs(nm) == -1 ) then
                a(nm) = 0.0
                b(nm) = 1.0
                c(nm) = 0.0
                d(nm) = s1(nm)
             endif
          enddo
          call timer_stop(timer_sud_solve, gdp)
       else
          !
          ! Domain decomposition:
          !
          ! Wang solver for subdomains
          !
          ! METHOD OF WANG
          !
          ! First part of Wang's algoritm: pre-elimination:
          !
          call timer_start(timer_sud_wangpre, gdp)
          call wangp1(s1        ,kcs       ,irocol    ,norow     ,icx       , &
                    & icy       ,j         ,nmmaxj    ,a         ,b         , &
                    & c         ,d         ,gdp       )
          call timer_stop(timer_sud_wangpre, gdp)
          !
          ! Now second part of Wang's algoritm: elimination of reduced
          ! system. This is carried out by a global mapper process.
          ! At end of global mapper process:
          ! at coupling point at the left side and
          ! at last computational point at the right side we have that
          ! a=c=0, b=1 and d=zeta=new water elevation
          !
          ! Now Mapper builds and solves the reduced system of equation
          !
          if (icx==1) then
             intdir = 0
          else
             !
             ! intdir = 1 corresponds to left_to_Right direction
             ! (see GAWS routines)
             !
             intdir = 1
          endif
          call timer_start(timer_sud_gwsslv, gdp)
          call gwsslv(intdir    )
          call timer_stop(timer_sud_gwsslv, gdp)
          !
          ! Now third part of Wang's algoritm: back substitution
          !
          call timer_start(timer_sud_wangback, gdp)
          call wangp3(s1        ,kcs       ,irocol    ,norow     ,icx       , &
                    & icy       ,j         ,nmmaxj    ,a         ,b         , &
                    & c         ,d         ,gdp       )
          call timer_stop(timer_sud_wangback, gdp)
          !
          ! in case of Hydra (DD method, Wang approach) array d does
          ! not contain the solution, but the right-hand side evaluation.
          ! Since array d (and also s1) is used in the remainder of SUD
          ! for the new water elevation, we have to copy s1 to d.
          !
          do nm = 1, nmmax
             d(nm) = s1(nm)
          enddo
       endif
       !
       ! TOTAL WATERDEPTH DRYING IN SUD
       !
       call timer_start(timer_sud_rest, gdp)
       nmu = +icx
       do nm = 1, nmmax
          nmu = nmu + 1
          if (kfu(nm)==1) then
             !
             ! Special approach for 2D weir points:
             ! - depth value DPU is not corrected like in general (3D) weir case
             ! - crest height is explicitly taken into account in drying check
             !
             hucres = 1E9
             if (abs(kspu(nm, 0))==9) then
                if (umean(nm)>=0.001) then
                   hucres = d(nm) + hkru(nm)
                elseif (umean(nm)<= - 0.001) then
                   hucres = d(nm + icx) + hkru(nm)
                else
                   hucres = max(d(nm + icx), d(nm)) + hkru(nm)
                endif
             endif
             hnm = tetau(nm)*d(nm) + (1. - tetau(nm))*d(nmu) + dpu(nm)
             !
             ! CHECK FOR DRYING
             !
             if (min(hnm, hucres)<=drytrsh) then
                aa(nm)  = 0.0
                bb(nm)  = 1.0
                cc(nm)  = 0.0
                dd(nm)  = 0.0
                kfu(nm) = 0
                itr = 1
             else
                if (momsol=='flood ') then
                   bb(nm) = 1.0
                else
                   bb(nm) = hu(nm)/hnm
                endif
                hu(nm) = hnm
             endif
          endif
       enddo
       call timer_stop(timer_sud_rest, gdp)
       !
       ! determine global maximum of 'itr' over all nodes
       ! Note: this enables to synchronize the repeating computation
       !
       call dfreduce_gdp( itr, 1, dfint, dfmax, gdp )
       !
       ! REPEAT COMPUTATION IF POINT IS SET DRY
       !       FIRST RESET HU
       !
       ! Domain decomposition:
       !    Synchronize on drying before finishing solve zeta
       !
       ! Note that if iter<iter1 flow goes always back to Build step
       ! (either because of DD flag or because of iter loop until iter=iter1)
       ! Since no mapping occurs for check_sud_dry, this communication
       ! step could be skipped for iter<iter1.
       !
       nhystp = nxtdry(d3dflow_check_sud_dry, itr, gdp)
       !
       ! repeat computation if point is set dry
       !
       if (nhystp==d3dflow_build_adi_zeta .or. &
         & (nhystp==noneighbors .and. itr==1)) then
          !
          ! End of Domain decomposition addition
          !
          !
          do nm = 1, nmmax
             hu(nm) = hu0(nm)
          enddo
          goto 9999
       endif
    enddo
    !
    ! exchange kfu with neighbours for parallel runs
    !
    call dfexchg ( kfu, 1, 1, dfint, nm_pos, gdp )
    !
    ! adapt discharge boundary conditions
    !
    call timer_start(timer_sud_cucdp, gdp)
    call cucdp(kfu       ,irocol    ,norow     ,j         ,nmmaxj    , &
             & icx       ,icy       ,bb        ,gdp       )
    call timer_stop(timer_sud_cucdp, gdp)
    !
    ! COMPUTATION (LAYER & DEPTH AVERAGED) VELOCITIES AND DISCHARGES
    !
    call timer_start(timer_sud_veldisch, gdp)
    do k = 1, kmax
       do nm = 1, nmmax
          if (kfu(nm)==1) then
             pr        = aak(nm, k)*d(nm) + cck(nm, k)*d(nm + icx)
             u1(nm, k) = (ddk(nm, k) - bb(nm)*pr)/bbk(nm, k)
          else
             u1(nm, k) = 0.0
          endif
       enddo
    enddo
    !
    ! exchange u1 with neighbours for parallel runs
    !
    call dfexchg ( u1, 1, kmax, dfloat, nm_pos, gdp )
    !
    ! compute horizontal discharge
    !
    do k = 1, kmax
       do nm = 1, nmmax
          if (momsol == 'flood') then
             qxk(nm, k) = guu(nm)*hu0(nm)*thick(k)*u1(nm, k)*porosu(nm,k)
          else
             qxk(nm, k) = guu(nm)*hu (nm)*thick(k)*u1(nm, k)*porosu(nm,k)
          endif
       enddo
    enddo
    !
    ! adapt waterlevels and velocities at coupling boundaries to
    ! prevent mass closure error in case of parallel runs
    ! Note: with 1 or 2 processors there is no need for this adaption
    !
    if ( nproc > 2 ) call dfmassc (s1        ,u1        ,qxk       ,hu        ,d0        , &
                                 & dpu       ,porosu    ,gsqs      ,guu       ,tetau     , &
                                 & kcs       ,kcu       ,kfu       ,thick     ,nmmax     , &
                                 & kmax      ,icx       ,gdp )
    !
    ! compute depth-averaged velocity
    !
    do nm = 1, nmmax
       umean(nm) = 0.0
    enddo
    do k = 1, kmax
       do nm = 1, nmmax
          umean(nm) = umean(nm) + thick(k)*u1(nm, k)
       enddo
    enddo
    !
    ! Domain decomposition:
    !
    nhystp = nxtstp(d3dflow_finish_wang, gdp)
    !
    ! End of Domain decomposition addition
    !
    if (kmax>1) then
       !
       ! COMPUTATION VERTICAL VELOCITIES AND DISCHARGES
       !
       ! Initialise arrays qzk and w1 for all (nm,k)
       !
       qzk = 0.0
       w1  = 0.0
       !
       do k = 1, kmax
          do nm = 1, nmmax
             if (kcs(nm)==1) then
                w1(nm, k) = w1(nm, k - 1) + thick(k)*s1(nm)*hdti        &
                          & + (qxk(nm, k) - qxk(nm - icx, k)            &
                          &    - d0k(nm, k)                 ) / gsqs(nm)
                qzk(nm, k) = w1(nm, k)*gsqs(nm)
             endif
          enddo
       enddo
       !
       ! exchange w1 with neighbours for parallel runs
       !
       call dfexchg ( w1, 0, kmax, dfloat, nm_pos, gdp )
       !
       ! compute vertical discharge
       !
       !
       epsomb = max(eps, eps*hdt)
       !
       error = .false.
       do nm = 1, nmmax
          if (abs(w1(nm, kmax))>epsomb) then
             error = .true.
             w1(nm, kmax) = 0.0
          endif
       enddo
       ierror = 0
       if (error) ierror = 1
       call dfreduce_gdp( ierror, 1, dfint, dfmax, gdp )
       error = ierror==1
       if (error) then
          write (errtxt, '(a,e12.3,a,i0,a)') 'Mass closure error exceeds ', &
               & epsomb, ' after ', ntstep, ' timesteps.'
          call prterr(lundia, 'U190', trim(errtxt))
       endif
    endif
    call timer_stop(timer_sud_veldisch, gdp)
    !
    ! Optionally compute individual momentum terms for output
    !
    if (gdp%gdflwpar%flwoutput%momentum) then
       ! pass hu or hu0 ?
       call timer_start(timer_sud_cucnp, gdp)
       call cucnp(dischy    ,icreep    ,dpdksi    ,s0        ,u0        , &
                & v1        ,w1        ,hu        ,hv        ,dps       ,dpu       , &
                & umean     ,guu       ,gvv       ,gvu       ,gsqs      , &
                & gvd       ,gud       ,gvz       ,gsqiu     ,qxk       , &
                & qyk       ,disch     ,umdis     ,mnksrc    ,dismmt    ,j         , &
                & nmmaxj    ,nmmax     ,kmax      ,icx       ,icy       , &
                & nsrc      ,lsecfl    ,lstsci    ,betac     ,aak       , &
                & bbk       ,cck       ,ddk       ,bbka      ,bbkc      , &
                & thick     ,sig       ,rho       ,sumrho    ,vicuv     , &
                & vnu2d     ,vicww     ,wsu       ,fxw       ,wsbodyu   , &
                & rxx       ,rxy       ,kcs       ,kcu       ,kfu       ,kfv       , &
                & kfs       ,kspu      ,kadu      ,kadv      ,dfu       ,deltau    , &
                & tp        ,rlabda    ,cfurou    ,cfvrou    ,rttfu     , &
                & r0        ,diapl     ,rnpl      ,taubpu    ,taubsu    , &
                & windsu    ,patm      ,fcorio    ,ubrlsu    ,uwtypu    , &
                & hkru      ,pship     ,tgfsep    ,dteu      ,ua        , &
                & ub        ,ustokes   ,.true.    ,u1        ,s1        ,gdp       )
       call timer_stop(timer_sud_cucnp, gdp)
    endif
end subroutine sud

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


module two_point_flux_mod

use globals_mod

private

public :: xflux, yflux

contains

pure elemental function logmean(a,b) result(r)

    real(dp), intent(in)    :: a,b
    real(dp)                :: x,u,r

    real(dp), parameter :: eps = 1e-5
    real(dp), parameter :: c1 = 1.0/6.0
    real(dp), parameter :: c2 = 2.0/45.0
    real(dp), parameter :: c3 = 22.0/945.0

    x = a/b

    if (abs(a-b) < eps) then
        u  = (x*(x-2.0)+1.0)/(x*(x+2.0)+1.0)
        r = (a+b)*(0.5 - u*(c1 - u*(c2 - c3*u)))
    else
        r = (a-b)/log(x)
    endif

end function

pure function xflux(u,v) result(f)

    use equations_mod
    use config_mod, only: kappa

    real(dp), intent(in)    :: u(N_VARS)
    real(dp), intent(in)    :: v(N_VARS)
    real(dp)                :: f(N_VARS)

    real(dp)                :: dens,pres,velx,vely,enth
    real(dp)                :: bu,bv !! beta
    real(dp)                :: pu,pv !! pressure
    real(dp)                :: vu, wu, vv, wv !! velocity

    pu = pressure(u)
    pv = pressure(v)

    bu = 0.5*u(DENS_VAR)/pu
    bv = 0.5*v(DENS_VAR)/pv

    vu = u(MOMX_VAR)/u(DENS_VAR)
    wu = u(MOMY_VAR)/u(DENS_VAR)

    vv = v(MOMX_VAR)/v(DENS_VAR)
    wv = v(MOMY_VAR)/v(DENS_VAR)

    dens = logmean(u(DENS_VAR),v(DENS_VAR))
    velx = 0.5*(vu + vv)
    vely = 0.5*(wu + wv)
    pres = 0.5*(u(DENS_VAR)+v(DENS_VAR))/(bu+bv)

    enth = 0.5/(kappa-1)/logmean(bu,bv) &
         + pres/dens + velx*velx + vely*vely &
         - 0.25*(vu*vu + wu*wu + vv*vv + wv*wv)

    f(DENS_VAR) = dens * velx
    f(MOMX_VAR) = f(DENS_VAR) * velx + pres
    f(MOMY_VAR) = f(DENS_VAR) * vely
    f(ENER_VAR) = f(DENS_VAR) * enth

end function

pure function yflux(uL,uR) result(f)

    use equations_mod, only: N_VARS
    use equations_mod, only: DENS_VAR
    use equations_mod, only: MOMX_VAR
    use equations_mod, only: MOMY_VAR
    use equations_mod, only: ENER_VAR

    real(dp), intent(in) :: uL(N_VARS), uR(N_VARS)
    real(dp)             :: f(N_VARS)

    real(dp) :: ruL(N_VARS), ruR(N_VARS), rf(N_VARS) !! rotated

    ruL(DENS_VAR) = uL(DENS_VAR)
    ruL(MOMX_VAR) = uL(MOMY_VAR)
    ruL(MOMY_VAR) = uL(MOMX_VAR)
    ruL(ENER_VAR) = uL(ENER_VAR)

    ruR(DENS_VAR) = uR(DENS_VAR)
    ruR(MOMX_VAR) = uR(MOMY_VAR)
    ruR(MOMY_VAR) = uR(MOMX_VAR)
    ruR(ENER_VAR) = uR(ENER_VAR)

    rf = xflux(ruL,ruR)

    f(DENS_VAR) = rf(DENS_VAR)
    f(MOMX_VAR) = rf(MOMY_VAR)
    f(MOMY_VAR) = rf(MOMX_VAR)
    f(ENER_VAR) = rf(ENER_VAR)

end function

end module


text:
      SUBROUTINE COST (N,X,WSAVE)
      IMPLICIT NONE
      INTEGER :: N, I, K, KC, MODN, NM1, NP1, NS2
      REAL :: X(*), WSAVE(*), C1, T1, T2, TX2, X1H, X1P3, XI, XIM2
      NM1 = N-1
      NP1 = N+1
      NS2 = N/2
      IF (N.lt.2) GO TO 106
      IF (N.eq.2) GO TO 101
      GO TO 102
  101 X1H = X(1)+X(2)
      X(2) = X(1)-X(2)
      X(1) = X1H
      RETURN
  102 IF (N .GT. 3) GO TO 103
      X1P3 = X(1)+X(3)
      TX2 = X(2)+X(2)
      X(2) = X(1)-X(3)
      X(1) = X1P3+TX2
      X(3) = X1P3-TX2
      RETURN
  103 C1 = X(1)-X(N)
      X(1) = X(1)+X(N)
      DO 104 K=2,NS2
         KC = NP1-K
         T1 = X(K)+X(KC)
         T2 = X(K)-X(KC)
         C1 = C1+WSAVE(KC)*T2
         T2 = WSAVE(K)*T2
         X(K) = T1-T2
         X(KC) = T1+T2
  104 CONTINUE
      MODN = MOD(N,2)
      IF (MODN .NE. 0) X(NS2+1) = X(NS2+1)+X(NS2+1)
      CALL RFFTF (NM1,X,WSAVE(N+1))
      XIM2 = X(2)
      X(2) = C1
      DO 105 I=4,N,2
         XI = X(I)
         X(I) = X(I-2)-X(I-1)
         X(I-1) = XIM2
         XIM2 = XI
  105 CONTINUE
      IF (MODN .NE. 0) X(N) = XIM2
  106 RETURN
      END

max_stars_repo_name: magnusja/scipy

max_issues_repo_name: magnusja/scipy

max_forks_repo_name: magnusja/scipy


text:
      MODULE setmcp_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.3E  06:40:03  12/28/06
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  11/02/17
      SUBROUTINE setmcp (AVAIL)
      LOGICAL, INTENT(OUT) :: AVAIL
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
*DECK SXLCAL
      SUBROUTINE SXLCAL (N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM,
     +   WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR, NELT, IA, JA, A,
     +   ISYM)
C***BEGIN PROLOGUE  SXLCAL
C***SUBSIDIARY
C***PURPOSE  Internal routine for SGMRES.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2A4, D2B4
C***TYPE      SINGLE PRECISION (SXLCAL-S, DXLCAL-D)
C***KEYWORDS  GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
C             NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
C***AUTHOR  Brown, Peter, (LLNL), pnbrown@llnl.gov
C           Hindmarsh, Alan, (LLNL), alanh@llnl.gov
C           Seager, Mark K., (LLNL), seager@llnl.gov
C             Lawrence Livermore National Laboratory
C             PO Box 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C***DESCRIPTION
C        This  routine computes the solution  XL,  the current SGMRES
C        iterate, given the  V(I)'s and  the  QR factorization of the
C        Hessenberg  matrix HES.   This routine  is  only called when
C        ITOL=11.
C
C *Usage:
C      INTEGER N, LGMR, MAXLP1, JSCAL, JPRE, NMSL, IPAR(USER DEFINED)
C      INTEGER NELT, IA(NELT), JA(NELT), ISYM
C      REAL X(N), XL(N), ZL(N), HES(MAXLP1,MAXL), Q(2*MAXL),
C     $     V(N,MAXLP1), R0NRM, WK(N), SZ(N), RPAR(USER DEFINED),
C     $     A(NELT)
C      EXTERNAL MSOLVE
C
C      CALL SXLCAL(N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM,
C     $     WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR,
C     $     NELT, IA, JA, A, ISYM)
C
C *Arguments:
C N      :IN       Integer
C         The order of the matrix A, and the lengths
C         of the vectors SR, SZ, R0 and Z.
C LGMR   :IN       Integer
C         The number of iterations performed and
C         the current order of the upper Hessenberg
C         matrix HES.
C X      :IN       Real X(N)
C         The current approximate solution as of the last restart.
C XL     :OUT      Real XL(N)
C         An array of length N used to hold the approximate
C         solution X(L).
C         Warning: XL and ZL are the same array in the calling routine.
C ZL     :IN       Real ZL(N)
C         An array of length N used to hold the approximate
C         solution Z(L).
C HES    :IN       Real HES(MAXLP1,MAXL)
C         The upper triangular factor of the QR decomposition
C         of the (LGMR+1) by LGMR upper Hessenberg matrix whose
C         entries are the scaled inner-products of A*V(*,i) and V(*,k).
C MAXLP1 :IN       Integer
C         MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
C         MAXL is the maximum allowable order of the matrix HES.
C Q      :IN       Real Q(2*MAXL)
C         A real array of length 2*MAXL containing the components
C         of the Givens rotations used in the QR decomposition
C         of HES.  It is loaded in SHEQR.
C V      :IN       Real V(N,MAXLP1)
C         The N by(LGMR+1) array containing the LGMR
C         orthogonal vectors V(*,1) to V(*,LGMR).
C R0NRM  :IN       Real
C         The scaled norm of the initial residual for the
C         current call to SPIGMR.
C WK     :IN       Real WK(N)
C         A real work array of length N.
C SZ     :IN       Real SZ(N)
C         A vector of length N containing the non-zero
C         elements of the diagonal scaling matrix for Z.
C JSCAL  :IN       Integer
C         A flag indicating whether arrays SR and SZ are used.
C         JSCAL=0 means SR and SZ are not used and the
C                 algorithm will perform as if all
C                 SR(i) = 1 and SZ(i) = 1.
C         JSCAL=1 means only SZ is used, and the algorithm
C                 performs as if all SR(i) = 1.
C         JSCAL=2 means only SR is used, and the algorithm
C                 performs as if all SZ(i) = 1.
C         JSCAL=3 means both SR and SZ are used.
C JPRE   :IN       Integer
C         The preconditioner type flag.
C MSOLVE :EXT      External.
C         Name of the routine which solves a linear system Mz = r for
C         z given r with the preconditioning matrix M (M is supplied via
C         RPAR and IPAR arrays.  The name of the MSOLVE routine must
C         be declared external in the calling program.  The calling
C         sequence to MSOLVE is:
C             CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RPAR, IPAR)
C         Where N is the number of unknowns, R is the right-hand side
C         vector and Z is the solution upon return.  NELT, IA, JA, A and
C         ISYM are defined as below.  RPAR is a real array that can be
C         used to pass necessary preconditioning information and/or
C         workspace to MSOLVE.  IPAR is an integer work array for the
C         same purpose as RPAR.
C NMSL   :IN       Integer
C         The number of calls to MSOLVE.
C RPAR   :IN       Real RPAR(USER DEFINED)
C         Real workspace passed directly to the MSOLVE routine.
C IPAR   :IN       Integer IPAR(USER DEFINED)
C         Integer workspace passed directly to the MSOLVE routine.
C NELT   :IN       Integer
C         The length of arrays IA, JA and A.
C IA     :IN       Integer IA(NELT)
C         An integer array of length NELT containing matrix data.
C         It is passed directly to the MATVEC and MSOLVE routines.
C JA     :IN       Integer JA(NELT)
C         An integer array of length NELT containing matrix data.
C         It is passed directly to the MATVEC and MSOLVE routines.
C A      :IN       Real A(NELT)
C         A real array of length NELT containing matrix data.
C         It is passed directly to the MATVEC and MSOLVE routines.
C ISYM   :IN       Integer
C         A flag to indicate symmetric matrix storage.
C         If ISYM=0, all non-zero entries of the matrix are
C         stored.  If ISYM=1, the matrix is symmetric and
C         only the upper or lower triangular part is stored.
C
C***SEE ALSO  SGMRES
C***ROUTINES CALLED  SAXPY, SCOPY, SHELS
C***REVISION HISTORY  (YYMMDD)
C   871001  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   910502  Removed MSOLVE from ROUTINES CALLED list.  (FNF)
C   910506  Made subsidiary to SGMRES.  (FNF)
C   920511  Added complete declaration section.  (WRB)
C***END PROLOGUE  SXLCAL
C         The following is for optimized compilation on LLNL/LTSS Crays.
CLLL. OPTIMIZE
C     .. Scalar Arguments ..
      REAL R0NRM
      INTEGER ISYM, JPRE, JSCAL, LGMR, MAXLP1, N, NELT, NMSL
C     .. Array Arguments ..
      REAL A(NELT), HES(MAXLP1,*), Q(*), RPAR(*), SZ(*), V(N,*), WK(N),
     +     X(N), XL(N), ZL(N)
      INTEGER IA(NELT), IPAR(*), JA(NELT)
C     .. Subroutine Arguments ..
      EXTERNAL MSOLVE
C     .. Local Scalars ..
      INTEGER I, K, LL, LLP1
C     .. External Subroutines ..
      EXTERNAL SAXPY, SCOPY, SHELS
C***FIRST EXECUTABLE STATEMENT  SXLCAL
      LL = LGMR
      LLP1 = LL + 1
      DO 10 K = 1,LLP1
         WK(K) = 0
 10   CONTINUE
      WK(1) = R0NRM
      CALL SHELS(HES, MAXLP1, LL, Q, WK)
      DO 20 K = 1,N
         ZL(K) = 0
 20   CONTINUE
      DO 30 I = 1,LL
         CALL SAXPY(N, WK(I), V(1,I), 1, ZL, 1)
 30   CONTINUE
      IF ((JSCAL .EQ. 1) .OR.(JSCAL .EQ. 3)) THEN
         DO 40 K = 1,N
            ZL(K) = ZL(K)/SZ(K)
 40      CONTINUE
      ENDIF
      IF (JPRE .GT. 0) THEN
         CALL SCOPY(N, ZL, 1, WK, 1)
         CALL MSOLVE(N, WK, ZL, NELT, IA, JA, A, ISYM, RPAR, IPAR)
         NMSL = NMSL + 1
      ENDIF
C         calculate XL from X and ZL.
      DO 50 K = 1,N
         XL(K) = X(K) + ZL(K)
 50   CONTINUE
      RETURN
C------------- LAST LINE OF SXLCAL FOLLOWS ----------------------------
      END


text:
    subroutine getrusanovflux(Qfl,Qfr,Fnl,Fnr,Az,normf,sign_l,sign_r,eigv)
    use setup2d,only: gam,eigvmax
    implicit none

    integer,intent(in) :: sign_l,sign_r
    double precision :: Qfl(9),Qfr(9),Fnl(8),Fnr(8),eigv,normf(2)
    double precision :: ul,vl,wl,ur,vr,wr,Vn_l,Vn_r,pl,pr,&
    Bxl,Byl,Bzl,Bxr,Byr,Bzr,NORM_Bl,NORM_Br,Bx,By,Bz,Bn_l,Bn_r,Az
    double precision :: rhoi,c_a,c_d,c_n,c_tmp,c_f
    double precision :: unit_x,unit_y,magnorm

    ! get the magnitude of normal vector |J| div(\xi) or |J| div(\eta) from left cell
    magnorm = dsqrt(normf(1)**2+normf(2)**2)

    ! unit normal vector pointing from left cell to right cell
    unit_x = normf(1)/magnorm*sign_l
    unit_y = normf(2)/magnorm*sign_l

    ! compute the Rosanov flux on physical domain along vector (unit_x,unit_y)

    !for the left side
    rhoi  = 1.d0/Qfl(1)
    ul    = rhoi*Qfl(2)
    vl    = rhoi*Qfl(3)
    wl    = rhoi*Qfl(4)
    Bxl   = Qfl(6)
    Byl   = Qfl(7)
    Bzl   = Qfl(8)

    NORM_Bl = Bxl**2+Byl**2+Bzl**2

    pl    = (gam-1.d0)*(Qfl(5) - 0.5d0*Qfl(1)*(ul**2+vl**2+wl**2)-NORM_Bl*0.5d0)
    pl    = pl + 0.5d0*NORM_Bl  ! Total pressure = pressure + magnetic pressure
    Vn_l  = ul*unit_x+vl*unit_y
    Bn_l  = Bxl*unit_x + Byl*unit_y

    !for the right side
    rhoi  = 1.d0/Qfr(1)
    ur    = rhoi*Qfr(2)
    vr    = rhoi*Qfr(3)
    wr    = rhoi*Qfr(4)
    Bxr   = Qfr(6)
    Byr   = Qfr(7)
    Bzr   = Qfr(8)

    NORM_Br = Bxr**2 + Byr**2 + Bzr**2

    pr    = (gam-1.d0)*(Qfr(5) - 0.5d0*Qfr(1)*(ur**2+vr**2+wr**2)-NORM_Br*0.5d0)
    pr    = pr + NORM_Br*0.5d0  ! Total pressure = pressure + magnetic pressure
    Vn_r  = ur*unit_x+vr*unit_y
    Bn_r  = Bxr*unit_x + Byr*unit_y

    ! calculate the normal flux component at left side
    Fnl(1) = Qfl(1)*Vn_l
    Fnl(2) = Qfl(2)*Vn_l + pl*unit_x - Bxl*Bn_l
    Fnl(3) = Qfl(3)*Vn_l + pl*unit_y - Byl*Bn_l
    Fnl(4) = Qfl(4)*Vn_l -Bzl*Bn_l
    Fnl(5) = (Qfl(5)+pl)*Vn_l - Bn_l*(ul*Bxl+vl*Byl+wl*Bzl)
    Fnl(6) = (vl*Bxl-ul*Byl)*unit_y
    Fnl(7) = (ul*Byl-vl*Bxl)*unit_x
    Fnl(8) = Bzl*Vn_l - wl*Bn_l

    ! calculate the normal flux component at right side
    Fnr(1) = Qfr(1)*Vn_r
    Fnr(2) = Qfr(2)*Vn_r + pr*unit_x - Bxr*Bn_r
    Fnr(3) = Qfr(3)*Vn_r + pr*unit_y - Byr*Bn_r
    Fnr(4) = Qfr(4)*Vn_r - Bzr*Bn_r
    Fnr(5) = (Qfr(5)+pr)*Vn_r - Bn_r*(ur*Bxr+vr*Byr+wr*Bzr)
    Fnr(6) = (vr*Bxr-ur*Byr)*unit_y
    Fnr(7) = (ur*Byr-vr*Bxr)*unit_x
    Fnr(8) = Bzr*Vn_r - wr*Bn_r

    ! calculate the eignvalue (fast magneto-acoustic wave speed)
    pl = pl -  0.5d0*NORM_Bl
    pr = pr -  0.5d0*NORM_Br
    Bx = 0.5d0*(Bxl+Bxr)
    By = 0.5d0*(Byl+Byr)
    Bz = 0.5d0*(Bzl+Bzr)

    c_a  = gam*(pl+pr)/(Qfl(1)+Qfr(1))
    c_d  = (Bx**2+By**2+Bz**2)/(0.5d0*(Qfl(1)+Qfr(1)))
    c_n  = (Bx*unit_x+By*unit_y)**2/(0.5d0*(Qfl(1)+Qfr(1)))
    c_tmp= dsqrt((c_a+c_d)**2-4.0d0*c_a*c_n)

    c_f  = dsqrt(0.5d0*(c_a+c_d+c_tmp))

    eigv = 0.5d0*dabs(Vn_l+Vn_r) + c_f

    ! Rusanov flux for left cell on physical domain
    Fnl(1:8) = 0.5d0*(Fnl(1:8)+Fnr(1:8)-eigv*(Qfr(1:8)-Qfl(1:8)))

    ! Transfer to computational domain for left cell
    Fnl(1:8) = magnorm*sign_l*Fnl(1:8)

    ! Transfer to computational domain for right cell
    Fnr(1:8) = sign_l*sign_r*Fnl(1:8)

    ! upwind scheme for magnetic potential A
    if (Vn_l.gt.0.d0) then
        Az = Qfl(9)
    else
        Az = Qfr(9)
    end if

    if (eigv.gt.eigvmax) then
        eigvmax = eigv
    end if

    end subroutine getrusanovflux

    subroutine getIfluxvectors(Qi,Fi,Gi)
    use setup2d,only: gam
    implicit none

    double precision,intent(in) :: Qi(8)
    double precision,intent(out) :: Fi(8),Gi(8)
    double precision :: rho,rhou,rhov,rhow,rhoE,pr,Bx,By,Bz,NORM_B

    rho  = Qi(1)
    rhou = Qi(2)
    rhov = Qi(3)
    rhow = Qi(4)
    rhoE = Qi(5)
    Bx   = Qi(6)
    By   = Qi(7)
    Bz   = Qi(8)

    NORM_B = Bx**2+By**2+Bz**2

    !compute the pressure
    pr = (gam-1.d0)*(rhoE - 0.5d0*(rhou**2+rhov**2+rhow**2)/rho-0.5d0*NORM_B)

    !compute flux vector F
    Fi(1) = rhou
    Fi(2) = (rhou**2)/rho + pr + (0.5d0*NORM_B-Bx**2)
    Fi(3) = rhou*rhov/rho - Bx*By
    Fi(4) = rhou*rhow/rho - Bx*Bz
    Fi(5) = (rhou/rho)*(rhoE+pr+0.5d0*NORM_B) - (rhou/rho*Bx+rhov/rho*By+rhow/rho*Bz)*Bx
    Fi(6) = 0.d0
    Fi(7) = rhou/rho*By-rhov/rho*Bx
    Fi(8) = rhou/rho*Bz-rhow/rho*Bx

    !compute flux vector G
    Gi(1) = rhov
    Gi(2) = rhou*rhov / rho - Bx*By
    Gi(3) = (rhov**2)/rho + pr + (0.5d0*NORM_B-By**2)
    Gi(4) = rhow*rhov/rho - By*Bz
    Gi(5) = (rhov/rho)*(rhoE+pr+0.5d0*NORM_B) - (rhou/rho*Bx+rhov/rho*By+rhow/rho*Bz)*By
    Gi(6) = rhov/rho*Bx - rhou/rho*By
    Gi(7) = 0.d0
    Gi(8) = rhov/rho*Bz - rhow/rho*By

    end subroutine getIfluxvectors

    !subroutine to calculate viscous flux
    !Qv and nablaQ are all smoothed ones for the interface flux points
    SUBROUTINE getVfluxvectors(Qv,nablaQ,Fv,Gv)

    use setup2d
    IMPLICIT NONE

    double precision,intent(in) :: Qv(numv),nablaQ(numv,2)
    double precision,intent(out) :: Fv(numv),Gv(numv)
    double precision :: u, v, w, NORM_U, Bx, By, Bz, Jx, Jy, Jz, NORM_B, inte
    double precision :: rho, rhou, rhov, rhow, rhoE
    double precision :: ux, vx, uy, vy, wx, wy, ex, ey, rhox, rhoy, Bxy, Byx, Bzy, Bzx
    double precision :: du(2),dv(2),dw(2),drho(2),dBx(2),dBy(2),dBz(2),de(2),dKE(2),dpre_B(2)

    rho  = Qv(1)
    rhou = Qv(2)
    rhov = Qv(3)
    rhow = Qv(4)
    rhoE = Qv(5)
    Bx   = Qv(6)
    By   = Qv(7)
    Bz   = Qv(8)

    u = rhou/rho
    v = rhov/rho
    w = rhow/rho

    NORM_U = u**2  + v**2  + w**2
    NORM_B = Bx**2 + By**2 + Bz**2

    inte = (rhoE - 0.5d0*NORM_U*rho - 0.5d0*NORM_B)/rho

    drho(:) =  nablaQ(1,:)
    du(:)   = (nablaQ(2,:)-drho(:)*u)/rho
    dv(:)   = (nablaQ(3,:)-drho(:)*v)/rho
    dw(:)   = (nablaQ(4,:)-drho(:)*w)/rho
    dBx(:)  =  nablaQ(6,:)
    dBy(:)  =  nablaQ(7,:)
    dBz(:)  =  nablaQ(8,:)

    dKE(:)    = 0.5d0*NORM_U*drho(:) + rho*(u*du(:)+v*dv(:)+w*dw(:))
    dpre_B(:) = Bx*dBx(:) + By*dBy(:) + Bz*dBz(:)

    de(:)  = (nablaQ(5,:) - dKE(:) - dpre_B(:) - drho(:)*inte)/rho

    rhox = drho(1)
    rhoy = drho(2)

    ux = du(1)
    uy = du(2)
    vx = dv(1)
    vy = dv(2)
    wx = dw(1)
    wy = dw(2)
    ex = de(1)
    ey = de(2)

    Bxy = dBx(2)
    Byx = dBy(1)
    Bzx = dBz(1)
    Bzy = dBz(2)

    ! Viscous term from hydrodynamics
    Fv(1) = 0.d0
    Fv(2) = mu*(2.d0*ux+lambda*(ux+vy))
    Fv(3) = mu*(vx+uy)
    Fv(4) = mu*wx
    Fv(5) = u*Fv(2)+v*Fv(3)+w*Fv(4)+mu*gam*ex/prandt

    Gv(1) = 0.d0
    Gv(2) = mu*(vx+uy)
    Gv(3) = mu*(2.d0*vy+lambda*(ux+vy))
    Gv(4) = mu*wy
    Gv(5) = Gv(2)*u+Gv(3)*v+Gv(4)*w+mu*gam*ey/prandt

    ! Resistive term from magnetic

    Jx = Bzy   - 0.0d0
    Jy = 0.0d0 - Bzx
    Jz = Byx   - Bxy

    Fv(6) = 0.0d0
    Fv(7) = eta0*Jz
    Fv(8) = -eta0*Jy
    Fv(5) = Fv(5) + Fv(7)*By + Fv(8)*Bz

    Gv(6) = -eta0*Jz
    Gv(7) = 0.0d0
    Gv(8) = eta0*Jx
    Gv(5) = Gv(5) + Gv(6)*Bx + Gv(8)*Bz

    Fv(9) = 0.d0
    Gv(9) = 0.d0

    END SUBROUTINE getVfluxvectors

    SUBROUTINE getVfluxvectors_AV(Qv,nablaQ,Fv,Gv,fam,ifp,jfp,ic)
    use setup2d, only: numv,muAVfi,muAVfj,Xs,Xf,prandt,gam,lambda

    double precision :: Qv(numv),nablaQ(numv,2),Fv(numv),Gv(numv)
    integer,intent(in) :: ifp,jfp,ic,fam
    double precision :: epsilon_loc(numv),xsi,eta,mu_art,eta_art
    double precision :: u, v, w, NORM_U, Bx, By, Bz, NORM_B, inte
    double precision :: rho, rhou, rhov, rhow, rhoE
    double precision :: ux, vx, uy, vy, wx, wy, ex, ey, rhox, rhoy, Bxy, Byx, Bzy, Bzx
    double precision :: Jx, Jy, Jz
    double precision :: du(2),dv(2),dw(2),drho(2),dBx(2),dBy(2),dBz(2),de(2),dKE(2),dpre_B(2)

    if (fam==1) then
        xsi = Xf(ifp)
        eta = Xs(jfp)
    else
        xsi = Xs(ifp)
        eta = Xf(jfp)
    end if

    if (fam==1) then
        epsilon_loc(1:numv) = muAVfi(1:numv,1,ic) * (1.d0-xsi) + &
        muAVfi(1:numv,2,ic) * xsi
    else
        epsilon_loc(1:numv) = muAVfj(1:numv,1,ic) * (1.d0-eta) + &
        muAVfj(1:numv,2,ic) * eta
    end if

    mu_art = epsilon_loc(1)
    eta_art = 0.d0

    rho  = Qv(1)
    rhou = Qv(2)
    rhov = Qv(3)
    rhow = Qv(4)
    rhoE = Qv(5)
    Bx   = Qv(6)
    By   = Qv(7)
    Bz   = Qv(8)

    u = rhou/rho
    v = rhov/rho
    w = rhow/rho

    NORM_U = u**2  + v**2  + w**2
    NORM_B = Bx**2 + By**2 + Bz**2

    inte = (rhoE - 0.5d0*NORM_U*rho - 0.5d0*NORM_B)/rho

    drho(:) =  nablaQ(1,:)
    du(:)   = (nablaQ(2,:)-drho(:)*u)/rho
    dv(:)   = (nablaQ(3,:)-drho(:)*v)/rho
    dw(:)   = (nablaQ(4,:)-drho(:)*w)/rho
    dBx(:)  =  nablaQ(6,:)
    dBy(:)  =  nablaQ(7,:)
    dBz(:)  =  nablaQ(8,:)

    dKE(:)    = 0.5d0*NORM_U*drho(:) + rho*(u*du(:)+v*dv(:)+w*dw(:))
    dpre_B(:) = Bx*dBx(:) + By*dBy(:) + Bz*dBz(:)

    de(:)  = (nablaQ(5,:) - dKE(:) - dpre_B(:) - drho(:)*inte)/rho

    rhox = drho(1)
    rhoy = drho(2)

    ux = du(1)
    uy = du(2)
    vx = dv(1)
    vy = dv(2)
    wx = dw(1)
    wy = dw(2)
    ex = de(1)
    ey = de(2)

    Bxy = dBx(2)
    Byx = dBy(1)
    Bzx = dBz(1)
    Bzy = dBz(2)

    ! Viscous term from hydrodynamics
    Fv(1) = 0.d0
    Fv(2) = mu_art*(2.d0*ux+lambda*(ux+vy))
    Fv(3) = mu_art*(vx+uy)
    Fv(4) = mu_art*wx
    Fv(5) = u*Fv(2)+v*Fv(3)+w*Fv(4)+mu_art*gam*ex/prandt

    Gv(1) = 0.d0
    Gv(2) = mu_art*(vx+uy)
    Gv(3) = mu_art*(2.d0*vy+lambda*(ux+vy))
    Gv(4) = mu_art*wy
    Gv(5) = Gv(2)*u+Gv(3)*v+Gv(4)*w+mu_art*gam*ey/prandt

    ! Resistive term from magnetic
    Jx = Bzy   - 0.0d0
    Jy = 0.0d0 - Bzx
    Jz = Byx   - Bxy

    Fv(6) = 0.0d0
    Fv(7) = eta_art*Jz
    Fv(8) = -eta_art*Jy
    Fv(5) = Fv(5) + Fv(7)*By + Fv(8)*Bz

    Gv(6) = -eta_art*Jz
    Gv(7) = 0.0d0
    Gv(8) = eta_art*Jx
    Gv(5) = Gv(5) + Gv(6)*Bx + Gv(8)*Bz

    Fv(numv) = 0.d0
    Gv(numv) = 0.d0

    END SUBROUTINE getVfluxvectors_AV

max_stars_repo_name: chenkuangxu/CHORUS-MHD

max_issues_repo_name: chenkuangxu/CHORUS-MHD

max_forks_repo_name: chenkuangxu/CHORUS-MHD


text:
      subroutine kostia(
     i                  rwdth)
c
c     + + + PURPOSE + + +
c     This subprogram calibrates modified Kostiakov infiltration
c     equation constants for a series of time-infiltration rate pairs
c     where infiltration rate is calculated using the Green and Ampt
c     infiltration function.
c
c     Called from FURROW
c     Author(s): E. R. Kottwitz
c     Reference in User Guide: Chapter 12
c
c     Changes:
c
c     Version: This module recoded from WEPP Version 93.06.
c     Date recoded: 07/07/93.
c     Recoded by: Charles R. Meyer.
c
c     + + + PARAMETERS + + +
      include 'pmxcrp.inc'
      include 'pmxhil.inc'
      include 'pmxnsl.inc'
      include 'pmxpln.inc'
      include 'pmxpnd.inc'
      include 'pmxres.inc'
      include 'pmxsrg.inc'
      include 'pmxtil.inc'
      include 'pmxtls.inc'
      include 'pntype.inc'
      include 'pxstep.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      real      rwdth
c
c     + + + ARGUMENT DEFINITIONS + + +
c     rwdth  - assumed row width of current overland flow element (m)
c
c     + + + COMMON BLOCKS + + +
      include 'cavepar.inc'
c     modify: aveks(mxplan),avesm(mxplan)
c
      include 'ccrpout.inc'
c       read: rh(mxplan)
c
      include 'cdist2.inc'
c       read: slplen(mxplan)
c
      include 'ciraflo.inc'
c       read: aqcnst,aqexp,botwid,sidslp
c
      include 'cirfur2.inc'
c       read: dtheta(mxplan),qspply(mxplan),surge,timtot
c     modify: srg
c
      include 'cirinfl.inc'
c     modify: kostf
c      write: kosta,kostk
c
      include 'cirriga.inc'
c       read: irofe,noirr
c
      include 'cparame.inc'
c       read: ks(mxplan),sm(mxplan)
c
      include 'cstruc.inc'
c       read: iplane
c
c     + + + LOCAL VARIABLES + + +
      integer   furint,i,nsteps
      real      avgdp2,deltai,denom,diydt,dlower,dupper,horcon,hormax,
     1          infadj,infil2,infilh(50),infilt,infilv(0:50),infold,
     1          inftmx,irsum1,irsum2,irsum3,irsum4,irsum5,kostfd,natlni,
     1          natlnt,piovr2,qmaxl,qmaxu,testf,time(50),timhmx,vertmx
      real      tmpvr1,tmpvr2,tmpvr3,tmpvr4
c
c     + + + LOCAL DEFINITIONS + + +
c     furint - irrigated furrow interval (1: all furrows irrigated, 2:
c              every other furrow irrigated)
c     i      - do loop counter
c     nsteps - number of infiltration increments
c     avgdp2 - two times the average of dlower and dupper (m)
c     deltai - 1-D vertical distance to the wetting front increment (m)
c     denom  - coefficient used as a denominator in calculating kosta
c              kostk (m**4)
c     diydt  - time derivative if vertical advance distance (m/s)
c     dlower - flow depth at lower end of current OFE for maximum of
c              furrow supply rates assuming steady state infiltration
c              (m)
c     dupper - flow depth at upper end of current OFE for maximum of
c              furrow supply rates (m)
c     horcon - horizontal infiltration constant equal to
c              (2*Ke*Ns)**.5/(change in soil moisture content) (m/s**.5)
c     hormax - maximum horizontal distance to the wetting front (m)
c     infadj - 2-D cumulative infiltration minus Ks*(infiltration
c              opportunity time)*(row spacing) (m**2)
c     infil2 - infilt*infilt (m**4)
c     infilh - 1-D horizontal distance to the wetting front (m)
c     infilt - 2-D cumulative infiltration (volume of water infiltrated
c              per unit furrow length) (m**2)
c     infilv - 1-D vertical distance to the wetting front as determined
c              by SR NEWTON (m)
c     infold - estimate of infiltration sent to SR NEWTON (m)
c     inftmx - cumulative 2-D infiltration at t=timhmx (m**2)
c     irsum1 - summation of infil2 (m**4)
c     irsum2 - summation of infil2*natlnt (m**4)
c     irsum3 - summation of infil2*natlnt**2 (m**4)
c     irsum4 - summation of infil2*natlni (m**4)
c     irsum5 - summation of infil2*natlnt*natlni (m**4)
c     kostfd - estimate of kostf using derivative of 2-D cumulative
c              infiltration function (m**2/s)
c     natlni - natural logarithm of infilt
c     natlnt - natural logarithm of time
c     piovr2 - one-half the ratio of circle area to square of its radius
c     qmaxl  - maximum flow rate at the lower end of the OFE (m**3/s)
c     qmaxu  - maximum furrow supply rate (m**3/s)
c     testf  - estimate of new value of kostf if old value was too large
c              (m**2/s)
c     time   - infiltration opportunity time (s)
c     timhmx - time when horizontal wetting fronts meet (s)
c     vertmx - vertical distance to the wetting front at t=timhmx (m)
c     tmpvr1 - variable to reduce number of calculations (m/s)
c     tmpvr2 - variable to reduce number of calculations (m)
c     tmpvr3 - variable to reduce number of calculations (m)
c
c     + + + SUBROUTINES CALLED + + +
c     newton
c
c     + + + DATA INITIALIZATIONS + + +
      data      nsteps/50/,piovr2/1.57080/
c
c     + + + OUTPUT FORMATS + + +
 2001 format(' *** WARNING ***',/,
     1       'Furrow irrigation flow depth (',f5.1,' mm) exceeds',
     1       ' ridge height (',f5.1,' mm).',/,
     1       ' *** WARNING ***',/)
c
c     + + + END SPECIFICATIONS + + +
c
c     Assume furint = 1.  User will eventually be required to provide
c     this information.
c
      furint = 1
c
c     Initialize variables
c
      infilv(0)=0.0
      inftmx = 0.0
      vertmx = 0.0
      kostfd = 1.0
      aveks(iplane) = ks(iplane)
      avesm(iplane) = sm(iplane)
      tmpvr1 = aveks(iplane)/dtheta(iplane)
      tmpvr2 = avesm(iplane)/dtheta(iplane)
      tmpvr3 = dtheta(iplane)*float(furint)*rwdth
c
c     Determine maximum furrow inflow rate
c
      qmaxu = qspply(1)
      if(noirr.eq.2 .and. surge.ge.2)then
        do 100  srg = 2, surge
          qmaxu = max(qmaxu,qspply(srg))
  100   continue
      endif
c
c     Determine flow depth at the upper end of the OFE
c
      if(sidslp.gt.0.001)then
        dupper = (sqrt(botwid**2+4.0*sidslp*(qmaxu/aqcnst)**(1.0/aqexp))
     1          -botwid)/2.0/sidslp
      else
        dupper = (qmaxu/aqcnst)**(1.0/aqexp)/botwid
      endif
c
c     Print warning to user if dupper exceeds ridge height
c
      if(iplane.eq.irofe .and. dupper.gt.rh(iplane))
     1       write (6,2001)dupper*1000.0,rh(iplane)*1000.0
c
c     Determine horizontal infiltration constant
c
c      ---- The value in parenthesis is from that portion of WEPP
c           Equation 12.3 also in parenthesis.  HORCON is the square
c           root of it.
c           Note that avesm(iplane) = (total head loss)*dtheta(iplane).
c
      horcon = sqrt(2.0*tmpvr1*tmpvr2)
c
c     Determine maximum horizontal wetting distance and time when
c     horizontal wetting fronts meet
c
c      ---- (WEPP Equation 12.5, with W = float(furint)*rwdth)
      hormax = (float(furint)*rwdth-botwid)*0.5
      timhmx = (hormax/horcon)**2
c
c     Determine vertical advance distance at t=timhmx
c
      infold = timhmx*aveks(iplane)
      call newton(timhmx,infold,aveks(iplane),avesm(iplane),vertmx)
      vertmx = vertmx/dtheta(iplane)
c
c     Determine 1-D vert. distance to the wetting front for t=timtot
c
      infold = timtot*aveks(iplane)
      call newton(timtot,infold,aveks(iplane),avesm(iplane),
     1            infilv(nsteps))
      infilv(nsteps) = infilv(nsteps)/dtheta(iplane)
c
c     Determine deltai
c
      deltai = infilv(nsteps)/float(nsteps)
c
c     Determine 1-D horizontal distance to the wetting front for
c     t=timtot
c
      if(timtot.lt.timhmx)then
c
c       ---- (WEPP Equation 12.3)
        infilh(nsteps) = horcon*sqrt(timtot)
      else
        infilh(nsteps) = hormax
      endif
c
c     Determine infiltration rate at t=timtot using 2-D infiltration
c     function
c
c      ---- (WEPP Equation 4.2.3) divided by change in soil water
c           content
      diydt = tmpvr1*(1.0+tmpvr2/infilv(nsteps))
c
c     If the following is true then horizontal wetting fronts have not
c     met.
c
      if(timtot.le.timhmx)then
c
c        ---- Time derivative of WEPP Equation 12.1
        kostfd = (infilh(nsteps)*(dupper/timtot+piovr2*infilv(nsteps)/
     1           (2.0*timtot)+piovr2*diydt)+botwid*diydt)*dtheta(iplane)
      else
        kostfd = diydt*tmpvr3
      endif
c
c     Calculate values for elements of infilh, infilv, and time arrays.
c     Also determine the minimum infiltration rate.
c
      do 400  i = 1, nsteps-1
c
c       Determine 1-D vertical distance to the wetting front
c
        infilv(i) = infilv(i-1)+deltai
c
c       Calculate time required for infiltration depth infilv(i)
c
c        ---- (WEPP Equation 12.2, solved for time (tau in 12.2))
        time(i) = infilv(i)/tmpvr1-avesm(iplane)/aveks(iplane)*
     1            log(1.0+infilv(i)/tmpvr2)
c
c       Determine 1-D horizontal advance distance to the wetting front
c
        if(time(i).lt.timhmx)then
c
c         ---- (WEPP Equation 12.3)
          infilh(i) = horcon*sqrt(time(i))
        else
          infilh(i) = hormax
        endif
c
c       Determine kostf using derivative of 2-D cumulative infiltration
c       function
c
c        ---- (WEPP Equation 4.2.3) divided by change in soil water
c             content
        diydt = tmpvr1*(1.0+tmpvr2/infilv(i))
        if(timtot.le.timhmx)then
c
c          ---- Time derivative of WEPP Equation 12.1
          kostf = (infilh(i)*(dupper/time(i)+piovr2*infilv(i)/(2.0*
     1            time(i))+piovr2*diydt)+botwid*diydt)*dtheta(iplane)
        else
          kostf = diydt*tmpvr3
        endif
        kostfd = min(kostf,kostfd)
  400 continue
c
c     Determine kostf using wetted perimeter at the upper end of the
c     furrow
c
      kostf = aveks(iplane)*(botwid+2.0*dupper*sqrt(1.0+sidslp**2))
c
c     Set kostf to minimum of possible values
c
      kostf = min(kostf,kostfd)
c
c     Determine maximum flow rate at the lower end of the OFE
c
   10 qmaxl = qmaxu-slplen(iplane)*kostf
c
c     Determine 2*(average depth of flow) for steady state conditions
c
      if(qmaxl.gt.0.0)then
        if(sidslp.gt.0.001)then
          dlower = (sqrt(botwid**2+4.0*sidslp*(qmaxl/aqcnst)**(1.0/
     1             aqexp))-botwid)/2.0/sidslp
        else
          dlower = (qmaxl/aqcnst)**(1.0/aqexp)/botwid
        endif
        avgdp2 = dupper+dlower
      else
        avgdp2 = dupper
      endif
c
c     Determine 2-D cumulative infiltration for t=timtot
c
      if(timtot.le.timhmx)then
c
c        ---- (WEPP Equation 12.1)
        infilt = (avgdp2*infilh(nsteps)+botwid*infilv(nsteps)+
     1           piovr2*infilh(nsteps)*infilv(nsteps))*dtheta(iplane)
      else
        inftmx = (avgdp2*hormax+botwid*vertmx+piovr2*hormax*vertmx)*
     1           dtheta(iplane)
        infilt = inftmx+(infilv(nsteps)-vertmx)*float(furint)*rwdth*
     1           dtheta(iplane)
      endif
c
c     Adjust 2-D cumulative infiltration to eliminate steady state
c     infiltration rate
c
      infadj = infilt-kostf*timtot
c
c     Adjust kostf if infadj is less than 0.000001 m**2
c
      if(infadj.lt.0.000001)then
        infadj = 0.000001
        testf = (infilt-infadj)/timtot
        if(testf.lt.0.999*kostf)then
          kostf = testf
        else
          kostf = 0.999*kostf
        endif
        goto 10
      endif
c
c     Determine values for functions of cumulative infiltration and time
c
      natlni = log(infadj)
      infil2 = infadj**2
      natlnt = log(timtot)
c
c     Calculate summations
c
      irsum1 = infil2
      irsum2 = infil2*natlnt
      irsum3 = irsum2*natlnt
      irsum4 = infil2*natlni
      irsum5 = irsum2*natlni
c
c     Initialize variables for calculation loop
c
      i = 0
c
c     Calculation loop
c
  200 i = i+1
c
c       Determine 2-D cumulative infiltration
c
        if(time(i).le.timhmx)then
c
c          ---- (WEPP Equation 12.1)
          infilt = (avgdp2*infilh(i)+botwid*infilv(i)+piovr2*
     1             infilh(i)*infilv(i))*dtheta(iplane)
        else
          infilt = inftmx+(infilv(i)-vertmx)*float(furint)*rwdth*
     1             dtheta(iplane)
        endif
c
c       Adjust 2-D cumulative infiltration to eliminate steady state
c       infiltration rate
c
        infadj = infilt-kostf*time(i)
c
c       Adjust kostf if infadj is less than 0.0001 m**2
c
        if(infadj.lt.0.000001)then
          testf = (infilt-infadj)/time(i)
          if(testf.lt.0.999*kostf)then
            kostf = testf
          else
            kostf = 0.999*kostf
          endif
          goto 10
        endif
c
c       Determine values for functions of infadj and time
c
        natlni = log(infadj)
        infil2 = infadj*infadj
        natlnt = log(time(i))
        tmpvr4 = infil2*natlnt
c
c       Perform summations
c
        irsum1 = irsum1+infil2
        irsum2 = irsum2+tmpvr4
        irsum3 = irsum3+tmpvr4*natlnt
        irsum4 = irsum4+infil2*natlni
        irsum5 = irsum5+tmpvr4*natlni
c
c     Check for completion of calculation loop
c
      if(i.lt.nsteps-1)goto 200
c
c     Determine Kostiakov-Lewis infiltration equation constants
c
      denom = irsum2*irsum2-irsum1*irsum3
      kosta = (irsum2*irsum4-irsum1*irsum5)/denom
      if(kosta.gt.1.0)then
        kosta = 1.0
        kostk = infadj/timtot
      else
        kostk = exp((irsum2*irsum5-irsum3*irsum4)/denom)
      endif
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
great(c1,v1).
great(jj1,aa1).
great(jj1,n1).
great(cc1,n1).
great(u1,ll1).
great(bb1,t1).
great(f1,l1).
great(jj1,o1).
great(d1,y1).
great(j1,ee1).
great(bb1,p1).
great(c1,d1).
great(hh1,jj1).
great(i1,t1).
great(hh1,h1).
great(e1,t1).
great(k1,s1).
great(g1,bb1).
great(f1,hh1).
great(j1,t1).
great(c1,dd1).
great(f1,ll1).
great(v1,x1).
great(d1,o1).
great(d1,aa1).
great(e1,n1).
great(f1,q1).
great(i1,ll1).
great(d1,kk1).
great(u1,ee1).
great(f1,r1).
great(r1,o1).
great(a1,t1).
great(k1,v1).
great(l1,u1).
great(c1,ll1).
great(c1,x1).
great(y1,i1).
great(y1,ee1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp