vtf-logo

src/generic/Generic_CaptureDc.f90

MODULE Generic_CaptureDC

  ! ----  over-writes the flux interplotion (fxi) elements
  ! ----  with those caluclated by WENO in locations where
  ! ----  dcflag is 1

  INTERFACE CaptureDC
     MODULE PROCEDURE OneDCaptureDc
     MODULE PROCEDURE TwoDCaptureDc
     MODULE PROCEDURE ThreeDCaptureDc
  END INTERFACE


CONTAINS


  
  SUBROUTINE OneDCaptureDc(ux,vx,fx,fxi,dcflag,direction,ifilter)

    ! ----  Shared Variables
    USE mesh
    USE array_bounds
    USE method_parms

    ! ----  Shared Procedures
    USE OneDweno
    USE Generic_EvalsAndETAX

    IMPLICIT NONE

    include 'cles.i'
    
    DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
    DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
    DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ixlo:ixhi)

    INTEGER, INTENT(IN) :: dcflag(1:nx+1,1)
    INTEGER, INTENT(IN) :: direction, ifilter

    DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,1)
    
    ! ----- local eigen values and H-correction
    DOUBLE PRECISION, ALLOCATABLE :: locfxi(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: lami(:,:,:)
    DOUBLE PRECISION, ALLOCATABLE :: eta(:)

    INTEGER :: i, iret

    ! nothing to do
    if ( use_dcflag .eq. CLES_TRUE .and. SUM(dcflag(:,direction) ) .eq. 0 ) return
    
    if ( use_carbfix .eq. 0 ) then
       ALLOCATE (lami(2,3,1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'OneDCaptureDc/lami', iret)
       endif
       ALLOCATE (eta(1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'OneDCaptureDc/eta', iret)
       endif
    endif

    whichdirection: SELECT CASE (direction)
       
    CASE (1)
       
       ALLOCATE (locfxi(ncomps,ixlo:ixhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'OneDCaptureDc/locfxi', iret)
       endif

       ! ---- the x-direction
       if ( use_carbfix .eq. 1 ) then
          ALLOCATE (lami(2,3,nx+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'OneDCaptureDc/lami', iret)
          endif
          ALLOCATE (eta(nx+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'OneDCaptureDc/eta', iret)
          endif
          do i=1,nx+1
             CALL EvalsAndETAX(ux,vx,i-1,lami(:,:,i),eta(i))
          enddo
       endif
       
       locfxi = fxi(:,:,direction)

       CALL FluxInterpWENO(locfxi, fx(1,ixlo), &
               ux(1,ixlo), vx(1,ixlo), lami, eta, &
               dcflag(1,direction), ixlo, ixhi, nx)

       if ( ifilter .eq. CLES_TRUE .and. method .eq. CLES_METHOD_HYBRID ) then
          do i=1, nx+1
             if ( dcflag(i,direction) .eq. CLES_SWITCH_TCD ) then
                fxi(1:nvars,i,direction) = fxi(1:nvars,i,direction)*alpha_eta1 + locfxi(1:nvars,i)*alpha_eta2
             else
                fxi(1:nvars,i,direction) = locfxi(1:nvars,i)
             endif
          enddo
       else
          fxi(1:nvars,1:nx+1,direction) = locfxi(1:nvars,1:nx+1)
       endif
       
10     if ( use_carbfix .eq. 1 ) then
          DEALLOCATE (lami, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'OneDCaptureDc/lami', iret)
          endif
          DEALLOCATE (eta, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'OneDCaptureDc/eta', iret)
          endif

          DEALLOCATE(locfxi, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'OneDCaptureDc/locfxi', iret)
          endif

       endif

    END SELECT whichdirection

    if ( use_carbfix .eq. 0 ) then
       DEALLOCATE (lami, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'OneDCaptureDc/lami', iret)
       endif
       DEALLOCATE (eta, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'OneDCaptureDc/eta', iret)
       endif
    endif
    
  END SUBROUTINE OneDCaptureDc

  
  SUBROUTINE TwoDCaptureDc(ux,vx,fx,fxi,dcflag,direction,ifilter)

    ! ----  Shared Variables
    USE mesh
    USE array_bounds
    USE method_parms

    ! ----  Shared Procedures
    USE OneDweno
    USE Generic_EvalsAndETAX
    USE Generic_EvalsAndETAY
    
    IMPLICIT NONE
    
    include 'cles.i'

    DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
    DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi) 
    DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ixlo:ixhi,iylo:iyhi)
    
    INTEGER, INTENT(IN) :: dcflag(1:nx+1,1:ny+1,2)
    INTEGER, INTENT(IN) :: direction, ifilter

    DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2)
    
    ! ----  local copies of the data used in WENO
    ! ----  these are a fixed size: no problems with AMR
    DOUBLE PRECISION, ALLOCATABLE :: locfxi(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: locfx(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: locux(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: locvx(:,:)
    INTEGER, ALLOCATABLE :: locdc(:)
    
    ! ----- local eigenvalues and H-correction
    DOUBLE PRECISION, ALLOCATABLE :: lami(:,:,:)
    DOUBLE PRECISION, ALLOCATABLE :: eta(:)
    DOUBLE PRECISION :: dtmp
    
    INTEGER :: i, j, iret

    ! nothing to do
    if ( use_dcflag .eq. CLES_TRUE .and. SUM(dcflag(:,:,direction) ) .eq. 0 ) return
    
    if ( use_carbfix .eq. 0 ) then
       ALLOCATE (lami(2,3,1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/lami', iret)
       endif
       ALLOCATE (eta(1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/eta', iret)
       endif
    endif

    whichdirection: SELECT CASE (direction)
       
    CASE (1)

       ALLOCATE (locfxi(ncomps,ixlo:ixhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/locfxi', iret)
       endif

       ! ---- the x-direction
       if ( use_carbfix .eq. 1 ) then
          ALLOCATE (lami(2,3,nx+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'TwoDCaptureDc/lami', iret)
          endif
          ALLOCATE (eta(nx+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'TwoDCaptureDc/eta', iret)
          endif
       endif

       DO j=1,ny, 1
         
          if ( use_carbfix .eq. 1 ) then
             do i=1,nx+1
                CALL EvalsAndETAX(ux,vx,i-1,j,lami(:,:,i),eta(i))
             enddo
          endif

          locfxi = fxi(:,:,j,direction)

          CALL FluxInterpWENO(locfxi, fx(1,ixlo,j), &
               ux(1,ixlo,j), vx(1,ixlo,j), lami, eta, &
               dcflag(1,j,direction), ixlo, ixhi, nx)

          if ( ifilter .eq. CLES_TRUE .and. method .eq. CLES_METHOD_HYBRID ) then
             do i=1, nx+1
                if ( dcflag(i,j,direction) .eq. CLES_SWITCH_TCD ) then
                   fxi(1:nvars,i,j,direction) = fxi(1:nvars,i,j,direction)*alpha_eta1 + locfxi(1:nvars,i)*alpha_eta2
                else
                   fxi(1:nvars,i,j,direction) = locfxi(1:nvars,i)
                endif
             enddo
          else
             fxi(1:nvars,1:nx+1,j,direction) = locfxi(1:nvars,1:nx+1)
          endif
          
       END DO

10     if ( use_carbfix .eq. 1 ) then
          DEALLOCATE (lami, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'TwoDCaptureDc/lami', iret)
          endif
          DEALLOCATE (eta, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'TwoDCaptureDc/eta', iret)
          endif
       endif

       DEALLOCATE(locfxi, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/locfxi', iret)
       endif
       
    CASE (2)
       
       ALLOCATE (locfxi(ncomps,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/locfxi', iret)
       endif
       ALLOCATE (locfx(nvars,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/locfx', iret)
       endif
       ALLOCATE (locux(ncomps,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/locux', iret)
       endif
       ALLOCATE (locvx(nvars,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/locvx', iret)
       endif
       ALLOCATE (locdc(1:ny+1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'TwoDCaptureDc/locdc', iret)
       endif

       if ( use_carbfix .eq. 1 ) then
          ALLOCATE (lami(2,3,ny+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'TwoDCaptureDc/lami', iret)
          endif
          ALLOCATE (eta(ny+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'TwoDCaptureDc/eta', iret)
          endif
       endif

       ! ---- the y-direction
       DO i=1, nx, 1
          
          if ( use_carbfix .eq. 1 ) then
             do j=1,ny+1
                CALL EvalsAndETAY(ux,vx,i,j-1,lami(:,:,j),eta(j))
             enddo
          endif

          locfxi(:,:) = fxi(:,i,:,direction)
          locfx(:,:) = fx(:,i,:)
          locux(:,:) = ux(:,i,:)
          locvx(:,:) = vx(:,i,:)
          locdc(:) = dcflag(i,:,direction)

          ! Swap u for v in fxi, fx, ux, vx
          do j=iylo, iyhi
             dtmp = locfxi(2,j)
             locfxi(2,j) = locfxi(3,j)
             locfxi(3,j) = dtmp
          enddo
          do j=iylo, iyhi
             dtmp = locfx(2,j)
             locfx(2,j) = locfx(3,j)
             locfx(3,j) = dtmp
          enddo
          do j=iylo, iyhi
             dtmp = locux(2,j)
             locux(2,j) = locux(3,j)
             locux(3,j) = dtmp
          enddo
          do j=iylo, iyhi
             dtmp = locvx(2,j)
             locvx(2,j) = locvx(3,j)
             locvx(3,j) = dtmp
          enddo

          ! ---- Do the WENO interpolation at this point
          CALL FluxInterpWENO(locfxi, locfx, locux, locvx, lami, eta, &
               locdc, iylo, iyhi, ny)

          do j=iylo, iyhi
             dtmp = locfxi(2,j)
             locfxi(2,j) = locfxi(3,j)
             locfxi(3,j) = dtmp
          enddo

          if ( ifilter .eq. CLES_TRUE .and. method .eq. CLES_METHOD_HYBRID ) then
             do j=1, ny+1
                if ( locdc(j) .eq. CLES_SWITCH_TCD ) then
                   fxi(1:nvars,i,j,direction) = fxi(1:nvars,i,j,direction)*alpha_eta1 + locfxi(1:nvars,j)*alpha_eta2
                else
                   fxi(1:nvars,i,j,direction) = locfxi(1:nvars,j)
                endif
             enddo
          else
             fxi(1:nvars,i,1:ny+1,direction) = locfxi(1:nvars,1:ny+1)
          endif

       ENDDO

20     DEALLOCATE(locfxi, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/locfxi', iret)
       endif
       DEALLOCATE(locfx, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/locfx', iret)
       endif
       DEALLOCATE(locux, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/locux', iret)
       endif
       DEALLOCATE(locvx, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/locvx', iret)
       endif
       DEALLOCATE(locdc, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/locdc', iret)
       endif

       if ( use_carbfix .eq. 1 ) then
          DEALLOCATE (lami, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'TwoDCaptureDc/lami', iret)
          endif
          DEALLOCATE (eta, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'TwoDCaptureDc/eta', iret)
          endif
       endif
       
    END SELECT whichdirection
  
30  if ( use_carbfix .eq. 0 ) then
       DEALLOCATE (lami, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/lami', iret)
       endif
       DEALLOCATE (eta, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'TwoDCaptureDc/eta', iret)
       endif
    endif

  END SUBROUTINE TwoDCaptureDc
  
  
  SUBROUTINE ThreeDCaptureDc(ux,vx,fx,fxi,dcflag,direction,ifilter)

    ! ----  Shared Variables
    USE mesh
    USE array_bounds
    USE method_parms
     
    ! ----  Shared Procedures
    USE OneDweno
    USE Generic_EvalsAndETAX
    USE Generic_EvalsAndETAY
    USE Generic_EvalsAndETAZ

    IMPLICIT NONE
    
    include 'cles.i'

    DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
    DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi) 
    DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)

    INTEGER, INTENT(IN) :: dcflag(1:nx+1,1:ny+1,1:nz+1,3)
    INTEGER, INTENT(IN) :: direction, ifilter

    DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3)

    ! ----  local copies of the data used in WENO
    ! ----  these are a fixed size: no problems with AMR
    DOUBLE PRECISION, ALLOCATABLE :: locfxi(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: locfx(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: locux(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: locvx(:,:)
    INTEGER, ALLOCATABLE :: locdc(:)
    
    ! ----- local eigenvalues and H-correction
    DOUBLE PRECISION, ALLOCATABLE :: lami(:,:,:)
    DOUBLE PRECISION, ALLOCATABLE :: eta(:)
    DOUBLE PRECISION :: dtmp

    INTEGER :: i,j,k, iret

    call cleslog_log_enter('ThreeDCaptureDc')

    ! nothing to do
    if ( use_dcflag .eq. CLES_TRUE .and. SUM(dcflag(:,:,:,direction) ) .eq. 0 ) goto 11

    if ( use_carbfix .eq. 0 ) then
       ALLOCATE (lami(2,3,1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/lami', iret)
       endif
       ALLOCATE (eta(1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/eta', iret)
       endif
    endif

    whichdirection: SELECT CASE (direction)
       
    CASE (1)
       
       ALLOCATE (locfxi(ncomps,ixlo:ixhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locfxi', iret)
       endif

       ! ---- the x-direction
       if ( use_carbfix .eq. 1 ) then
          ALLOCATE (lami(2,3,nx+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'ThreeDCaptureDc/lami', iret)
          endif
          ALLOCATE (eta(nx+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'ThreeDCaptureDc/eta', iret)
          endif
       endif
       
       DO k = 1,nz,1
          DO j = 1,ny, 1
             
             if ( use_carbfix .eq. 1 ) then
                do i=1,nx+1
                   CALL EvalsAndETAX(ux,vx,i-1,j,k,lami(:,:,i), eta(i))
                enddo
             endif

             locfxi = fxi(:,:,j,k,direction)

             CALL FluxInterpWENO(locfxi, fx(1,ixlo,j,k), &
                  ux(1,ixlo,j,k), vx(1,ixlo,j,k), lami, eta, &
                  dcflag(1,j,k,direction), ixlo, ixhi, nx)
             
             if ( ifilter .eq. CLES_TRUE .and. method .eq. CLES_METHOD_HYBRID ) then
                do i=1, nx+1
                   if ( dcflag(i,j,k,direction) .eq. CLES_SWITCH_TCD ) then
                      fxi(1:nvars,i,j,k,direction) = fxi(1:nvars,i,j,k,direction)*alpha_eta1 + locfxi(1:nvars,i)*alpha_eta2
                   else
                      fxi(1:nvars,i,j,k,direction) = locfxi(1:nvars,i)
                   endif
                enddo
             else
                fxi(1:nvars,1:nx+1,j,k,direction) = locfxi(1:nvars,1:nx+1)
             endif

          END DO
       END DO

10     if ( use_carbfix .eq. 1 ) then
          DEALLOCATE (lami, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'ThreeDCaptureDc/lami', iret)
          endif
          DEALLOCATE (eta, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'ThreeDCaptureDc/eta', iret)
          endif
       endif

       DEALLOCATE(locfxi, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locfxi', iret)
       endif
       
    CASE (2)

       ALLOCATE (locfxi(ncomps,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locfxi', iret)
       endif
       ALLOCATE (locfx(nvars,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locfx', iret)
       endif
       ALLOCATE (locux(ncomps,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locux', iret)
       endif
       ALLOCATE (locvx(nvars,iylo:iyhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locvx', iret)
       endif
       ALLOCATE (locdc(1:ny+1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locdc', iret)
       endif

       if ( use_carbfix .eq. 1 ) then
          ALLOCATE (lami(2,3,ny+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'ThreeDCaptureDc/lami', iret)
          endif
          ALLOCATE (eta(ny+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'ThreeDCaptureDc/eta', iret)
          endif
       endif

       ! ---- the y-direction
       DO k=1, nz, 1
          DO i=1, nx, 1

             if ( use_carbfix .eq. 1 ) then
                do j=1,ny+1
                   CALL EvalsAndETAY(ux,vx,i,j-1,k,lami(:,:,j), eta(j))
                enddo
             endif

             locfxi(:,:) = fxi(:,i,:,k,direction)
             locfx(:,:) = fx(:,i,:,k)
             locux(:,:) = ux(:,i,:,k)
             locvx(:,:) = vx(:,i,:,k)
             locdc(:) = dcflag(i,:,k,direction)

             ! Swap u for v in fx, ux, vx
             do j=iylo, iyhi
                dtmp = locfxi(2,j)
                locfxi(2,j) = locfxi(3,j)
                locfxi(3,j) = dtmp
             enddo
             do j=iylo, iyhi
                dtmp = locfx(2,j)
                locfx(2,j) = locfx(3,j)
                locfx(3,j) = dtmp
             enddo
             do j=iylo, iyhi
                dtmp = locux(2,j)
                locux(2,j) = locux(3,j)
                locux(3,j) = dtmp
             enddo
             do j=iylo, iyhi
                dtmp = locvx(2,j)
                locvx(2,j) = locvx(3,j)
                locvx(3,j) = dtmp
             enddo
             
             ! ---- Do the WENO interpolation at this point
             CALL FluxInterpWENO(locfxi, locfx, locux, locvx, lami, eta,&
                  locdc, iylo, iyhi, ny)
             
             do j=iylo, iyhi
                dtmp = locfxi(2,j)
                locfxi(2,j) = locfxi(3,j)
                locfxi(3,j) = dtmp
             enddo
             
             if ( ifilter .eq. CLES_TRUE .and. method .eq. CLES_METHOD_HYBRID ) then
                do j=1, ny+1
                   if ( locdc(j) .eq. CLES_SWITCH_TCD ) then
                      fxi(1:nvars,i,j,k,direction) = fxi(1:nvars,i,j,k,direction)*alpha_eta1 + locfxi(1:nvars,j)*alpha_eta2
                   else
                      fxi(1:nvars,i,j,k,direction) = locfxi(1:nvars,j)
                   endif
                enddo
             else
                fxi(1:nvars,i,1:ny+1,k,direction) = locfxi(1:nvars,1:ny+1)
             endif
             
          END DO
       END DO 

20     DEALLOCATE(locfxi, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locfxi', iret)
       endif
       DEALLOCATE(locfx, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locfx', iret)
       endif
       DEALLOCATE(locux, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locux', iret)
       endif
       DEALLOCATE(locvx, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locvx', iret)
       endif
       DEALLOCATE(locdc, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locdc', iret)
       endif

       if ( use_carbfix .eq. 1 ) then
          DEALLOCATE (lami, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'ThreeDCaptureDc/lami', iret)
          endif
          DEALLOCATE (eta, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'ThreeDCaptureDc/eta', iret)
          endif
       endif

    CASE (3)

       ALLOCATE (locfxi(ncomps,izlo:izhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locfxi', iret)
       endif
       ALLOCATE (locfx(nvars,izlo:izhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locfx', iret)
       endif
       ALLOCATE (locux(ncomps,izlo:izhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locux', iret)
       endif
       ALLOCATE (locvx(nvars,izlo:izhi), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locvx', iret)
       endif
       ALLOCATE (locdc(1:nz+1), stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
               'ThreeDCaptureDc/locdc', iret)
       endif

       if ( use_carbfix .eq. 1 ) then
          ALLOCATE (lami(2,3,nz+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'ThreeDCaptureDc/lami', iret)
          endif
          ALLOCATE (eta(nz+1), stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
                  'ThreeDCaptureDc/eta', iret)
          endif
       endif

       ! ---- the z-direction
       DO j = 1, ny, 1   
          DO i= 1, nx, 1
            
             if ( use_carbfix .eq. 1 ) then
                do k=1,nz+1
                   CALL EvalsAndETAZ(ux,vx,i,j,k-1,lami(:,:,k), eta(k))
                enddo
             endif

             locfxi(:,:) = fxi(:,i,j,:,direction)
             locfx(:,:) = fx(:,i,j,:)
             locux(:,:) = ux(:,i,j,:)
             locvx(:,:) = vx(:,i,j,:)
             locdc(:) = dcflag(i,j,:,direction)
             
             ! Swap u for v in fx, ux, vx
             do k=izlo, izhi
                dtmp = locfxi(2,k)
                locfxi(2,k) = locfxi(4,k)
                locfxi(4,k) = dtmp
             enddo
             do k=izlo, izhi
                dtmp = locfx(2,k)
                locfx(2,k) = locfx(4,k)
                locfx(4,k) = dtmp
             enddo
             do k=izlo, izhi
                dtmp = locux(2,k)
                locux(2,k) = locux(4,k)
                locux(4,k) = dtmp
             enddo
             do k=izlo, izhi
                dtmp = locvx(2,k)
                locvx(2,k) = locvx(4,k)
                locvx(4,k) = dtmp
             enddo

             ! ---- Do the WENO interpolation at this point
             CALL FluxInterpWENO(locfxi, locfx, locux, locvx, lami, eta,&
                  locdc, izlo, izhi, nz)

             do k=izlo, izhi
                dtmp = locfxi(2,k)
                locfxi(2,k) = locfxi(4,k)
                locfxi(4,k) = dtmp
             enddo
             
             if ( ifilter .eq. CLES_TRUE .and. method .eq. CLES_METHOD_HYBRID ) then
                do k=1, nz+1
                   if ( locdc(k) .eq. CLES_SWITCH_TCD ) then
                      fxi(1:nvars,i,j,k,direction) = fxi(1:nvars,i,j,k,direction)*alpha_eta1 + locfxi(1:nvars,k)*alpha_eta2
                   else
                      fxi(1:nvars,i,j,k,direction) = locfxi(1:nvars,k)
                   endif
                enddo
             else
                fxi(1:nvars,i,j,1:nz+1,direction) = locfxi(1:nvars,1:nz+1)
             endif

          END DO
       END DO

30     DEALLOCATE(locfxi, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locfxi', iret)
       endif
       DEALLOCATE(locfx, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locfx', iret)
       endif
       DEALLOCATE(locux, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locux', iret)
       endif
       DEALLOCATE(locvx, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locvx', iret)
       endif
       DEALLOCATE(locdc, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/locdc', iret)
       endif

       if ( use_carbfix .eq. 1 ) then
          DEALLOCATE (lami, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'ThreeDCaptureDc/lami', iret)
          endif
          DEALLOCATE (eta, stat=iret)
          if ( iret .ne. 0 ) then
             call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
                  'ThreeDCaptureDc/eta', iret)
          endif
       endif

    END SELECT whichdirection

40  if ( use_carbfix .eq. 0 ) then
       DEALLOCATE (lami, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/lami', iret)
       endif
       DEALLOCATE (eta, stat=iret)
       if ( iret .ne. 0 ) then
          call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
               'ThreeDCaptureDc/eta', iret)
       endif
    endif

11    call cleslog_log_exit('ThreeDCaptureDc')
  
  END SUBROUTINE ThreeDCaptureDc

END MODULE Generic_CaptureDC

<