subs1.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
subroutine setdt
  use kinds, only : dp
  use size_m, only : nid
  use input, only : param, ifflow, ifprint
  use soln, only : vx, vy, vz
  use tstep, only : dt, dtinit, time, fintim, lastep, courno, ctarg, istep
  use tstep, only : re_cell
  implicit none

  real(DP), save :: umax = 0._dp
  real(DP), save :: uxmax = 0._dp
  REAL(DP), save :: DTOLD = 0._dp
  REAL(DP), save :: DTOpf = 0._dp
  logical, save :: iffxdt = .FALSE.

  real(DP) :: dtcfl, dtfs, dtmin, dtmax
  real(DP), external :: uglmin

  if (param(12) < 0 .OR. iffxdt) then
      iffxdt    = .true. 
      param(12) = abs(param(12))
      dt        = param(12)
      dtopf     = dt
      call compute_cfl(umax,vx,vy,vz,1.0)
      goto 200
  else IF (param(84) /= 0.0) THEN
      if (istep < 6) then
          dt   =param(84)
          dtold=param(84)
          dtopf=param(84)
          return
      else
          dtold=dt
          dtopf=dt
          dt=dtopf*param(85)
          dt=min(dt,param(12))
      endif
  endif

!   Find DT=DTCFL based on CFL-condition (if applicable)

  CALL setdtc(umax, uxmax)
  re_cell = uxmax / param(2)
  dtcfl = dt

!   Find DTFS based on surface tension (if applicable)

!  CALL SETDTFS (DTFS)

!   Select appropriate DT
  dtfs = 0.
  IF ((dt == 0.) .AND. (dtfs > 0.)) THEN
      dt = dtfs
  ELSEIF ((dt > 0.) .AND. (dtfs > 0.)) THEN
      dt = min(dt,dtfs)
  ELSEIF ((dt == 0.) .AND. (dtfs == 0.)) THEN
      dt = 0.
      IF (ifflow .AND. nid == 0 .AND. ifprint) THEN
          WRITE (6,*) 'WARNING: CFL-condition & surface tension'
          WRITE (6,*) '         are not applicable'
      endif
  ELSEIF ((dt > 0.) .AND. (dtfs == 0.)) THEN
      dt = dt
  ELSE
      dt = 0.
      IF (nid == 0) WRITE (6,*) 'WARNING: DT<0 or DTFS<0'
      IF (nid == 0) WRITE (6,*) '         Reset DT      '
  endif

!   Check DT against user-specified input, DTINIT=PARAM(12).

  IF ((dt > 0.) .AND. (dtinit > 0.)) THEN
      dt = min(dt,dtinit)
  ELSEIF ((dt == 0.) .AND. (dtinit > 0.)) THEN
      dt = dtinit
  ELSEIF ((dt > 0.) .AND. (dtinit == 0.)) THEN
      dt = dt
  ELSEIF ( .NOT. iffxdt) THEN
      dt = 0.001
      IF(nid == 0)WRITE (6,*) 'WARNING: Set DT=0.001 (arbitrarily)'
  endif

!   Check if final time (user specified) has been reached.

  200 IF (time+dt >= fintim .AND. fintim /= 0.0) THEN
  !        Last step
      lastep = 1
      dt = fintim-time
      IF (nid == 0) WRITE (6,*) 'Final time step = ',dt
  endif

  courno = dt*umax
  IF (nid == 0 .AND. ifprint .AND. dt /= dtold) &
  WRITE (6,100) dt,dtcfl,dtfs,dtinit
  100 FORMAT(5x,'DT/DTCFL/DTFS/DTINIT',4e12.3)

  ! only change dt if its off by more than 10%
  if (abs(dt - dtold)/dt < 0.1) dt = dtold

  ! Put limits on how much DT can change.
  IF (dtold /= 0.0) THEN
      dtmin=0.8*dtold
      dtmax=1.2*dtold
      dt = min(dtmax,dt)
      dt = max(dtmin,dt)
  endif
  dtold=dt

!    IF (PARAM(84).NE.0.0) THEN
!          dt=dtopf*param(85)
!          dt=min(dt,param(12))
!    endif

  if (iffxdt) dt=dtopf
  courno = dt*umax

! synchronize time step for multiple sessions
!max  if (ifneknek) dt=uglmin(dt,1)

  if (iffxdt .AND. abs(courno) > 10.*abs(ctarg)) then
      if (nid == 0) write(6,*) 'CFL, Ctarg!',courno,ctarg
