xref: /phasta/phSolver/incompressible/itrdrv.f (revision bd36043d6c86bfee311b9eaff539555a870d94ea)

      subroutine itrdrv (y,         ac,
     &                   uold,      x,
     &                   iBC,       BC,
     &                   iper,      ilwork,     shp,
     &                   shgl,      shpb,       shglb,
     &                   ifath,     velbar,     nsons )
c
c----------------------------------------------------------------------
c
c This iterative driver is the semi-discrete, predictor multi-corrector
c algorithm. It contains the Hulbert Generalized Alpha method which
c is 2nd order accurate for Rho_inf from 0 to 1.  The method can be
c made  first-order accurate by setting Rho_inf=-1. It uses CGP and
c GMRES iterative solvers.
c
c working arrays:
c  y      (nshg,ndof)           : Y variables
c  x      (nshg,nsd)            : node coordinates
c  iBC    (nshg)                : BC codes
c  BC     (nshg,ndofBC)         : BC constraint parameters
c  iper   (nshg)                : periodicity table
c
c
c Zdenek Johan,  Winter 1991.  (Fortran 90)
c Alberto Figueroa, Winter 2004.  CMM-FSI
c Irene Vignon, Fall 2004. Impedance BC
c----------------------------------------------------------------------
c
      use pvsQbi     !gives us splag (the spmass at the end of this run
      use specialBC !gives us itvn
      use timedata   !allows collection of time series
      use convolImpFlow !for Imp bc
      use convolRCRFlow !for RCR bc
!MR CHANGE
      use turbsa          ! used to access d2wall
!MR CHANGE END
      use fncorpmod
      use iso_c_binding

c      use readarrays !reads in uold and acold

        include "common.h"
        include "mpif.h"
        include "auxmpi.h"
#ifdef HAVE_SVLS
        include "svLS.h"
#endif
#if !defined(HAVE_SVLS) && !defined(HAVE_ACUSOLVE)
#error "You must enable a linear solver"
#endif

c


        real*8    y(nshg,ndof),              ac(nshg,ndof),
     &            yold(nshg,ndof),           acold(nshg,ndof),
     &            u(nshg,nsd),               uold(nshg,nsd),
     &            x(numnp,nsd),              solinc(nshg,ndof),
     &            BC(nshg,ndofBC),           tf(nshg,ndof),
     &            GradV(nshg,nsdsq)

c
        real*8    res(nshg,ndof)
c
        real*8    shp(MAXTOP,maxsh,MAXQPT),
     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &            shpb(MAXTOP,maxsh,MAXQPT),
     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
        integer   rowp(nshg,nnz),         colm(nshg+1),
     &            iBC(nshg),
     &            ilwork(nlwork),
     &            iper(nshg),            ifuncs(6)

        real*8 vbc_prof(nshg,3)

        integer stopjob
        character*10 cname2
        character*5  cname
c
c  stuff for dynamic model s.w.avg and wall model
c
        dimension ifath(numnp),    velbar(nfath,ndof),  nsons(nfath)

        dimension wallubar(2),walltot(2)
c
c.... For linear solver Library
c
        integer eqnType, prjFlag, presPrjFlag, verbose
c
        real*8, allocatable, dimension(:,:) :: aperm,  atemp, atempS
        real*8, allocatable, dimension(:,:,:) :: apermS

        real*8, allocatable, dimension(:,:) :: lhsP, lhsK, lhsS
        real*8   almit, alfit, gamit
c
        character*1024    servername
        character*20    fname1,fmt1
        character*20    fname2,fmt2
        character*60    fnamepold, fvarts
        character*4     fname4c ! 4 characters
        integer         iarray(50) ! integers for headers
        integer         isgn(ndof), isgng(ndof)

!MR CHANGE
!        real*8 rerr(nshg,10), ybar(nshg,13) ! including 7 sq. terms with 3 cross terms of uv, uw and vw
!        real*8 rerr(nshg,10), ybar(nshg,12) ! including 7 sq. terms with 3 cross terms of uv, uw and vw
        real*8 rerr(nshg,10)
        real*8, allocatable, dimension(:,:) :: ybar, strain, vorticity
        real*8, allocatable, dimension(:,:) :: wallssVec, wallssVecbar
!MR CHANGE

        real*8 tcorecp(2), tcorecpscal(2)

        integer, allocatable, dimension(:) :: ivarts
        integer, allocatable, dimension(:) :: ivartsg
        integer, allocatable, dimension(:) :: iv_rank
        real*8, allocatable, dimension(:) :: vartssoln
        real*8, allocatable, dimension(:) :: vartssolng
        real*8, allocatable, dimension(:,:,:) :: yphbar
        real*8 CFLworst(numel)

!MR CHANGE
        integer :: iv_rankpernode, iv_totnodes, iv_totcores
        integer :: iv_node, iv_core, iv_thread
!MR CHANGE
!--------------------------------------------------------------------
!     Setting up svLS
#ifdef HAVE_SVLS
      INTEGER svLS_nFaces, gnNo, nNo, faIn, facenNo
      INTEGER, ALLOCATABLE :: gNodes(:)
      REAL*8, ALLOCATABLE :: sV(:,:)

      TYPE(svLS_commuType) communicator
      TYPE(svLS_lhsType) svLS_lhs
      TYPE(svLS_lsType) svLS_ls
! repeat for scalar solves (up to 4 at this time which is consistent with rest of PHASTA)
      TYPE(svLS_commuType) communicator_S(4)
      TYPE(svLS_lhsType) svLS_lhs_S(4)
      TYPE(svLS_lsType) svLS_ls_S(4)
#endif

        impistat = 0
        impistat2 = 0
        iISend = 0
        iISendScal = 0
        iIRecv = 0
        iIRecvScal = 0
        iWaitAll = 0
        iWaitAllScal = 0
        iAllR = 0
        iAllRScal = 0
        rISend = zero
        rISendScal = zero
        rIRecv = zero
        rIRecvScal = zero
        rWaitAll = zero
        rWaitAllScal = zero
        rAllR = zero
        rAllRScal = zero
        rCommu = zero
        rCommuScal = zero

        call initmemstat()

c
c  hack SA variable
c
cHack        BC(:,7)=BC(:,7)*0.001
cHack        if(lstep.eq.0) y(:,6)=y(:,6)*0.001
!--------------------------------------------------------------------
!     Setting up svLS Moved down for better org

#ifdef HAVE_ACUSOLVE
      IF (svLSFlag .EQ. 0) THEN  !When we get a PETSc option it also could block this or a positive leslib
        call SolverLicenseServer(servername)
      ENDIF
#endif
c
c only master should be verbose
c

        if(numpe.gt.0 .and. myrank.ne.master)iverbose=0
c

        lskeep=lstep

        inquire(file='xyzts.dat',exist=exts)
        if(exts) then

           open(unit=626,file='xyzts.dat',status='old')
           read(626,*) ntspts, freq, tolpt, iterat, varcod
           call sTD             ! sets data structures

           do jj=1,ntspts       ! read coordinate data where solution desired
              read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3)
           enddo
           close(626)

           statptts(:,:) = 0
           parptts(:,:) = zero
           varts(:,:) = zero

           allocate (ivarts(ntspts*ndof))
           allocate (ivartsg(ntspts*ndof))
           allocate (iv_rank(ntspts))
           allocate (vartssoln(ntspts*ndof))
           allocate (vartssolng(ntspts*ndof))

           iv_rankpernode = iv_rankpercore*iv_corepernode
           iv_totnodes = numpe/iv_rankpernode
           iv_totcores = iv_corepernode*iv_totnodes
           if (myrank .eq. 0) then
             write(*,*) 'Info for probes:'
             write(*,*) '  Ranks per core:',iv_rankpercore
             write(*,*) '  Cores per node:',iv_corepernode
             write(*,*) '  Ranks per node:',iv_rankpernode
             write(*,*) '  Total number of nodes:',iv_totnodes
             write(*,*) '  Total number of cores',iv_totcores
           endif

!           if (myrank .eq. numpe-1) then
            do jj=1,ntspts

               ! Compute the adequate rank which will take care of probe jj
               jjm1 = jj-1
               iv_node = (iv_totnodes-1)-mod(jjm1,iv_totnodes)
               iv_core = (iv_corepernode-1) - mod((jjm1 -
     &              mod(jjm1,iv_totnodes))/iv_totnodes,iv_corepernode)
               iv_thread = (iv_rankpercore-1) - mod((jjm1-
     &              (mod(jjm1,iv_totcores)))/iv_totcores,iv_rankpercore)
               iv_rank(jj) = iv_node*iv_rankpernode
     &                     + iv_core*iv_rankpercore
     &                     + iv_thread

               if(myrank == 0) then
                 write(*,*) '  Probe', jj, 'handled by rank',
     &                         iv_rank(jj), ' on node', iv_node
               endif

               ! Verification just in case
               if(iv_rank(jj) .lt.0 .or. iv_rank(jj) .ge. numpe) then
                 write(*,*) 'WARNING: iv_rank(',jj,') is ', iv_rank(jj),
     &                      ' and reset to numpe-1'
                 iv_rank(jj) = numpe-1
               endif

               ! Open the varts files
               if(myrank == iv_rank(jj)) then
                 fvarts='varts/varts'
                 fvarts=trim(fvarts)//trim(cname2(jj))
                 fvarts=trim(fvarts)//trim(cname2(lstep))
                 fvarts=trim(fvarts)//'.dat'
                 fvarts=trim(fvarts)
                 open(unit=1000+jj, file=fvarts, status='unknown')
               endif
            enddo
