masscorr Subroutine

  subroutine masscorr
    !> correct the velocities to get prescribed flow rate

    use modglobal, only : ib,ie,jb,je,ih,jh,kb,ke,kh,dzf,dxf,dy,zh,dt,rk3step,&
                          uflowrate,vflowrate,linoutflow,&
                          luoutflowr,lvoutflowr,luvolflowr,lvvolflowr
    use modfields, only : um,up,vm,vp,uouttot,udef,vouttot,vdef,&
                          uoutarea,voutarea,fluidvol,IIu,IIv,IIus,IIvs
    use modmpi,    only : myid,comm3d,mpierr,nprocs,MY_REAL,sumy_ibm,sumx_ibm,avexy_ibm

    real, dimension(kb:ke+kh)     :: uvol
    real, dimension(kb:ke+kh)     :: vvol
    real, dimension(kb:ke+kh)     :: uvolold
    real, dimension(kb:ke+kh)     :: vvolold
    real, dimension(kb:ke)        :: uout
    real, dimension(kb:ke)        :: vout
    real, dimension(kb:ke)        :: uoutold
    real, dimension(kb:ke)        :: voutold

    real                          rk3coef,rk3coefi,&
                                  uoutflow,voutflow,&
                                  uflowrateold,vflowrateold
    integer                       i,j,k

    if ((.not.linoutflow) .and. (luoutflowr)) then
      rk3coef = dt / (4. - dble(rk3step))
      rk3coefi = 1 / rk3coef

      ! Assumes ie=itot
      udef = 0.
      uout = 0.
      uoutflow = 0.
      uoutold = 0.

      ! integrate u fixed at outlet ie along y
      call sumy_ibm(uout,up(ie,jb:je,kb:ke)*dy,ie,ie,jb,je,kb,ke,IIu(ie,jb:je,kb:ke))  ! u tendency at previous time step
      call sumy_ibm(uoutold,um(ie,jb:je,kb:ke)*dy,ie,ie,jb,je,kb,ke,IIu(ie,jb:je,kb:ke))  ! u at previous time step

      ! integrate u in z
      do k=kb,ke
         uout(k) = rk3coef*uout(k)*dzf(k)
         uoutold(k) = uoutold(k)*dzf(k)
      end do

      uoutflow = sum(uout(kb:ke))
      uflowrateold = sum(uoutold(kb:ke))

      ! average over outflow area
      uoutflow = uoutflow/uoutarea
      uflowrateold = uflowrateold/uoutarea

      ! flow correction to match outflow rate
      udef = uflowrate - (uoutflow + uflowrateold)

      do k = kb,ke
        do j = jb,je
            do i = ib,ie
              up(i,j,k) = up(i,j,k) + udef*rk3coefi
            end do
        end do
      end do

      ! bss116 calculate uouttot which is used in modboundary.
      ! this really should be in the routine directly!
      uouttot = sum(uout(kb:ke))  ! mass flow rate at outlet

    elseif ((.not.linoutflow) .and. (luvolflowr)) then
      rk3coef = dt / (4. - dble(rk3step))
      rk3coefi = 1 / rk3coef

      udef = 0.
      uoutflow = 0.
      uvol = 0.
      uvolold = 0.

      ! Assumes equidistant grid
      call avexy_ibm(uvol(kb:ke+kh),up(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIu(ib:ie,jb:je,kb:ke+kh),IIus(kb:ke+kh),.false.)
      call avexy_ibm(uvolold(kb:ke+kh),um(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIu(ib:ie,jb:je,kb:ke+kh),IIus(kb:ke+kh),.false.)

      ! average over fluid volume
      uoutflow = rk3coef*sum(uvol(kb:ke)*dzf(kb:ke)) / zh(ke+1)
      uflowrateold =  sum(uvolold(kb:ke)*dzf(kb:ke)) / zh(ke+1)

      ! flow correction to match outflow rate
      udef = uflowrate - (uoutflow + uflowrateold)

      do k = kb,ke
        do j = jb,je
            do i = ib,ie
              up(i,j,k) = up(i,j,k) + udef*rk3coefi
            end do
        end do
      end do

    end if

    if ((.not.linoutflow) .and. (lvoutflowr)) then
      rk3coef = dt / (4. - dble(rk3step))
      rk3coefi = 1 / rk3coef

      ! Assumes je=jtot
      vdef = 0.
      vout = 0.
      voutflow = 0.
      voutold = 0.

      ! integrate v fixed at outlet je along x
      ! if (myid==nprocs-1) then
      !    do k=kb,ke
      !       vout(k) = sum(vp(ib:ie,je,k)*IIv(ib:ie,je,k)*dxf(ib:ie))  ! v tendency at previous time step
      !       voutold(k) = sum(vm(ib:ie,je,k)*IIv(ib:ie,je,k)*dxf(ib:ie))  ! v at previous time step
      !    end do
      ! end if
      call sumy_ibm(vout,vp(ib:je,je,kb:ke)*dxf(1),ib,ie,je,je,kb,ke,IIv(ib:ie,je,kb:ke))  ! v tendency at previous time step
      call sumy_ibm(voutold,vm(ib:ie,je,kb:ke)*dxf(1),ib,ie,je,je,kb,ke,IIv(ib:ie,je,kb:ke))  ! v at previous time step

      ! integrate v in z
      do k=kb,ke
         vout(k) = rk3coef*vout(k)*dzf(k)
         voutold(k) = voutold(k)*dzf(k)
      end do

      voutflow = sum(vout(kb:ke))
      vflowrateold = sum(voutold(kb:ke))

      ! average over outflow area
      voutflow = voutflow/voutarea
      vflowrateold = vflowrateold/voutarea

      ! flow correction to match outflow rate
      vdef = vflowrate - (voutflow + vflowrateold)

      do k = kb,ke
        do j = jb,je
            do i = ib,ie
              vp(i,j,k) = vp(i,j,k) + vdef*rk3coefi
            end do
        end do
      end do

    elseif ((.not.linoutflow) .and. (lvvolflowr)) then
      rk3coef = dt / (4. - dble(rk3step))
      rk3coefi = 1 / rk3coef

      vdef = 0.
      voutflow = 0.
      vvol = 0.
      vvolold = 0.

      ! Assumes equidistant grid
      call avexy_ibm(vvol(kb:ke+kh),vp(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIv(ib:ie,jb:je,kb:ke+kh),IIvs(kb:ke+kh),.false.)
      call avexy_ibm(vvolold(kb:ke+kh),vm(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIv(ib:ie,jb:je,kb:ke+kh),IIvs(kb:ke+kh),.false.)

      ! average over fluid volume
      voutflow = rk3coef*sum(vvol(kb:ke)*dzf(kb:ke)) / zh(ke+1)
      vflowrateold =  sum(vvolold(kb:ke)*dzf(kb:ke)) / zh(ke+1)

      ! flow correction to match outflow rate
      vdef = vflowrate - (voutflow + vflowrateold)
      do k = kb,ke
        do j = jb,je
            do i = ib,ie
              vp(i,j,k) = vp(i,j,k) + vdef*rk3coefi
            end do
        end do
      end do

    end if

  end subroutine masscorr