!max      call emerxit
  endif

  return
end subroutine setdt


!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine cvgnlps (ifconv)
  use kinds, only : dp
  use input, only : ifnonl
  use tstep, only : ifield, tnrmh1, tolnl
  implicit none

  LOGICAL ::  IFCONV
  real(DP) :: tnorm1, tnorm2, eps

  IF (ifnonl(ifield)) THEN
      ifconv = .false. 
  ELSE
      ifconv = .true. 
      return
  endif

  tnorm1 = tnrmh1(ifield-1)
  CALL unorm
  tnorm2 = tnrmh1(ifield-1)
  eps = abs((tnorm2-tnorm1)/tnorm2)
  IF (eps < tolnl) ifconv = .true. 

  return
end subroutine cvgnlps

!---------------------------------------------------------------------
!---------------------------------------------------------------------
subroutine unorm
  use kinds, only : dp
  use soln, only : vx, vy, vz, t
  use tstep, only : ifield, imesh, tnrml8, tnrmsm, tnrmh1, time, dt, tnrml2
  use tstep, only : vnrmh1, vnrmsm, vnrml2, vnrml8, istep, dtinvm, vmean, tmean
  implicit none

  real(DP) :: tden, arg

  IF (ifield == 1) THEN
  
  !        Compute norms of the velocity.
  !        Compute time mean (L2) of the inverse of the time step.
  !        Compute L2 in time, H1 in space of the velocity.
  
      CALL normvc(vnrmh1,vnrmsm,vnrml2,vnrml8,vx,vy,vz)
      IF (istep == 0) return
      IF (istep == 1) THEN
          dtinvm = 1./dt
          vmean  = vnrml8
      ELSE
          tden   = time
          if (time <= 0) tden = abs(time)+1.e-9
          arg    = ((time-dt)*dtinvm**2+1./dt)/tden
          if (arg > 0) dtinvm = sqrt(arg)
          arg    = ((time-dt)*vmean**2+dt*vnrmh1**2)/tden
          if (arg > 0) vmean  = sqrt(arg)
      endif
  ELSE
  
  !     Compute norms of a passive scalar
  
      CALL normsc(tnrmh1(ifield-1),tnrmsm(ifield-1), &
      tnrml2(ifield-1),tnrml8(ifield-1), &
      t(1,1,1,1,ifield-1),imesh)
      tmean(ifield-1) = 0.
  endif

  return
end subroutine unorm

!--------------------------------------------------------------
!--------------------------------------------------------------
subroutine setdtc(umax, uxmax)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelv
  use size_m, only : nx1, ny1, nz1, nelv, ndim, nid
  use geom, only : ifwcno, xm1, ym1, zm1
  use input, only : param, iftran, ifflow, ifnav, ifheat, ifadvc
  use input, only : ipscal, npscal
  use geom, only : binvm1
  use soln, only : vx, vy, vz, v1mask, v2mask, v3mask, bfx, bfy, bfz
  use tstep, only : lastep, dt, ifield, courno, ctarg, avtran, dtinit
  implicit none

  real(DP), intent(out) :: umax, uxmax
  real(DP), allocatable :: u(:,:,:,:), v(:,:,:,:), w(:,:,:,:)

  REAL(DP), save :: VCOUR

  INTEGER, save :: IFIRST = 0
  integer :: irst, iconv, ntot, ntotl, ntotd, i
  real(DP) :: dtold, cold, cmax, cmin, fmax, density, amax, dxchar, vold, cpred
  real(DP) :: a, b, c, discr, dtlow, dthi
  real(DP), external :: glmax, glmin

!   Steady state => all done
  IF ( .NOT. iftran) THEN
      ifirst=1
      lastep=1
      return
  endif

  irst = int(param(46))
  if (irst > 0) ifirst=1

!   First time around

  IF (ifirst == 0) THEN
      dt     = dtinit
      IF (ifflow) THEN
          ifield = 1
          CALL bcdirvc(vx,vy,vz,v1mask,v2mask,v3mask)
      endif
  endif
  ifirst=ifirst+1

  dtold = dt

!   Convection ?