!           endif

        endif
c
c.... open history and aerodynamic forces files
c
        if (myrank .eq. master) then
           open (unit=ihist,  file=fhist,  status='unknown')
           open (unit=iforce, file=fforce, status='unknown')
           open (unit=76, file="fort.76", status='unknown')
           if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then
              fnamepold = 'pold'
              fnamepold = trim(fnamepold)//trim(cname2(lstep))
              fnamepold = trim(fnamepold)//'.dat'
              fnamepold = trim(fnamepold)
              open (unit=8177, file=fnamepold, status='unknown')
           endif
        endif
c
c.... initialize
c
        ifuncs(:)  = 0              ! func. evaluation counter
        istep  = 0
        yold   = y
        acold  = ac

        rerr = zero

        if(ierrcalc.eq.1 .or. ioybar.eq.1) then ! we need ybar for error too
          if (ivort == 1) then
            allocate(ybar(nshg,17)) ! more space for vorticity if requested
          else
            allocate(ybar(nshg,13))
          endif
          ybar = zero ! Initialize ybar to zero, which is essential
        endif

        if(ivort == 1) then
          allocate(strain(nshg,6))
          allocate(vorticity(nshg,5))
        endif

        if(abs(itwmod).ne.1 .and. iowflux.eq.1) then
          allocate(wallssVec(nshg,3))
          if (ioybar .eq. 1) then
            allocate(wallssVecbar(nshg,3))
            wallssVecbar = zero ! Initialization important if mean wss computed
          endif
        endif

! both nstepsincycle and nphasesincycle needs to be set
        if(nstepsincycle.eq.0) nphasesincycle = 0
        if(nphasesincycle.ne.0) then
!     &     allocate(yphbar(nshg,5,nphasesincycle))
          if (ivort == 1) then
            allocate(yphbar(nshg,15,nphasesincycle)) ! more space for vorticity
          else
            allocate(yphbar(nshg,11,nphasesincycle))
          endif
          yphbar = zero
        endif

!MR CHANGE END

        vbc_prof(:,1:3) = BC(:,3:5)
        if(iramp.eq.1) then
          call BCprofileInit(vbc_prof,x)
        endif

c
c.... ---------------> initialize LesLib Library <---------------
c
c.... assign parameter values
c
        do i = 1, 100
           numeqns(i) = i
        enddo
        nKvecs       = Kspace
        prjFlag      = iprjFlag
        presPrjFlag  = ipresPrjFlag
        verbose      = iverbose
c
c.... determine how many scalar equations we are going to need to solve
c
      nsolt=mod(impl(1),2)      ! 1 if solving temperature
      nsclrsol=nsolt+nsclr      ! total number of scalars solved At
                                ! some point we probably want to create
                                ! a map, considering stepseq(), to find
                                ! what is actually solved and only
                                ! dimension lhs to the appropriate
                                ! size. (see 1.6.1 and earlier for a
                                ! "failed" attempt at this).


      nsolflow=mod(impl(1),100)/10  ! 1 if solving flow

c
c.... Now, call lesNew routine to initialize
c     memory space
c
      call genadj(colm, rowp, icnt )  ! preprocess the adjacency list

      nnz_tot=icnt ! this is exactly the number of non-zero blocks on
                   ! this proc

      if (nsolflow.eq.1) then
         lesId   = numeqns(1)
         eqnType = 1
         nDofs   = 4

!--------------------------------------------------------------------
!     Rest of configuration of svLS is added here, where we have LHS
!     pointers

         nPermDims = 1
         nTmpDims = 1

         allocate (lhsP(4,nnz_tot))
         allocate (lhsK(9,nnz_tot))

!     Setting up svLS or leslib for flow
#ifdef HAVE_SVLS
      IF (svLSFlag .EQ. 1) THEN
         IF(nPrjs.eq.0) THEN
           svLSType=2  !GMRES if borrowed ACUSIM projection vectors variable set to zero
         ELSE
           svLSType=3 !NS solver
         ENDIF
!  reltol for the NSSOLVE is the stop criterion on the outer loop
!  reltolIn is (eps_GM, eps_CG) from the CompMech paper
!  for now we are using
!  Tolerance on ACUSIM Pressure Projection for CG and
!  Tolerance on Momentum Equations for GMRES
! also using Kspaceand maxIters from setup for ACUSIM
!
         eps_outer=40.0*epstol(1)  !following papers soggestion for now
         CALL svLS_LS_CREATE(svLS_ls, svLSType, dimKry=Kspace,
     2      relTol=eps_outer, relTolIn=(/epstol(1),prestol/),
     3      maxItr=maxIters,
     4      maxItrIn=(/maxIters,maxIters/))

         CALL svLS_COMMU_CREATE(communicator, MPI_COMM_WORLD)
            nNo=nshg
            gnNo=nshgt
            IF  (ipvsq .GE. 2) THEN

#if((VER_CORONARY == 1)&&(VER_CLOSEDLOOP == 1))
               svLS_nFaces = 1 + numResistSrfs + numNeumannSrfs
     2            + numImpSrfs + numRCRSrfs + numCORSrfs
#elif((VER_CORONARY == 1)&&(VER_CLOSEDLOOP == 0))
               svLS_nFaces = 1 + numResistSrfs
     2            + numImpSrfs + numRCRSrfs + numCORSrfs
#elif((VER_CORONARY == 0)&&(VER_CLOSEDLOOP == 1))
               svLS_nFaces = 1 + numResistSrfs + numNeumannSrfs
     2            + numImpSrfs + numRCRSrfs
#else
               svLS_nFaces = 1 + numResistSrfs
     2            + numImpSrfs + numRCRSrfs
#endif

            ELSE
               svLS_nFaces = 1   !not sure about this...looks like 1 means 0 for array size issues
            END IF

            CALL svLS_LHS_CREATE(svLS_lhs, communicator, gnNo, nNo,
     2         nnz_tot, ltg, colm, rowp, svLS_nFaces)

            faIn = 1
            facenNo = 0
            DO i=1, nshg
               IF (IBITS(iBC(i),3,3) .NE. 0)  facenNo = facenNo + 1
            END DO
            ALLOCATE(gNodes(facenNo), sV(nsd,facenNo))
            sV = 0D0
            j = 0
            DO i=1, nshg
               IF (IBITS(iBC(i),3,3) .NE. 0) THEN
                  j = j + 1
                  gNodes(j) = i
                  IF (.NOT.BTEST(iBC(i),3)) sV(1,j) = 1D0
                  IF (.NOT.BTEST(iBC(i),4)) sV(2,j) = 1D0
                  IF (.NOT.BTEST(iBC(i),5)) sV(3,j) = 1D0
               END IF
            END DO
            CALL svLS_BC_CREATE(svLS_lhs, faIn, facenNo,
     2         nsd, BC_TYPE_Dir, gNodes, sV)
            DEALLOCATE(gNodes)
            DEALLOCATE(sV)
#endif
#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE)
         ELSE ! leslib solve
#elif defined(HAVE_SVLS)
         ENDIF
#endif
#ifdef HAVE_ACUSOLVE
!--------------------------------------------------------------------
           call myfLesNew( lesId,   41994,
     &                 eqnType,
     &                 nDofs,          minIters,       maxIters,
     &                 nKvecs,         prjFlag,        nPrjs,
     &                 presPrjFlag,    nPresPrjs,      epstol(1),
     &                 prestol,        verbose,        statsflow,
     &                 nPermDims,      nTmpDims,      servername  )

           allocate (aperm(nshg,nPermDims))
           allocate (atemp(nshg,nTmpDims))
           call readLesRestart( lesId,  aperm, nshg, myrank, lstep,
     &                        nPermDims )
#endif
#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE)
         ENDIF !flow solver selector
#endif
      else   ! not solving flow just scalar
         nPermDims = 0
         nTempDims = 0
      endif


      if(nsclrsol.gt.0) then
       do isolsc=1,nsclrsol
         lesId       = numeqns(isolsc+1)
         eqnType     = 2
         nDofs       = 1
         presPrjflag = 0
         nPresPrjs   = 0
         prjFlag     = 1
         indx=isolsc+2-nsolt ! complicated to keep epstol(2) for
                             ! temperature followed by scalars