!   Don't enforce Courant condition if there is no convection.


  iconv=0
  IF (ifflow .AND. ifnav) iconv=1
  IF (ifwcno)             iconv=1
  IF (ifheat) THEN
      DO 10 ipscal=0,npscal
          IF (ifadvc(ipscal+2)) iconv=1
      10 END DO
  endif
  IF (iconv == 0) THEN
      dt=0.
      return
  endif


!   Find Courant and Umax


  ntot   = nx1*ny1*nz1*nelv
  ntotl  = lx1*ly1*lz1*lelv
  ntotd  = ntotl*ndim
  cold   = courno
  cmax   = 1.2*ctarg
  cmin   = 0.8*ctarg
  allocate(u(lx1,ly1,lz1,lelv), v(lx1,ly1,lz1,lelv), w(lx1,ly1,lz1,lelv))
  CALL cumax(vx,vy,vz,u,v,w,umax, uxmax)

!   Zero DT

  IF (dt == 0.0) THEN
  
      IF (umax /= 0.0) THEN
          dt = ctarg/umax
          vcour = umax
      ELSEIF (ifflow) THEN
      
      !           We'll use the body force to predict max velocity
      
          CALL setprop
          ifield = 1
      
          CALL makeuf
          CALL opdssum(bfx,bfy,bfz)
          bfx = bfx * binvm1; bfy = bfy * binvm1; bfz = bfz * binvm1
          fmax=0.0
          u = 0._dp
          DO 600 i=1,ntot
              u(i,1,1,1) = abs(bfx(i,1,1,1))
              v(i,1,1,1) = abs(bfy(i,1,1,1))
              w(i,1,1,1) = abs(bfz(i,1,1,1))
          600 END DO
          fmax    = glmax(u,ntotd)
          density = avtran(1)
          amax    = fmax/density
          dxchar  = sqrt( (xm1(1,1,1,1)-xm1(2,1,1,1))**2 + &
          (ym1(1,1,1,1)-ym1(2,1,1,1))**2 + &
          (zm1(1,1,1,1)-zm1(2,1,1,1))**2 )
          dxchar  = glmin((/dxchar/),1)
          IF (amax /= 0.) THEN
              dt = sqrt(ctarg*dxchar/amax)
          ELSE
              IF (nid == 0) &
              WRITE (6,*) 'CFL: Zero velocity and body force'
              dt = 0.0
              return
          endif
      ELSEIF (ifwcno) THEN
          IF (nid == 0) &
          WRITE (6,*) ' Stefan problem with no fluid flow'
          dt = 0.0
          return
      endif
  
  ELSEIF ((dt > 0.0) .AND. (umax /= 0.0)) THEN
  
  
  !     Nonzero DT & nonzero velocity
  
  
      courno = dt*umax
      vold   = vcour
      vcour  = umax
      IF (ifirst == 1) THEN
          cold = courno
          vold = vcour
      endif
      cpred  = 2.*courno-cold
  
  !     Change DT if it is too big or if it is too small
  
  !     if (nid.eq.0)
  !    $write(6,917) dt,umax,vold,vcour,cpred,cmax,courno,cmin
  ! 917 format(' dt',4f9.5,4f10.6)
      IF(courno > cmax .OR. cpred > cmax .OR. courno < cmin) THEN
      
          a=(vcour-vold)/dt
          b=vcour
      !           -C IS Target Courant number
          c=-ctarg
          discr=b**2-4*a*c
          dtold=dt
          IF(discr <= 0.0)THEN
              if (nid == 0) &
              print*,'Problem calculating new DT Discriminant=',discr
              dt=dt*(ctarg/courno)
          !               IF(DT.GT.DTOLD) DT=DTOLD
          ELSE IF(abs((vcour-vold)/vcour) < 0.001)THEN
          !              Easy: same v as before (LINEARIZED)
              dt=dt*(ctarg/courno)
          !     if (nid.eq.0)
          !    $write(6,918) dt,dthi,dtlow,discr,a,b,c
          ! 918 format(' d2',4f9.5,4f10.6)
          ELSE
              dtlow=(-b+sqrt(discr) )/(2.0*a)
              dthi =(-b-sqrt(discr) )/(2.0*a)
              IF(dthi > 0.0 .AND. dtlow > 0.0)THEN
                  dt = min(dthi,dtlow)
              !     if (nid.eq.0)
              !    $write(6,919) dt,dthi,dtlow,discr,a,b,c
              ! 919 format(' d3',4f9.5,4f10.6)
              ELSE IF(dthi <= 0.0 .AND. dtlow <= 0.0)THEN
              !                 PRINT*,'DTLOW,DTHI',DTLOW,DTHI
              !                 PRINT*,'WARNING: Abnormal DT from CFL-condition'
              !                 PRINT*,'         Keep going'
                  dt=dt*(ctarg/courno)
              ELSE
              !                 Normal case; 1 positive root, one negative root
                  dt = max(dthi,dtlow)
              !     if (nid.eq.0)
              !    $write(6,929) dt,dthi,dtlow,discr,a,b,c
              ! 929 format(' d4',4f9.5,4f10.6)
              endif
          endif
      !           We'll increase gradually-- make it the geometric mean between
      !     if (nid.eq.0)
      !    $write(6,939) dt,dtold
      ! 939 format(' d5',4f9.5,4f10.6)
          IF (dtold/dt < 0.2) dt = dtold*5
      endif
  
  endif

  return
end subroutine setdtc

subroutine cumax (v1,v2,v3,u,v,w,umax, uxmax)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelv, nx1, ny1, nz1, nelv, ndim
  use geom, only : rxm1, rym1, sxm1, sym1, rzm1, szm1, txm1, tym1, tzm1
  use input, only : ifaxis
  use wz_m, only : zgm1
  use mesh, only : if_ortho
  implicit none

  real(DP), intent(in)  :: V1(lx1,ly1,lz1,lelv)
  real(DP), intent(in)  :: V2(lx1,ly1,lz1,lelv)
  real(DP), intent(in)  :: V3(lx1,ly1,lz1,lelv)
  real(DP), intent(out) :: u(lx1,ly1,lz1,lelv) ! scratch
  real(DP), intent(out) :: v(lx1,ly1,lz1,lelv) ! scratch
  real(DP), intent(out) :: w(lx1,ly1,lz1,lelv) ! scratch
  real(DP), intent(out) :: umax
  real(DP), intent(out) :: uxmax

  real(DP), allocatable, dimension(:,:,:,:) :: &
    xrm1, xsm1, xtm1, yrm1, ysm1, ytm1, zrm1, zsm1, ztm1

  real(DP), allocatable :: x(:,:,:,:), r(:,:,:,:)
  real(DP), allocatable :: tmp(:,:,:,:)

  real(DP), save :: drst(lx1), drsti(lx1)

  real(DP) :: U3(3)
  INTEGER, save :: ICALLD = 0

  integer :: ntot, ntotl, ntotd
  integer :: i, ie, ix, iy, iz
  real(DP), external :: glmax


  ntot  = nx1*ny1*nz1*nelv
  ntotl = lx1*ly1*lz1*lelv
  ntotd = ntotl*ndim

!   Compute isoparametric partials.
  allocate(xrm1(nx1,ny1,nz1,nelv), ysm1(nx1,ny1,nz1,nelv), ztm1(nx1,ny1,nz1,nelv))
  if (.not. if_ortho) then
    allocate(xsm1(nx1,ny1,nz1,nelv), xtm1(nx1,ny1,nz1,nelv))
    allocate(yrm1(nx1,ny1,nz1,nelv), ytm1(nx1,ny1,nz1,nelv))
    allocate(zrm1(nx1,ny1,nz1,nelv), zsm1(nx1,ny1,nz1,nelv))
  else
    allocate(xsm1(nx1,ny1,nz1,1), xtm1(nx1,ny1,nz1,1))
    allocate(yrm1(nx1,ny1,nz1,1), ytm1(nx1,ny1,nz1,1))
    allocate(zrm1(nx1,ny1,nz1,1), zsm1(nx1,ny1,nz1,1))
  endif
  CALL xyzrst(xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,xtm1,ytm1,ztm1, ifaxis)

!   Compute maximum U/DX

  IF (icalld == 0) THEN
      icalld=1
      drst(1)=abs(zgm1(2,1)-zgm1(1,1))
      drsti(1)=1.0/drst(1)
      DO 400 i=2,nx1-1
          drst(i)=abs(zgm1(i+1,1)-zgm1(i-1,1))/2.0
          drsti(i)=1.0/drst(i)
      400 END DO
      drst(nx1)=drst(1)
      drsti(nx1)=1.0/drst(nx1)
  endif

!   Zero out scratch arrays U,V,W for ALL declared elements...

  u = 0_dp; v = 0_dp; w = 0_dp

  allocate(x(lx1,ly1,lz1,lelv), r(lx1,ly1,lz1,lelv))

  IF (ndim == 2) THEN