!     Setting up svLS or leslib for scalar
#ifdef HAVE_SVLS
      IF (svLSFlag .EQ. 1) THEN
           svLSType=2  !only option for scalars
!  reltol for the GMRES is the stop criterion
! also using Kspaceand maxIters from setup for ACUSIM
!
         CALL svLS_LS_CREATE(svLS_ls_S(isolsc), svLSType, dimKry=Kspace,
     2      relTol=epstol(indx),
     3      maxItr=maxIters
     4      )

         CALL svLS_COMMU_CREATE(communicator_S(isolsc), MPI_COMM_WORLD)

               svLS_nFaces = 1   !not sure about this...should try it with zero

            CALL svLS_LHS_CREATE(svLS_lhs_S(isolsc), communicator_S(isolsc), gnNo, nNo,
     2         nnz_tot, ltg, colm, rowp, svLS_nFaces)

              faIn = 1
              facenNo = 0
              ib=5+isolsc
              DO i=1, nshg
                 IF (btest(iBC(i),ib))  facenNo = facenNo + 1
              END DO
              ALLOCATE(gNodes(facenNo), sV(1,facenNo))
              sV = 0D0
              j = 0
              DO i=1, nshg
               IF (btest(iBC(i),ib)) THEN
                  j = j + 1
                  gNodes(j) = i
               END IF
              END DO

            CALL svLS_BC_CREATE(svLS_lhs_S(isolsc), faIn, facenNo,
     2         1, BC_TYPE_Dir, gNodes, sV(1,:))
            DEALLOCATE(gNodes)
            DEALLOCATE(sV)

            if( isolsc.eq.1) then ! if multiple scalars make sure done once
              allocate (apermS(1,1,1))
              allocate (atempS(1,1))  !they can all share this
            endif
#endif
#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE)
         ELSE ! leslib solve of scalar
#elif defined(HAVE_SVLS)
         ENDIF
#endif
#ifdef HAVE_ACUSOLVE
         call myfLesNew( lesId,            41994,
     &                 eqnType,
     &                 nDofs,          minIters,       maxIters,
     &                 nKvecs,         prjFlag,        nPrjs,
     &                 presPrjFlag,    nPresPrjs,      epstol(indx),
     &                 prestol,        verbose,        statssclr,
     &                 nPermDimsS,     nTmpDimsS,   servername )
#endif
#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE)
        ENDIF
#endif
       enddo  !loop over scalars to solve  (not yet worked out for multiple svLS solves
       allocate (lhsS(nnz_tot,nsclrsol))
#ifdef HAVE_ACUSOLVE
       if(svLSFlag.eq.0) then  ! we just prepped scalar solves for leslib so allocate arrays
c
c  Assume all scalars have the same size needs
c
       allocate (apermS(nshg,nPermDimsS,nsclrsol))
       allocate (atempS(nshg,nTmpDimsS))  !they can all share this
       endif
#endif
c
c actually they could even share with atemp but leave that for later
c
      else !no scalar solves at all so zero dims not used
         nPermDimsS = 0
         nTmpDimsS  = 0
      endif
c
c...  prepare lumped mass if needed
c
      if((flmpr.ne.0).or.(flmpl.ne.0)) call genlmass(x, shp,shgl)
c
c.... -----------------> End of initialization <-----------------
c
c.....open the necessary files to gather time series
c
      lstep0 = lstep+1
      nsteprcr = nstep(1)+lstep
c
c.... loop through the time sequences
c


      do 3000 itsq = 1, ntseq
         itseq = itsq

CAD         tcorecp1 = second(0)
CAD         tcorewc1 = second(-1)
c
c.... set up the time integration parameters
c
         nstp   = nstep(itseq)
         nitr   = niter(itseq)
         LCtime = loctim(itseq)
         dtol(:)= deltol(itseq,:)

         call itrSetup ( y, acold )
c
c...initialize the coefficients for the impedance convolution,
c   which are functions of alphaf so need to do it after itrSetup
         if(numImpSrfs.gt.zero) then
            call calcImpConvCoef (numImpSrfs, ntimeptpT)
         endif
c
c...initialize the initial condition P(0)-RQ(0)-Pd(0) for RCR BC
c   need ndsurf so should be after initNABI
         if(numRCRSrfs.gt.zero) then
            call calcRCRic(y,nsrflistRCR,numRCRSrfs)
         endif
c
c  find the last solve of the flow in the step sequence so that we will
c         know when we are at/near end of step
c
c         ilast=0
         nitr=0  ! count number of flow solves in a step (# of iterations)
         do i=1,seqsize
            if(stepseq(i).eq.0) nitr=nitr+1
         enddo

         if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         tcorecp(:) = zero ! used in solfar.f (solflow)
         tcorecpscal(:) = zero ! used in solfar.f (solflow)
         if(myrank.eq.0)  then
            tcorecp1 = TMRC()
         endif

c
c.... loop through the time steps
c
         istop=0
         rmub=datmat(1,2,1)
         if(rmutarget.gt.0) then
            rmue=rmutarget
         else
            rmue=datmat(1,2,1) ! keep constant
         endif

        if(iramp.eq.1) then
            call BCprofileScale(vbc_prof,BC,yold) ! fix the yold values to the reset BC
            isclr=1 ! fix scalar
            do isclr=1,nsclr
               call itrBCSclr(yold,ac,iBC,BC,iper,ilwork)
            enddo
        endif

         do 2000 istp = 1, nstp
           if(iramp.eq.1)
     &        call BCprofileScale(vbc_prof,BC,yold)

           call rerun_check(stopjob)
           if(myrank.eq.master) write(*,*)
     &         'stopjob,lstep,istep', stopjob,lstep,istep
           if(stopjob.eq.lstep) then
              stopjob=-2 ! this is the code to finish
             if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then
                if(myrank.eq.master) write(*,*)
     &         'line 473 says last step written so exit'
                goto 2002  ! the step was written last step so just exit
             else
                if(myrank.eq.master)
     &         write(*,*) 'line 473 says last step not written'
                istep=nstp  ! have to do this so that solution will be written
                goto 2001
             endif
           endif

            xi=istp*1.0/nstp
            datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue
c            write(*,*) "current mol. visc = ", datmat(1,2,1)
c.... if we have time varying boundary conditions update the values of BC.
c     these will be for time step n+1 so use lstep+1
c
            if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl,
     &                               shpb, shglb, x, BC, iBC)

c
c ... calculate the pressure contribution that depends on the history for the Imp. BC
c
            if(numImpSrfs.gt.0) then
               call pHist(poldImp,QHistImp,ImpConvCoef,
     &                    ntimeptpT,numImpSrfs)
               if (myrank.eq.master)
     &             write(8177,*) (poldImp(n), n=1,numImpSrfs)
            endif
c
c ... calc the pressure contribution that depends on the history for the RCR BC
c
            if(numRCRSrfs.gt.0) then
               call CalcHopRCR (Delt(itseq), lstep, numRCRSrfs)
               call CalcRCRConvCoef(lstep,numRCRSrfs)
               call pHist(poldRCR,QHistRCR,RCRConvCoef,nsteprcr,
     &              numRCRSrfs)
               if (myrank.eq.master)
     &             write(8177,*) (poldRCR(n), n=1,numRCRSrfs)
            endif
c
c Decay of scalars
c
           if(nsclr.gt.0 .and. tdecay.ne.1) then
              yold(:,6:ndof)=y(:,6:ndof)*tdecay
              BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay
           endif

           if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8


            if(iLES.gt.0) then  !complicated stuff has moved to
                                        !routine below
               call lesmodels(yold,  acold,     shgl,      shp,
     &                        iper,  ilwork,    rowp,      colm,
     &                        nsons, ifath,     x,
     &                        iBC,   BC)


            endif

c.... set traction BCs for modeled walls
c
            if (itwmod.ne.0) then
               call asbwmod(yold,   acold,   x,      BC,     iBC,
     &                      iper,   ilwork,  ifath,  velbar)
            endif

!MR CHANGE
c
c.... Determine whether the vorticity field needs to be computed for this time step or not
c
            icomputevort = 0
            if (ivort == 1) then ! Print vorticity = True in solver.inp
              ! We then compute the vorticity only if we
              ! 1) we write an intermediate checkpoint
              ! 2) we reach the last time step and write the last checkpoint
              ! 3) we accumulate statistics in ybar for every time step
              ! BEWARE: we need here lstep+1 and istep+1 because the lstep and
              ! istep gets incremened after the flowsolve, further below
              if (((irs .ge. 1) .and. (mod(lstep+1, ntout) .eq. 0)) .or.
     &                   istep+1.eq.nstep(itseq) .or. ioybar == 1) then
                icomputevort = 1
              endif
            endif