#if 0
      CALL vdot2(u,v1  ,v2  ,rxm1,rym1,ntot)
      CALL vdot2(r,rxm1,rym1,rxm1,rym1,ntot)
      CALL vdot2(x,xrm1,yrm1,xrm1,yrm1,ntot)
      r = r * x
      r = sqrt(r) 
      u = u / r
  
      CALL vdot2(v,v1  ,v2  ,sxm1,sym1,ntot)
      CALL vdot2(r,sxm1,sym1,sxm1,sym1,ntot)
      CALL vdot2(x,xsm1,ysm1,xsm1,ysm1,ntot)
      r = r * x
      r = sqrt(r) 
      v = v / r
#endif 
  ELSE

      if (if_ortho) then 
        u = v1 * rxm1 
        r = rxm1 * rxm1 
        x = xrm1 * xrm1 
      else
        u = v1 * rxm1 + v2 * rym1 + v3 * rzm1 
        r = rxm1 * rxm1 + rym1 * rym1 + rzm1 * rzm1
        x = xrm1 * xrm1 + yrm1 * yrm1 + zrm1 * zrm1
      endif

      r = r * x
      r = sqrt(r) 
      u = u / r

      if (if_ortho) then 
        v = v2*sym1
        r = sym1 * sym1
        x = ysm1 * ysm1
      else
        v = v1*sxm1     + v2*sym1     + v3*szm1
        r = sxm1 * sxm1 + sym1 * sym1 + szm1 * szm1
        x = xsm1 * xsm1 + ysm1 * ysm1 + zsm1 * zsm1 
      endif

      r = r * x
      r = sqrt(r) 
      v = v / r

      if (if_ortho) then
        w = v3*tzm1
        r = tzm1 * tzm1
        x = ztm1 * ztm1
      else
        w = v1*txm1     + v2*tym1     + v3*tzm1
        r = tzm1 * tzm1
        x = ztm1 * ztm1
      endif
 
      r = r * x
      r = sqrt(r) 
      w = w / r
  
  endif

  deallocate(xsm1, xtm1, yrm1, ytm1, zrm1, zsm1)
  deallocate(x,r)

  allocate(tmp(lx1,ly1,lz1,lelv))

  DO ie=1,nelv
      DO ix=1,nx1
          DO iy=1,ny1
              DO iz=1,nz1
                  tmp(ix,iy,iz,ie)=abs( v1(ix,iy,iz,ie)*xrm1(ix,iy,iz,ie)*drst(ix) )
              enddo
          enddo
      enddo
  END DO
  u3(1)   = maxval(tmp)

  DO ie=1,nelv
      DO ix=1,nx1
          DO iy=1,ny1
              DO iz=1,nz1
                  tmp(ix,iy,iz,ie)=abs( v2(ix,iy,iz,ie)*ysm1(ix,iy,iz,ie)*drst(iy) )
              enddo
          enddo
      enddo
  END DO
  u3(2)   = maxval(tmp)

  DO ie=1,nelv
      DO ix=1,nx1
          DO iy=1,ny1
              DO iz=1,nz1
                  tmp(ix,iy,iz,ie)=abs( v3(ix,iy,iz,ie)*ztm1(ix,iy,iz,ie)*drst(iz) )
              enddo
          enddo
      enddo
  END DO
  u3(3)   = maxval(tmp)
  uxmax    = glmax(u3,3)
  deallocate(xrm1, ysm1, ztm1)

  DO ie=1,nelv
      DO ix=1,nx1
          DO iy=1,ny1
              DO iz=1,nz1
                  u(ix,iy,iz,ie)=abs( u(ix,iy,iz,ie)*drsti(ix) )
                  v(ix,iy,iz,ie)=abs( v(ix,iy,iz,ie)*drsti(iy) )
                  w(ix,iy,iz,ie)=abs( w(ix,iy,iz,ie)*drsti(iz) )
              enddo
          enddo
      enddo
  END DO

  u3(1)   = maxval(u)
  u3(2)   = maxval(v)
  u3(3)   = maxval(w)
  umax    = glmax(u3,3)

  return
end subroutine cumax