!            write(*,*) 'icomputevort: ',icomputevort, ' - istep: ',
!     &                istep,' - nstep(itseq):',nstep(itseq),'- lstep:',
!     &                lstep, '- ntout:', ntout
!MR CHANGE

c
c.... -----------------------> predictor phase <-----------------------
c
            call itrPredict(yold, y,   acold,  ac ,  uold,  u, iBC)
            call itrBC (y,  ac,  iBC,  BC,  iper,ilwork)

            if(nsolt.eq.1) then
               isclr=0
               call itrBCSclr (y, ac,  iBC, BC, iper, ilwork)
            endif
            do isclr=1,nsclr
               call itrBCSclr (y, ac,  iBC, BC, iper, ilwork)
            enddo
            iter=0
            ilss=0  ! this is a switch thrown on first solve of LS redistance
            do istepc=1,seqsize
               icode=stepseq(istepc)
               if(mod(icode,5).eq.0) then ! this is a solve
                  isolve=icode/10
                  if(icode.eq.0) then ! flow solve (encoded as 0)
c
                     iter   = iter+1
                     ifuncs(1)  = ifuncs(1) + 1
c
                     Force(1) = zero
                     Force(2) = zero
                     Force(3) = zero
                     HFlux    = zero
                     lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1)))

                     call SolFlow(y,          ac,        u,
     &                         yold,          acold,     uold,
     &                         x,             iBC,
     &                         BC,            res,
     &                         nPermDims,     nTmpDims,  aperm,
     &                         atemp,         iper,
     &                         ilwork,        shp,       shgl,
     &                         shpb,          shglb,     rowp,
     &                         colm,          lhsK,      lhsP,
     &                         solinc,        rerr,      tcorecp,
     &                         GradV,      sumtime
#ifdef HAVE_SVLS
     &                         ,svLS_lhs,     svLS_ls,  svLS_nFaces)
#else
     &                         )
#endif

                  else          ! scalar type solve
                     if (icode.eq.5) then ! Solve for Temperature
                                ! (encoded as (nsclr+1)*10)
                        isclr=0
                        ifuncs(2)  = ifuncs(2) + 1
                        j=1
                     else       ! solve a scalar  (encoded at isclr*10)
                        isclr=isolve
                        ifuncs(isclr+2)  = ifuncs(isclr+2) + 1
                        j=isclr+nsolt
                        if((iLSet.eq.2).and.(ilss.eq.0)
     &                       .and.(isclr.eq.2)) then
                           ilss=1 ! throw switch (once per step)
                           y(:,7)=y(:,6) ! redistance field initialized
                           ac(:,7)   = zero
                           call itrBCSclr (  y,  ac,  iBC,  BC, iper,
     &                          ilwork)
c
c....store the flow alpha, gamma parameter values and assigm them the
c....Backward Euler parameters to solve the second levelset scalar
c
                           alfit=alfi
                           gamit=gami
                           almit=almi
                           Deltt=Delt(1)
                           Dtglt=Dtgl
                           alfi = 1
                           gami = 1
                           almi = 1
c     Delt(1)= Deltt ! Give a pseudo time step
                           Dtgl = one / Delt(1)
                        endif  ! level set eq. 2
                     endif ! deciding between temperature and scalar

                     lhs = 1 - min(1,mod(ifuncs(isclr+2)-1,
     &                                   LHSupd(isclr+2)))

                     call SolSclr(y,          ac,        u,
     &                         yold,          acold,     uold,
     &                         x,             iBC,
     &                         BC,            nPermDimsS,nTmpDimsS,
     &                         apermS(1,1,j), atempS,    iper,
     &                         ilwork,        shp,       shgl,
     &                         shpb,          shglb,     rowp,
     &                         colm,          lhsS(1,j),
     &                         solinc(1,isclr+5), tcorecpscal
#ifdef HAVE_SVLS
     &                         ,svLS_lhs_S(isclr),   svLS_ls_S(isclr), svls_nfaces)
#else
     &                         )
#endif


                  endif         ! end of scalar type solve

               else ! this is an update  (mod did not equal zero)
                  iupdate=icode/10  ! what to update
                  if(icode.eq.1) then !update flow
                     call itrCorrect ( y,    ac,    u,   solinc, iBC)
                     call itrBC (y,  ac,  iBC,  BC, iper, ilwork)
                  else  ! update scalar
                     isclr=iupdate  !unless
                     if(icode.eq.6) isclr=0
                     if(iRANS.lt.-100)then  ! RANS
                        call itrCorrectSclrPos(y,ac,solinc(1,isclr+5))
                     else
                        call itrCorrectSclr (y, ac, solinc(1,isclr+5))
                     endif
                     if (ilset.eq.2 .and. isclr.eq.2)  then
                        if (ivconstraint .eq. 1) then
                           call itrBCSclr (  y,  ac,  iBC,  BC, iper,
     &                          ilwork)
c
c ... applying the volume constraint on second level set scalar
c
                           call solvecon (y,    x,      iBC,  BC,
     &                          iper, ilwork, shp,  shgl)
c
                        endif   ! end of volume constraint calculations
                     endif      ! end of redistance calculations
c
                        call itrBCSclr (  y,  ac,  iBC,  BC, iper,
     &                       ilwork)
                     endif      ! end of flow or scalar update
                  endif         ! end of switch between solve or update
               enddo            ! loop over sequence in step
c
c
c.... obtain the time average statistics
c
            if (ioform .eq. 2) then

               call stsGetStats( y,      yold,     ac,     acold,
     &                           u,      uold,     x,
     &                           shp,    shgl,     shpb,   shglb,
     &                           iBC,    BC,       iper,   ilwork,
     &                           rowp,   colm,     lhsK,   lhsP )
            endif

c
c  Find the solution at the end of the timestep and move it to old
c
c
c ...First to reassign the parameters for the original time integrator scheme
c
            if((iLSet.eq.2).and.(ilss.eq.1)) then
               alfi =alfit
               gami =gamit
               almi =almit
               Delt(1)=Deltt
               Dtgl =Dtglt
            endif
            call itrUpdate( yold,  acold,   uold,  y,    ac,   u)
            call itrBC (yold, acold,  iBC,  BC,  iper,ilwork)

            istep = istep + 1
            lstep = lstep + 1
c
c ..  Print memory consumption on BGQ
c
            call printmeminfo("itrdrv"//char(0))

c
c ..  Compute vorticity
c
            if ( icomputevort == 1) then

              ! vorticity components and magnitude
              vorticity(:,1) = GradV(:,8)-GradV(:,6) !omega_x
              vorticity(:,2) = GradV(:,3)-GradV(:,7) !omega_y
              vorticity(:,3) = GradV(:,4)-GradV(:,2) !omega_z
              vorticity(:,4) = sqrt(   vorticity(:,1)*vorticity(:,1)
     &                               + vorticity(:,2)*vorticity(:,2)
     &                               + vorticity(:,3)*vorticity(:,3) )
              ! Q
              strain(:,1) = GradV(:,1)                  !S11
              strain(:,2) = 0.5*(GradV(:,2)+GradV(:,4)) !S12
              strain(:,3) = 0.5*(GradV(:,3)+GradV(:,7)) !S13
              strain(:,4) = GradV(:,5)                  !S22
              strain(:,5) = 0.5*(GradV(:,6)+GradV(:,8)) !S23
              strain(:,6) = GradV(:,9)                  !S33

              vorticity(:,5) = 0.25*( vorticity(:,4)*vorticity(:,4)  !Q
     &                            - 2.0*(      strain(:,1)*strain(:,1)
     &                                    + 2* strain(:,2)*strain(:,2)
     &                                    + 2* strain(:,3)*strain(:,3)
     &                                    +    strain(:,4)*strain(:,4)
     &                                    + 2* strain(:,5)*strain(:,5)
     &                                    +    strain(:,6)*strain(:,6)))

            endif
c
c .. write out the instantaneous solution
c
2001    continue  ! we could get here by 2001 label if user requested stop
        if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or.
     &      istep.eq.nstep(itseq)) then

!so that we can see progress in force file close it so that it flushes
!and  then reopen in append mode

           close(iforce)
           open (unit=iforce, file=fforce, position='append')

!              Call to restar() will open restart file in write mode (and not append mode)
!              that is needed as other fields are written in append mode

           call restar ('out ',  yold  ,ac)
           if(ideformwall == 1) then
             if(myrank.eq.master) then
               write(*,*) 'ideformwall is 1 and is a dead code path... exiting'
             endif
             stop
           endif

           if(ivort == 1) then
             call write_field(myrank,'a','vorticity',9,vorticity,
     &                       'd',nshg,5,lstep)
           endif

           call printmeminfo("itrdrv after checkpoint"//char(0))
         else if(stopjob.eq.-2) then
           if(myrank.eq.master) then
             write(*,*) 'line 755 says no write before stopping'
             write(*,*) 'istep,nstep,irs',istep,nstep(itseq),irs
           endif
        endif  !just the instantaneous stuff for videos
c
c.... compute the consistent boundary flux
c
            if(abs(itwmod).ne.1 .and. iowflux.eq.1) then
               call Bflux ( yold,      acold,      uold,    x,
     &                      shp,       shgl,       shpb,
     &                      shglb,     ilwork,     iBC,
     &                      BC,        iper,       wallssVec)
            endif

           if(stopjob.eq.-2) goto 2003


c
c ... update the flow history for the impedance convolution, filter it and write it out
c
            if(numImpSrfs.gt.zero) then
               call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1)
            endif

c
c ... update the flow history for the RCR convolution
c
            if(numRCRSrfs.gt.zero) then
               call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep
            endif


c...  dump TIME SERIES

            if (exts) then
               if (mod(lstep-1,freq).eq.0) then

                  if (numpe > 1) then
                     do jj = 1, ntspts
                        vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:)
                        ivarts=zero
                     enddo
                     do k=1,ndof*ntspts
                        if(vartssoln(k).ne.zero) ivarts(k)=1
                     enddo

!                     call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts,
!     &                    MPI_DOUBLE_PRECISION, MPI_SUM, master,
!     &                    MPI_COMM_WORLD, ierr)

                     call MPI_BARRIER(MPI_COMM_WORLD, ierr)
                     call MPI_ALLREDUCE(vartssoln, vartssolng,
     &                    ndof*ntspts,
     &                    MPI_DOUBLE_PRECISION, MPI_SUM,
     &                    MPI_COMM_WORLD, ierr)

!                     call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts,
!     &                    MPI_INTEGER, MPI_SUM, master,
!     &                    MPI_COMM_WORLD, ierr)

                     call MPI_BARRIER(MPI_COMM_WORLD, ierr)
                     call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts,
     &                    MPI_INTEGER, MPI_SUM,
     &                    MPI_COMM_WORLD, ierr)

!                     if (myrank.eq.zero) then
                     do jj = 1, ntspts

                        if(myrank .eq. iv_rank(jj)) then
                           ! No need to update all varts components, only the one treated by the expected rank
                           ! Note: keep varts as a vector, as multiple probes could be treated by one rank
                           indxvarts = (jj-1)*ndof
                           do k=1,ndof
                              if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero
                                 varts(jj,k)=vartssolng(indxvarts+k)/
     &                                             ivartsg(indxvarts+k)
                              endif
                           enddo
                       endif !only if myrank eq iv_rank(jj)
                     enddo
!                     endif !only on master
                  endif !only if numpe > 1

!                  if (myrank.eq.zero) then
                  do jj = 1, ntspts
                     if(myrank .eq. iv_rank(jj)) then
                        ifile = 1000+jj
                        write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof
c                        call flush(ifile)
                        if (((irs .ge. 1) .and.
     &                       (mod(lstep, ntout) .eq. 0))) then
                           close(ifile)
                           fvarts='varts/varts'
                           fvarts=trim(fvarts)//trim(cname2(jj))
                           fvarts=trim(fvarts)//trim(cname2(lskeep))
                           fvarts=trim(fvarts)//'.dat'
                           fvarts=trim(fvarts)
                           open(unit=ifile, file=fvarts,
     &                          position='append')
                        endif !only when dumping restart
                     endif
                  enddo
!                  endif !only on master

                  varts(:,:) = zero ! reset the array for next step


!555              format(i6,5(2x,E12.5e2))
555               format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here

               endif
            endif

c
c.... update and the aerodynamic forces
c
            call forces ( yold,  ilwork )

            if((irscale.ge.0).or.(itwmod.gt.0))
     &           call getvel (yold,     ilwork, iBC,
     &                        nsons,    ifath, velbar)

            if((irscale.ge.0).and.(myrank.eq.master)) then
               call genscale(yold,       x,       iper,
     &                       iBC,     ifath,   velbar,
     &                       nsons)
            endif
c
c....  print out results.
c
            ntoutv=max(ntout,100)   ! velb is not needed so often
            if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then
               if( (mod(lstep, ntoutv) .eq. 0) .and.
     &              ((irscale.ge.0).or.(itwmod.gt.0) .or.
     &              ((nsonmax.eq.1).and.(iLES.gt.0))))
     &              call rwvelb  ('out ',  velbar  ,ifail)
            endif
c
c.... end of the NSTEP and NTSEQ loops
c


c
c.... -------------------> error calculation  <-----------------
c
            if(ierrcalc.eq.1 .or. ioybar.eq.1) then
c$$$c
c$$$c compute average
c$$$c
c$$$               tfact=one/istep
c$$$               ybar =tfact*yold + (one-tfact)*ybar

c compute average
c ybar(:,1:3) are average velocity components
c ybar(:,4) is average pressure
c ybar(:,5) is average speed
c ybar(:,6:8) is average of sq. of each vel. component
c ybar(:,9) is average of sq. of pressure
c ybar(:,10:12) is average of cross vel. components : uv, uw and vw
c averaging procedure justified only for identical time step sizes
c ybar(:,13) is average of eddy viscosity
c ybar(:,14:16) is average vorticity components
c ybar(:,17) is average vorticity magnitude
c istep is number of time step
c
               icollectybar = 0
          if(nphasesincycle.eq.0 .or.
     &            istep.gt.ncycles_startphaseavg*nstepsincycle) then
                 icollectybar = 1
                  if((istep-1).eq.ncycles_startphaseavg*nstepsincycle)
     &               istepsinybar = 0 ! init. to zero in first cycle in avg.
               endif

               if(icollectybar.eq.1) then
                  istepsinybar = istepsinybar+1
                  tfact=one/istepsinybar

                  if(myrank.eq.master .and. nphasesincycle.ne.0 .and.
     &               mod((istep-1),nstepsincycle).eq.0)
     &               write(*,*)'nsamples in phase average:',istepsinybar

c ybar to contain the averaged ((u,v,w),p)-fields
c and speed average, i.e., sqrt(u^2+v^2+w^2)
c and avg. of sq. terms including
c u^2, v^2, w^2, p^2 and cross terms of uv, uw and vw

                  ybar(:,1) = tfact*yold(:,1) + (one-tfact)*ybar(:,1)
                  ybar(:,2) = tfact*yold(:,2) + (one-tfact)*ybar(:,2)
                  ybar(:,3) = tfact*yold(:,3) + (one-tfact)*ybar(:,3)
                  ybar(:,4) = tfact*yold(:,4) + (one-tfact)*ybar(:,4)
                  ybar(:,5) = tfact*sqrt(yold(:,1)**2+yold(:,2)**2+
     &                        yold(:,3)**2) + (one-tfact)*ybar(:,5)
                  ybar(:,6) = tfact*yold(:,1)**2 +
     &                        (one-tfact)*ybar(:,6)
                  ybar(:,7) = tfact*yold(:,2)**2 +
     &                        (one-tfact)*ybar(:,7)
                  ybar(:,8) = tfact*yold(:,3)**2 +
     &                        (one-tfact)*ybar(:,8)
                  ybar(:,9) = tfact*yold(:,4)**2 +
     &                        (one-tfact)*ybar(:,9)
                  ybar(:,10) = tfact*yold(:,1)*yold(:,2) + !uv
     &                         (one-tfact)*ybar(:,10)
                  ybar(:,11) = tfact*yold(:,1)*yold(:,3) + !uw
     &                         (one-tfact)*ybar(:,11)
                  ybar(:,12) = tfact*yold(:,2)*yold(:,3) + !vw
     &                         (one-tfact)*ybar(:,12)
!MR CHANGE
                  if(nsclr.gt.0) !nut
     &             ybar(:,13) = tfact*yold(:,6) + (one-tfact)*ybar(:,13)

                  if(ivort == 1) then !vorticity
                    ybar(:,14) = tfact*vorticity(:,1) +
     &                           (one-tfact)*ybar(:,14)
                    ybar(:,15) = tfact*vorticity(:,2) +
     &                           (one-tfact)*ybar(:,15)
                    ybar(:,16) = tfact*vorticity(:,3) +
     &                           (one-tfact)*ybar(:,16)
                    ybar(:,17) = tfact*vorticity(:,4) +
     &                           (one-tfact)*ybar(:,17)
                  endif

                  if(abs(itwmod).ne.1 .and. iowflux.eq.1) then
                    wallssVecBar(:,1) = tfact*wallssVec(:,1)
     &                                  +(one-tfact)*wallssVecBar(:,1)
                    wallssVecBar(:,2) = tfact*wallssVec(:,2)
     &                                  +(one-tfact)*wallssVecBar(:,2)
                    wallssVecBar(:,3) = tfact*wallssVec(:,3)
     &                                  +(one-tfact)*wallssVecBar(:,3)
                  endif
!MR CHANGE END
               endif