subroutine faccl2(a,b,iface1)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, nx1, ny1, nz1
  use topol, only : eface1, skpdat
  implicit none

  real(DP) :: A(lx1,ly1,lz1),B(lx1,ly1)
  real(DP) :: af(lx1*ly1*lz1), bf(lx1*ly1)
  integer :: iface1

  integer :: j1, j2, i, jskip2, jf2, js2, jskip1, jf1, js1, iface

!   Set up counters

  CALL dsset(nx1,ny1,nz1)
  iface  = eface1(iface1)
  js1    = skpdat(1,iface)
  jf1    = skpdat(2,iface)
  jskip1 = skpdat(3,iface)
  js2    = skpdat(4,iface)
  jf2    = skpdat(5,iface)
  jskip2 = skpdat(6,iface)

  af = reshape(a, (/lx1*ly1*lz1/))
  bf = reshape(b, (/lx1*ly1/))
  i = 0
  DO j2=js2,jf2,jskip2
      DO j1=js1,jf1,jskip1
          i = i+1
          af(j1+lx1*(j2-1)) = af(j1+lx1*(j2-1))*bf(i)
      enddo
  END DO
  a = reshape(af, (/lx1, ly1, lz1/))

  return
end subroutine faccl2

subroutine faccl3(a,b,c,iface1)
!
! Collocate B with A on the surface IFACE1 of element IE.
!
!     A is a (NX,NY,NZ) data structure
!     B is a (NX,NY,IFACE) data structure
!     IFACE1 is in the preprocessor notation 
!     IFACE  is the dssum notation.
!     5 Jan 1989 15:12:22      PFF

  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, nx1, ny1, nz1
  use topol, only : eface1, skpdat

  real(DP) :: A(lx1,ly1,lz1),B(lx1,ly1,lz1),C(lx1,ly1)
  integer :: iface1, iface, js1, jf1, jskip1, js2, jf2, jskip2
  integer :: i, j2, j1
  real(DP) :: af(lx1*ly1*lz1), bf(lx1*ly1*lz1), cf(lx1*ly1)
  af = reshape(a, (/lx1*ly1*lz1/))
  bf = reshape(b, (/lx1*ly1*lz1/))
  cf = reshape(c, (/lx1*ly1/))

!
!     Set up counters
!
  CALL dsset(nx1,ny1,nz1)
  iface  = eface1(iface1)
  js1    = skpdat(1,iface)
  jf1    = skpdat(2,iface)
  jskip1 = skpdat(3,iface)
  js2    = skpdat(4,iface)
  jf2    = skpdat(5,iface)
  jskip2 = skpdat(6,iface)

  i = 0
  DO j2=js2,jf2,jskip2
    DO j1=js1,jf1,jskip1
     i = i+1
     !A(J1,J2,1) = B(J1,J2,1)*C(I,1)
     af(j1+lx1*(j2-1)) = bf(j1+lx1*(j2-1))*cf(i)
    enddo
  enddo
  a = reshape(af, (/lx1, ly1, lz1/))

  return
end


!-----------------------------------------------------------------------
subroutine sethlm (h1,h2,intloc)
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1
  use size_m, only : lx1, ly1, lz1, lelt
  use input, only : iftran, param
  use soln, only : vdiff, vtrans
  use tstep, only : ifield, nelfld, bd, dt
  implicit none

  real(DP) :: h1(lx1,ly1,lz1,lelt),h2(lx1,ly1,lz1,lelt)
  integer :: intloc

  integer :: nel, ntot1
  real(DP) :: dtbd

  nel   = nelfld(ifield)
  ntot1 = nx1*ny1*nz1*nel

  if (iftran) then
      dtbd = bd(1)/dt
      h1 = vdiff(:,:,:,:,ifield)
      if (intloc == 0) then
          h2 = 0._dp
      else
          if (ifield == 1 .OR. param(107) == 0) then
              h2 = vtrans(:,:,:,:,ifield) * dtbd
          else   ! unsteady reaction-diffusion type equation
              h2 = vtrans(:,:,:,:,ifield) * dtbd + param(107)
          endif

      endif

  !        if (ifield.eq.1 .and. ifanls) then   ! this should be replaced
  !           const = 2.                        ! with a correct stress
  !           call cmult (h1,const,ntot1)       ! formulation
  !        endif

  ELSE
      CALL copy(h1,vdiff(1,1,1,1,ifield),ntot1)
      h2 = 0._dp
      if (param(107) /= 0) then
          write(6,*) 'SPECIAL SETHLM!!',param(107)
      !           call cfill (h2,param(107),ntot1)
          call copy(h2,vtrans(1,1,1,1,ifield),ntot1)
      endif
  endif

  return