c
c compute phase average
c
               if(nphasesincycle.ne.0 .and.
     &            istep.gt.ncycles_startphaseavg*nstepsincycle) then

c beginning of cycle is considered as ncycles_startphaseavg*nstepsincycle+1
                  if((istep-1).eq.ncycles_startphaseavg*nstepsincycle)
     &               icyclesinavg = 0 ! init. to zero in first cycle in avg.

                  ! find number of steps between phases
                  nstepsbtwphase = nstepsincycle/nphasesincycle ! integer value
                  if(mod(istep-1,nstepsincycle).eq.0) then
                     iphase = 1 ! init. to one in beginning of every cycle
                     icyclesinavg = icyclesinavg + 1
                  endif

                  icollectphase = 0
                  istepincycle = mod(istep,nstepsincycle)
                  if(istepincycle.eq.0) istepincycle=nstepsincycle
                  if(istepincycle.eq.iphase*nstepsbtwphase) then
                     icollectphase = 1
                     iphase = iphase+1 ! use 'iphase-1' below
                  endif

                  if(icollectphase.eq.1) then
                     tfactphase = one/icyclesinavg

                     if(myrank.eq.master) then
                       write(*,*) 'nsamples in phase ',iphase-1,': ',
     &                             icyclesinavg
                     endif

                     yphbar(:,1,iphase-1) = tfactphase*yold(:,1) +
     &                          (one-tfactphase)*yphbar(:,1,iphase-1)
                     yphbar(:,2,iphase-1) = tfactphase*yold(:,2) +
     &                          (one-tfactphase)*yphbar(:,2,iphase-1)
                     yphbar(:,3,iphase-1) = tfactphase*yold(:,3) +
     &                          (one-tfactphase)*yphbar(:,3,iphase-1)
                     yphbar(:,4,iphase-1) = tfactphase*yold(:,4) +
     &                          (one-tfactphase)*yphbar(:,4,iphase-1)
                     yphbar(:,5,iphase-1) = tfactphase*sqrt(yold(:,1)**2
     &                          +yold(:,2)**2+yold(:,3)**2) +
     &                          (one-tfactphase)*yphbar(:,5,iphase-1)
!MR CHANGE
                     yphbar(:,6,iphase-1) =
     &                              tfactphase*yold(:,1)*yold(:,1)
     &                           +(one-tfactphase)*yphbar(:,6,iphase-1)

                     yphbar(:,7,iphase-1) =
     &                              tfactphase*yold(:,1)*yold(:,2)
     &                           +(one-tfactphase)*yphbar(:,7,iphase-1)

                     yphbar(:,8,iphase-1) =
     &                              tfactphase*yold(:,1)*yold(:,3)
     &                           +(one-tfactphase)*yphbar(:,8,iphase-1)

                     yphbar(:,9,iphase-1) =
     &                              tfactphase*yold(:,2)*yold(:,2)
     &                           +(one-tfactphase)*yphbar(:,9,iphase-1)

                     yphbar(:,10,iphase-1) =
     &                              tfactphase*yold(:,2)*yold(:,3)
     &                           +(one-tfactphase)*yphbar(:,10,iphase-1)

                     yphbar(:,11,iphase-1) =
     &                              tfactphase*yold(:,3)*yold(:,3)
     &                           +(one-tfactphase)*yphbar(:,11,iphase-1)

                     if(ivort == 1) then
                       yphbar(:,12,iphase-1) =
     &                              tfactphase*vorticity(:,1)
     &                           +(one-tfactphase)*yphbar(:,12,iphase-1)
                       yphbar(:,13,iphase-1) =
     &                              tfactphase*vorticity(:,2)
     &                           +(one-tfactphase)*yphbar(:,13,iphase-1)
                       yphbar(:,14,iphase-1) =
     &                              tfactphase*vorticity(:,3)
     &                           +(one-tfactphase)*yphbar(:,14,iphase-1)
                       yphbar(:,15,iphase-1) =
     &                              tfactphase*vorticity(:,4)
     &                           +(one-tfactphase)*yphbar(:,15,iphase-1)
                    endif
!MR CHANGE END
                  endif
               endif
c
c compute rms
c
               if(icollectybar.eq.1) then
                  rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2
                  rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2
                  rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2
                  rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2
               endif
            endif
 2003       continue ! we get here if stopjob equals lstep and this jumped over
!           the statistics computation because we have no new data to average in
!           rather we are just trying to output the last state that was not already
!           written
c
c.... ---------------------->  Complete Restart  Processing  <----------------------
c
!   for now it is the same frequency but need to change this
!   soon.... but don't forget to change the field counter in
!  new_interface.cc
!
        if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or.
     &      istep.eq.nstep(itseq)) then
          if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or.
     &         (nstp .eq. 0))) then
             if(
     &          ((irscale.ge.0).or.(itwmod.gt.0) .or.
     &          ((nsonmax.eq.1).and.iLES.gt.0)))
     &          call rwvelb  ('out ',  velbar  ,ifail)
          endif

          lesId   = numeqns(1)
          if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
          if(myrank.eq.0)  then
            tcormr1 = TMRC()
          endif
          if((nsolflow.eq.1).and.(ipresPrjFlag.eq.1)) then
           call saveLesRestart( lesId,  aperm , nshg, myrank, lstep,
     &                    nPermDims )
           if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
           if(myrank.eq.0)  then
            tcormr2 = TMRC()
            write(6,*) 'call saveLesRestart for projection and'//
     &           'pressure projection vectors', tcormr2-tcormr1
           endif
          endif

          if(ierrcalc.eq.1) then
c
c.....smooth the error indicators
c
            do i=1,ierrsmooth
              call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC )
            end do
            if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            if(myrank.eq.0)  then
              tcormr1 = TMRC()
            endif
            call write_error(myrank, lstep, nshg, 10, rerr )
            if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            if(myrank.eq.0)  then
              tcormr2 = TMRC()
              write(6,*) 'Time to write the error fields to the disks',
     &            tcormr2-tcormr1
            endif
          endif ! ierrcalc

          if(ioybar.eq.1) then
            if(ivort == 1) then
              call write_field(myrank,'a','ybar',4,
     &                  ybar,'d',nshg,17,lstep)
            else
              call write_field(myrank,'a','ybar',4,
     &                ybar,'d',nshg,13,lstep)
            endif

            if(abs(itwmod).ne.1 .and. iowflux.eq.1) then
              call write_field(myrank,'a','wssbar',6,
     &             wallssVecBar,'d',nshg,3,lstep)
            endif

            if(nphasesincycle .gt. 0) then
              if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
              if(myrank.eq.0)  then
                tcormr1 = TMRC()
              endif
              do iphase=1,nphasesincycle
                if(ivort == 1) then
                 call write_phavg2(myrank,'a','phase_average',13,iphase,
     &              nphasesincycle,yphbar(:,:,iphase),'d',nshg,15,lstep)
                else
                 call write_phavg2(myrank,'a','phase_average',13,iphase,
     &              nphasesincycle,yphbar(:,:,iphase),'d',nshg,11,lstep)
                endif
              end do
              if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
              if(myrank.eq.0)  then
                tcormr2 = TMRC()
                write(6,*) 'write all phase avg to the disks = ',
     &                tcormr2-tcormr1
              endif
            endif !nphasesincyle
          endif !ioybar

          if ( ( ihessian .eq. 1 ) .and. ( numpe < 2 )  )then
            uhess = zero
            gradu = zero
            tf = zero

            do ku=1,nshg
              tf(ku,1) = x(ku,1)**3
            end do
            call hessian( yold, x,     shp,  shgl,   iBC,
     &              shpb, shglb, iper, ilwork, uhess, gradu )

            call write_hessian( uhess, gradu, nshg )
          endif

          if(iRANS.lt.0) then
            if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            if(myrank.eq.0)  then
              tcormr1 = TMRC()
            endif
            call write_field(myrank,'a','dwal',4,d2wall,'d',
     &                       nshg,1,lstep)
            if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            if(myrank.eq.0)  then
              tcormr2 = TMRC()
              write(6,*) 'Time to write dwal to the disks = ',
     &        tcormr2-tcormr1
            endif
          endif !iRANS

        endif ! write out complete restart state
        !next 2 lines are two ways to end early
        if(stopjob.eq.-2) goto 2002
        if(istop.eq.1000) goto 2002 ! stop when delta small (see rstatic)
 2000 continue
 2002 continue

! done with time stepping so deallocate fields already written
!
          if(ioybar.eq.1) then
            deallocate(ybar)
            if(abs(itwmod).ne.1 .and. iowflux.eq.1) then
              deallocate(wallssVecbar)
            endif
            if(nphasesincycle .gt. 0) then
              deallocate(yphbar)
            endif !nphasesincyle
          endif !ioybar
          if(ivort == 1) then
            deallocate(strain,vorticity)
          endif
          if(abs(itwmod).ne.1 .and. iowflux.eq.1) then
            deallocate(wallssVec)
          endif
          if(iRANS.lt.0) then
            deallocate(d2wall)
          endif