end subroutine sethlm

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine vprops
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1
  use input, only : ifuservp, ifvarp, matype, igroup, iflomach
  use input, only : cpfld, cpgrp
  use soln, only : vdiff, vtrans
  use tstep, only : ifield, nelfld, istep
  implicit none

  integer :: nxyz1, nel, ntot1, iel, igrp, itype, itest
  real(DP) :: cdiff, ctrans
  integer, external :: iglmax

  nxyz1 = nx1*ny1*nz1
  nel   = nelfld(ifield)
  ntot1 = nxyz1*nel

  IF (istep == 0) THEN
  
  !        First time around, set defaults
  
      ifvarp(ifield) = .false. 
      if (iflomach) ifvarp(ifield) = .true. 

      if ( .NOT. ifvarp(ifield)) then ! check all groups
          do iel=1,nel
              igrp  = igroup(iel)
              itype = matype(igrp,ifield)
              if(itype /= 0) ifvarp(ifield) = .true. 
          enddo
      endif

      itest = 0                        ! test against all processors
      if (ifvarp(ifield)) itest = 1
      itest = iglmax((/itest/),1)
      if (itest > 0) ifvarp(ifield) = .true. 

  endif

!   Fill up property arrays every time step

!   First, check for turbulence models

#if 0
  IF (ifmodel .AND. ifkeps) THEN
      CALL turbfld(ifkfld,ifefld)
      IF (ifkfld)           CALL tpropk
      IF (ifefld)           CALL tprope
      IF (ifkfld .OR. ifefld) return
  endif
#endif

!...  No turbulence models, OR current field is not k or e.

  DO iel=1,nel
  
      igrp=igroup(iel)

      if (ifuservp) then
#if 0
      
      !           User specified fortran function   (pff 2/13/01)
          CALL nekuvp(iel)
          difmin = vlmin(vdiff(1,1,1,iel,ifield),nxyz1)
          IF (difmin <= 0.0) THEN
              WRITE (6,100) difmin,ifield,igrp
              CALL exitt
          endif