CAD         tcorecp2 = second(0)
CAD         tcorewc2 = second(-1)

CAD         write(6,*) 'T(core) cpu-wallclock = ',tcorecp2-tcorecp1,
CAD     &                                        tcorewc2-tcorewc1

         if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         if(myrank.eq.0)  then
            tcorecp2 = TMRC()
             write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1
             write(6,*) '(Elm. form.',tcorecp(1),
     &                    ',Lin. alg. sol.',tcorecp(2),')'
             write(6,*) '(Elm. form. Scal.',tcorecpscal(1),
     &                   ',Lin. alg. sol. Scal.',tcorecpscal(2),')'
             write(6,*) ''

         endif

         call print_system_stats(tcorecp, tcorecpscal)
         call print_mesh_stats()
         call print_mpi_stats()
         if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
!         return
c         call MPI_Finalize()
c         call MPI_ABORT(MPI_COMM_WORLD, ierr)

 3000 continue


c
c.... close history and aerodynamic forces files
c
      if (myrank .eq. master) then
!         close (ihist)
         close (iforce)
         close(76)
         if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then
            close (8177)
         endif
      endif
c
c.... close varts file for probes
c
      if(exts) then
        do jj=1,ntspts
          if (myrank == iv_rank(jj)) then
            close(1000+jj)
          endif
        enddo
        deallocate (ivarts)
        deallocate (ivartsg)
        deallocate (iv_rank)
        deallocate (vartssoln)
        deallocate (vartssolng)
      endif

!MR CHANGE
      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      if(myrank.eq.0)  then
          write(*,*) 'itrdrv - done with aerodynamic forces'
      endif
!MR CHANGE

      do isrf = 0,MAXSURF
!        if ( nsrflist(isrf).ne.0 ) then
        if ( nsrflist(isrf).ne.0 .and.
     &                     myrank.eq.irankfilesforce(isrf)) then
          iunit=60+isrf
          close(iunit)
        endif
      enddo

!MR CHANGE
      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      if(myrank.eq.0)  then
          write(*,*) 'itrdrv - done with MAXSURF'
      endif
!MR CHANGE


 5    format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10)
 444  format(6(2x,e14.7))
c
c.... end
c
      if(nsolflow.eq.1) then
         deallocate (lhsK)
         deallocate (lhsP)
         IF (svLSFlag .NE. 1) THEN
         deallocate (aperm)
         deallocate (atemp)
         ENDIF
      endif
      if(nsclrsol.gt.0) then
         deallocate (lhsS)
         deallocate (apermS)
         deallocate (atempS)
      endif

      if(iabc==1) deallocate(acs)

!MR CHANGE
      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      if(myrank.eq.0)  then
          write(*,*) 'itrdrv - done - BACK TO process.f'
      endif
!MR CHANGE



      return
      end

      subroutine lesmodels(y,     ac,        shgl,      shp,
     &                     iper,  ilwork,    rowp,      colm,
     &                     nsons, ifath,     x,
     &                     iBC,   BC)

      include "common.h"

      real*8    y(nshg,ndof),              ac(nshg,ndof),
     &            x(numnp,nsd),
     &            BC(nshg,ndofBC)
      real*8    shp(MAXTOP,maxsh,MAXQPT),
     &            shgl(MAXTOP,nsd,maxsh,MAXQPT)

c
      integer   rowp(nshg,nnz),         colm(nshg+1),
     &            iBC(nshg),
     &            ilwork(nlwork),
     &            iper(nshg)
      dimension ifath(numnp),    nsons(nfath)

      real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4
      real*8, allocatable, dimension(:) :: stabdis,cdelsq1
      real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3

      if( (iLES.gt.1) )   then ! Allocate Stuff for advanced LES models
         allocate (fwr2(nshg))
         allocate (fwr3(nshg))
         allocate (fwr4(nshg))
         allocate (xavegt(nfath,12))
         allocate (xavegt2(nfath,12))
         allocate (xavegt3(nfath,12))
         allocate (stabdis(nfath))
      endif

c.... get dynamic model coefficient
c
      ilesmod=iLES/10
c
c digit bit set filter rule, 10 bit set model
c
      if (ilesmod.eq.0) then    ! 0 < iLES< 10 => dyn. model calculated
                                ! at nodes based on discrete filtering


         if(isubmod.eq.2) then
            call SUPGdis(y,      ac,        shgl,      shp,
     &                   iper,   ilwork,
     &                   nsons,  ifath,     x,
     &                   iBC,    BC, stabdis, xavegt3)
         endif

         if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no
                                                     ! sub-model
                                                     ! or SUPG
                                                     ! model wanted

            if(i2filt.eq.0)then ! If simple filter

               if(modlstats .eq. 0) then ! If no model stats wanted
                  call getdmc (y,       shgl,      shp,
     &                         iper,       ilwork,    nsons,
     &                         ifath,      x)
               else             ! else get model stats
                  call stdfdmc (y,       shgl,      shp,
     &                          iper,       ilwork,    nsons,
     &                          ifath,      x)
               endif            ! end of stats if statement

            else                ! else if twice filtering

               call widefdmc(y,       shgl,      shp,
     &                       iper,       ilwork,    nsons,
     &                       ifath,      x)


            endif               ! end of simple filter if statement

         endif                  ! end of SUPG or no sub-model if statement


         if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted
            call cdelBHsq (y,       shgl,      shp,
     &                     iper,       ilwork,    nsons,
     &                     ifath,      x,         cdelsq1)
            call FiltRat (y,       shgl,      shp,
     &                    iper,       ilwork,    nsons,
     &                    ifath,      x,         cdelsq1,
     &                    fwr4,       fwr3)


            if (i2filt.eq.0) then ! If simple filter wanted
               call DFWRsfdmc(y,       shgl,      shp,
     &                        iper,       ilwork,    nsons,
     &                        ifath,      x,         fwr2, fwr3)
            else                ! else if twice filtering wanted
               call DFWRwfdmc(y,       shgl,      shp,
     &                        iper,       ilwork,    nsons,
     &                        ifath,      x,         fwr4, fwr4)
            endif               ! end of simple filter if statement

         endif                  ! end of DFWR sub-model if statement

         if( (isubmod.eq.2) )then ! If SUPG sub-model wanted
            call dmcSUPG (y,           ac,         shgl,
     &                    shp,         iper,       ilwork,
     &                    nsons,       ifath,      x,
     &                    iBC,    BC,  rowp,       colm,
     &                    xavegt2,    stabdis)
         endif

         if(idis.eq.1)then      ! If SUPG/Model dissipation wanted
            call ediss (y,        ac,      shgl,
     &                  shp,      iper,       ilwork,
     &                  nsons,    ifath,      x,
     &                  iBC,      BC,  xavegt)
         endif

      endif                     ! end of ilesmod

      if (ilesmod .eq. 1) then  ! 10 < iLES < 20 => dynamic-mixed
                                ! at nodes based on discrete filtering
         call bardmc (y,       shgl,      shp,
     &                iper,    ilwork,
     &                nsons,   ifath,     x)
      endif

      if (ilesmod .eq. 2) then  ! 20 < iLES < 30 => dynamic at quad
                                ! pts based on lumped projection filt.

         if(isubmod.eq.0)then
            call projdmc (y,       shgl,      shp,
     &                    iper,       ilwork,    x)
         else
            call cpjdmcnoi (y,      shgl,      shp,
     &                      iper,   ilwork,       x,
     &                      rowp,   colm,
     &                      iBC,    BC)
         endif

      endif

      if( (iLES.gt.1) )   then ! Deallocate Stuff for advanced LES models
         deallocate (fwr2)
         deallocate (fwr3)
         deallocate (fwr4)
         deallocate (xavegt)
         deallocate (xavegt2)
         deallocate (xavegt3)
         deallocate (stabdis)
      endif
      return
      end

c
c...initialize the coefficients for the impedance convolution
c
      subroutine CalcImpConvCoef (numISrfs, numTpoints)

      use convolImpFlow !uses flow history and impedance for convolution

      include "common.h" !for alfi

      integer numISrfs, numTpoints

      allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC
      do j=1,numTpoints+2
         ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt
         ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi)
         ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi))
         ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi
      enddo
      ConvCoef(1,2)=zero
      ConvCoef(1,3)=zero
      ConvCoef(2,3)=zero
      ConvCoef(numTpoints+1,1)=zero
      ConvCoef(numTpoints+2,2)=zero
      ConvCoef(numTpoints+2,1)=zero
c
c...calculate the coefficients for the impedance convolution
c
      allocate (ImpConvCoef(numTpoints+2,numISrfs))