#endif        
      ELSE IF(matype(igrp,ifield) == 1)THEN
      
      !           Constant property within groups of elements
      
          cdiff  = cpgrp(igrp,ifield,1)
          ctrans = cpgrp(igrp,ifield,2)
          vdiff(:,:,:,iel,ifield) = cdiff
          vtrans(:,:,:,iel,ifield) = ctrans
          IF (cdiff <= 0.0) THEN
              WRITE(6,100) cdiff,ifield,igrp
              100 FORMAT(2x,'ERROR:  Non-positive diffusivity (' &
              ,g12.3,') specified for field',i2,', group',i2 &
              ,' element',i4,'.' &
              ,/,'ABORTING in VPROPS',//)
              CALL exitt
          endif
      
      ELSE IF(matype(igrp,ifield) == 2)THEN
        write(*,*) "Oops: matype" 
#if 0
      !           User specified fortran function
      
          CALL nekuvp(iel)
      
          difmin = vlmin(vdiff(1,1,1,iel,ifield),nxyz1)
          IF (difmin <= 0.0) THEN
              WRITE (6,100) difmin,ifield,igrp
              CALL exitt
          endif
#endif        
      ELSE IF(matype(igrp,ifield) == 0)THEN
      
      !           Default constant property
      
          cdiff  = cpfld(ifield,1)
          ctrans = cpfld(ifield,2)
      !           write(6,*) 'vdiff:',ifield,cdiff,ctrans
          vdiff(:,:,:,iel,ifield) = cdiff
          vtrans(:,:,:,iel,ifield) = ctrans

          IF (cdiff <= 0.0) THEN
              WRITE(6,200) cdiff,ifield
              200 FORMAT(2x,'ERROR:  Non-positive diffusivity (' &
              ,g12.3,') specified for field',i2,'.',/ &
              ,'ABORTING in VPROPS',//)
              CALL exitt
          endif
      endif
  
  END DO

!     Turbulence models --- sum eddy viscosity/diffusivity
#if 0
  IF (ifmodel .AND. (ifield == 1 .OR. ifield == 2)) &
  CALL tviscos
#endif

  return
end subroutine vprops

!-----------------------------------------------------------------------
subroutine setsolv
  use mesh, only : ifsolv
  implicit none
    
  ifsolv = .false. 
  return
end subroutine setsolv

!-----------------------------------------------------------------------
subroutine flush_io
  return
end subroutine flush_io

!-----------------------------------------------------------------------
!!     IOP = 0
!!     Extract vector (A1,A2,A3) from (B1,B2,B3) on face IFACE1.
!!     IOP = 1
!!     Extract vector (B1,B2,B3) from (A1,A2,A3) on face IFACE1.
!!     A1, A2, A3 have the (NX,NY,NFACE) data structure
!!     B1, B2, B3 have the (NX,NY,NZ)    data structure
!!     IFACE1 is in the preprocessor notation
!!     IFACE  is the dssum notation.
SUBROUTINE facexv (A1,A2,A3,B1,B2,B3,IFACE1,IOP)
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, lx1, ly1, lz1
  use topol, only : eface1, skpdat
  implicit none

  real(DP) :: A1(lx1,ly1),    A2(lx1,ly1),    A3(lx1,ly1)
  real(DP) :: B1(lx1,ly1,lz1),B2(lx1,ly1,lz1),B3(lx1,ly1,lz1)
  real(DP) :: A1f(lx1*ly1), A2f(lx1*ly1), A3f(lx1*ly1)
  real(DP) :: B1f(lx1*ly1*lz1), b2f(lx1*ly1*lz1), b3f(lx1*ly1*lz1)
  integer :: iface1, iop

  integer :: j1, j2, i, jskip2, jf2, js2, jskip1, jf1, js1, iface

  CALL dsset(nx1,ny1,nz1)
  iface  = eface1(iface1)
  js1    = skpdat(1,iface)
  jf1    = skpdat(2,iface)
  jskip1 = skpdat(3,iface)
  js2    = skpdat(4,iface)
  jf2    = skpdat(5,iface)
  jskip2 = skpdat(6,iface)
  i = 0

  IF (iop == 0) THEN
      b1f = reshape(b1, (/lx1*ly1*lz1/))
      b2f = reshape(b2, (/lx1*ly1*lz1/))
      b3f = reshape(b3, (/lx1*ly1*lz1/))
      DO j2=js2,jf2,jskip2
          DO j1=js1,jf1,jskip1
              i = i+1
              a1f(i) = b1f(j1+lx1*(j2-1))
              a2f(i) = b2f(j1+lx1*(j2-1))
              a3f(i) = b3f(j1+lx1*(j2-1))
          enddo
      END DO
      a1 = reshape(a1f, (/lx1, ly1/))
      a2 = reshape(a2f, (/lx1, ly1/))
      a3 = reshape(a3f, (/lx1, ly1/))
  ELSE
      a1f = reshape(a1, (/lx1*ly1/))
      a2f = reshape(a2, (/lx1*ly1/))
      a3f = reshape(a3, (/lx1*ly1/))
      DO j2=js2,jf2,jskip2
          DO j1=js1,jf1,jskip1
              i = i+1
              b1f(j1+lx1*(j2-1)) = a1f(i)
              b2f(j1+lx1*(j2-1)) = a2f(i)
              b3f(j1+lx1*(j2-1)) = a3f(i)
          enddo
      END DO
      b1 = reshape(b1f, (/lx1, ly1, lz1/))
      b2 = reshape(b2f, (/lx1, ly1, lz1/))
      b3 = reshape(b3f, (/lx1, ly1, lz1/))
  ENDIF

  RETURN
END SUBROUTINE facexv

SUBROUTINE updmsys (IFLD)
  use size_m
  use geom, only : iflmsf
  implicit none

  integer, intent(in) :: ifld

  IF ( .NOT. iflmsf(ifld)) RETURN
#if 0
  nel  = nelfld(ifld)
  CALL sethmsk(hvmask,hfmask,ifld,nel)
  CALL setcsys(hvmask,hfmask,nel)
#endif

  RETURN
END SUBROUTINE updmsys

SUBROUTINE setcdof
  use size_m, only : ndim, nelt
  use input, only : cbc, cdof
  implicit none 

  integer :: nface, ifc, iel
  nface = 2*ndim

  DO iel=1,nelt
      DO ifc=1,nface
          cdof(ifc,iel)=cbc(ifc,iel,0)(1:1)
      enddo
  END DO

  RETURN
END SUBROUTINE setcdof