c..coefficients below assume Q linear in time step, Z constant
c            do j=3,numTpoints
c                ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3)
c     &                             + ValueListImp(j,:)*ConvCoef(j,2)
c     &                             + ValueListImp(j+1,:)*ConvCoef(j,1)
c            enddo
c            ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1)
c            ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2)
c     &                       + ValueListImp(3,:)*ConvCoef(2,1)
c            ImpConvCoef(numTpoints+1,:) =
c     &           ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3)
c     &         + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2)
c            ImpConvCoef(numTpoints+2,:) =
c     &           ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3)

c..try easiest convolution Q and Z constant per time step
      do j=3,numTpoints+1
         ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints
      enddo
      ImpConvCoef(1,:) =zero
      ImpConvCoef(2,:) =zero
      ImpConvCoef(numTpoints+2,:) =
     &           ValueListImp(numTpoints+1,:)/numTpoints
c compensate for yalpha passed not y in Elmgmr()
      ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:)
     &                  - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi
      ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi
      return
      end

c
c ... update the flow rate history for the impedance convolution, filter it and write it out
c
      subroutine UpdHistConv(y,nsrfIdList,numSrfs)

      use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs
      use convolRCRFlow !brings QHistRCR, numRCRSrfs

      include "common.h"

      integer   nsrfIdList(0:MAXSURF), numSrfs
      character*20 fname1
      character*10 cname2
      character*5 cname
      real*8    y(nshg,3) !velocity at time n+1
      real*8    NewQ(0:MAXSURF) !temporary unknown for the flow rate
                        !that needs to be added to the flow history

      call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1
c
c... for imp BC: shift QHist, add new constribution, filter and write out
c
      if(numImpSrfs.gt.zero) then
         do j=1, ntimeptpT
            QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs)
         enddo
         QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs)

c
c....filter the flow rate history
c
         cutfreq = 10           !hardcoded cutting frequency of the filter
         do j=1, ntimeptpT
            QHistTry(j,:)=QHistImp(j+1,:)
         enddo
         call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq)
c.... if no filter applied then uncomment next three lines
c         do j=1, ntimeptpT
c            QHistTryF(j,:)=QHistTry(j,:)
c         enddo

c         QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters
         do j=1, ntimeptpT
            QHistImp(j+1,:)=QHistTryF(j,:)
         enddo
c
c.... write out the new history of flow rates to Qhistor.dat
c
         if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or.
     &        (istep .eq. nstep(1)))) .and.
     &        (myrank .eq. master)) then
            open(unit=816, file='Qhistor.dat',status='replace')
            write(816,*) ntimeptpT
            do j=1,ntimeptpT+1
               write(816,*) (QHistImp(j,n),n=1, numSrfs)
            enddo
            close(816)
c... write out a copy with step number to be able to restart
            fname1 = 'Qhistor'
            fname1 = trim(fname1)//trim(cname2(lstep))//'.dat'
            open(unit=8166,file=trim(fname1),status='unknown')
            write(8166,*) ntimeptpT
            do j=1,ntimeptpT+1
               write(8166,*) (QHistImp(j,n),n=1, numSrfs)
            enddo
            close(8166)
         endif
      endif

c
c... for RCR bc just add the new contribution
c
      if(numRCRSrfs.gt.zero) then
         QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs)
c
c.... write out the new history of flow rates to Qhistor.dat
c
         if ((irs .ge. 1) .and. (myrank .eq. master)) then
            if(istep.eq.1) then
               open(unit=816,file='Qhistor.dat',status='unknown')
            else
               open(unit=816,file='Qhistor.dat',position='append')
            endif
            if(istep.eq.1) then
               do j=1,lstep
                  write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run
               enddo
            endif
            write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs)
            close(816)
c... write out a copy with step number to be able to restart
            if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or.
     &           (istep .eq. nstep(1)))) .and.
     &           (myrank .eq. master)) then
               fname1 = 'Qhistor'
               fname1 = trim(fname1)//trim(cname2(lstep))//'.dat'
               open(unit=8166,file=trim(fname1),status='unknown')
               write(8166,*) lstep+1
               do j=1,lstep+1
                  write(8166,*) (QHistRCR(j,n),n=1, numSrfs)
               enddo
               close(8166)
            endif
         endif
      endif

      return
      end

c
c...calculate the time varying coefficients for the RCR convolution
c
      subroutine CalcRCRConvCoef (stepn, numSrfs)

      use convolRCRFlow !brings in ValueListRCR, dtRCR

      include "common.h" !brings alfi

      integer numSrfs, stepn

      RCRConvCoef = zero
      if (stepn .eq. 0) then
        RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) +
     &   ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:)
     &     - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:)))
        RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi
     &     + ValueListRCR(3,:)
     &     *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:))
      endif
      if (stepn .ge. 1) then
        RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi))
     &        *(1 + (1 - exp(dtRCR(:)))/dtRCR(:))
        RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi)
     &     - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:)
     &     + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:))))
        RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi
     &     + ValueListRCR(3,:)
     &     *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:))
      endif
      if (stepn .ge. 2) then
        do j=2,stepn
         RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)*
     &        exp(-dtRCR(:)*(stepn + alfi + 2 - j))*
     &        (1 - exp(dtRCR(:)))**2
        enddo
      endif

c compensate for yalpha passed not y in Elmgmr()
      RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:)
     &                  - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi
      RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi

      return
      end

c
c...calculate the time dependent H operator for the RCR convolution
c
      subroutine CalcHopRCR (timestepRCR, stepn, numSrfs)

      use convolRCRFlow !brings in HopRCR, dtRCR

      include "common.h"

      integer numSrfs, stepn
      real*8  PdistCur(0:MAXSURF), timestepRCR

      HopRCR=zero
      call RCRint(timestepRCR*(stepn + alfi),PdistCur)
      HopRCR(1:numSrfs) = RCRic(1:numSrfs)
     &     *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs)
      return
      end
c
c ... initialize the influence of the initial conditions for the RCR BC
c
      subroutine calcRCRic(y,srfIdList,numSrfs)

      use convolRCRFlow    !brings RCRic, ValueListRCR, ValuePdist

      include "common.h"

      integer   srfIdList(0:MAXSURF), numSrfs, irankCoupled
      real*8    y(nshg,4) !need velocity and pressure
      real*8    Qini(0:MAXSURF) !initial flow rate
      real*8    PdistIni(0:MAXSURF) !initial distal pressure
      real*8    Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure
      real*8    VelOnly(nshg,3), POnly(nshg)

      allocate (RCRic(0:MAXSURF))

      if(lstep.eq.0) then
         VelOnly(:,1:3)=y(:,1:3)
         call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow
         QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR
         POnly(:)=y(:,4)        ! pressure
         call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral
         POnly(:)=one           ! one to get area
         call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area
         Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs)
      else
         Qini(1:numSrfs)=QHistRCR(1,1:numSrfs)
         Pini(1:numSrfs)=zero    ! hack
      endif
      call RCRint(istep,PdistIni) !get initial distal P (use istep)
      RCRic(1:numSrfs) = Pini(1:numSrfs)
     &          - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs)
      return
      end

c.........function that integrates a scalar over a boundary
      subroutine integrScalar(scalInt,scal,srfIdList,numSrfs)

      use pvsQbi !brings ndsurf, NASC

      include "common.h"
      include "mpif.h"

      integer   srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k
      real*8    scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF)

      scalIntProc = zero
      do i = 1,nshg
        if(numSrfs.gt.zero) then
          do k = 1,numSrfs
            irankCoupled = 0
            if (srfIdList(k).eq.ndsurf(i)) then
              irankCoupled=k
              scalIntProc(irankCoupled) = scalIntProc(irankCoupled)
     &                            + NASC(i)*scal(i)
              exit
            endif
          enddo
        endif
      enddo
c
c     at this point, each scalint has its "nodes" contributions to the scalar
c     accumulated into scalIntProc. Note, because NASC is on processor this
c     will NOT be the scalar for the surface yet
c
c.... reduce integrated scalar for each surface, push on scalInt
c
        npars=MAXSURF+1
       call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)

      return
      end

      subroutine writeTimingMessage(key,iomode,timing)
      use iso_c_binding
      use phstr
      implicit none

      character(len=*) :: key
      integer :: iomode
      real*8 :: timing
      character(len=1024) :: timing_msg
      character(len=*), parameter ::
     &  streamModeString = c_char_"stream"//c_null_char,
     &  fileModeString = c_char_"disk"//c_null_char

      timing_msg = c_char_"Time to write "//c_null_char
      call phstr_appendStr(timing_msg,key)
      if ( iomode .eq. -1 ) then
        call phstr_appendStr(timing_msg, streamModeString)
      else
        call phstr_appendStr(timing_msg, fileModeString)
      endif
      call phstr_appendStr(timing_msg, c_char_' = '//c_null_char)
      call phstr_appendDbl(timing_msg, timing)
      write(6,*) trim(timing_msg)
      return
      end subroutine