RADI_FlimDiffFlux.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
!******************************************************************************
!
! Purpose: compute FLDTRAN diffusion flux: c*lamda(Er)/Kr*d_j(Er)
!
! Description: this routine compute d_j(Er), while the parameters defining
!              the diffusion coefficient are known from previous computations.
!
! Input: region  = data of current region
!
! Output: region%levels%radi%diss = diffusion flux added to FLD dissipation.
!
! Notes: none.
!
!******************************************************************************
!
! $Id: RADI_FlimDiffFlux.F90,v 1.3 2008/12/06 08:44:37 mtcampbe Exp $
!
! Copyright: (c) 2003 by the University of Illinois
!
!******************************************************************************

SUBROUTINE radi_flimdiffflux( region )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
#ifdef RFLO
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
                            rflo_getnodeoffset, rflo_calcgradvector

#include "Indexing.h"
#endif
#ifdef RFLU
  USE rflu_moddifferentiation, ONLY: rflu_computegradfacesbak, & 
                                     rflu_computegradfacespatchesbak
#endif
  USE radi_modinterfaces, ONLY : radi_flimdifffluxpatch
  USE moderror
  USE modparameters
  USE radi_modparameters
  IMPLICIT NONE

! ... parameters
#ifdef RFLO
  TYPE(t_region) :: region
#endif
#ifdef RFLU
  TYPE(t_region), POINTER :: region
#endif

! ... loop variables
  INTEGER :: i, j, k, ic, ipatch

! ... local variables
  INTEGER :: ibegv, iendv, ibegg, iendg, ijkn, ijkc0, ijkc1
  REAL(RFREAL) :: beta, sounda, flima, coefa, diffcoef, modsf, fd
  REAL(RFREAL) :: sface(3), radenx, radeny, radenz
  REAL(RFREAL), POINTER :: dv(:,:), rcv(:,:), rdiss(:,:), flim(:), coef(:,:)
  REAL(RFREAL), POINTER :: qr(:)

#ifdef RFLO
  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ilev, icoff, ijcoff, inoff, ijnoff
  REAL(RFREAL), POINTER :: avgco(:,:), sf(:,:), grad(:,:)
#endif
#ifdef RFLU
  INTEGER, POINTER      :: f2c(:,:)
  REAL(RFREAL), POINTER :: fn(:,:), grad(:,:,:)
#endif

!******************************************************************************

  CALL registerfunction( region%global,'RADI_FlimDiffFlux',&
  'RADI_FlimDiffFlux.F90' )

! get dimensions and pointers ------------------------------------------------

#ifdef RFLO
  ilev  =  region%currLevel
  dv    => region%levels(ilev)%mixt%dv
  rcv   => region%levels(ilev)%radi%cv
  rdiss => region%levels(ilev)%radi%diss
  flim  => region%levels(ilev)%radi%fluxLim
  coef  => region%levels(ilev)%radi%radCoef
  
  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
#endif
#ifdef RFLU
  dv    => region%mixt%dv
  rcv   => region%radi%cv
  rdiss => region%radi%diss  
  flim  => region%radi%fluxLim
  coef  => region%radi%radCoef
#endif
  ibegv = cv_radi_ener
  iendv = cv_radi_ener
  ibegg = gr_radi_ex
  iendg = gr_radi_ez

! get needed quantities

  beta  = region%mixtInput%betrk(region%irkStep)

! get gradients of radiation Energy Er

#ifdef RFLO
  CALL rflo_calcgradvector( region,ibegv,iendv,ibegg,iendg, &
                                   region%levels(ilev)%radi%cv, &
                                   region%levels(ilev)%radi%gradi, &
                                   region%levels(ilev)%radi%gradj, &
                                   region%levels(ilev)%radi%gradk )
#endif
#ifdef RFLU
  CALL rflu_computegradfacesbak( region,ibegv,iendv,ibegg,iendg, &
                                 region%radi%cv,region%radi%gradi )
  CALL rflu_computegradfacespatchesbak( region,ibegv,iendv,ibegg,iendg, &
                                        region%radi%cv,region%radi%bGradi )
#endif

! interior fluxes -------------------------------------------------------------

  CALL computeflux( icoord )
#ifdef RFLO
  CALL computeflux( jcoord )
  CALL computeflux( kcoord )
#endif

! fluxes through boundaries ---------------------------------------------------
#ifdef RFLO  
  DO ipatch=1,region%nPatches
    CALL radi_flimdifffluxpatch( region,region%levels(ilev)%patches(ipatch) )
  ENDDO
#endif
#ifdef RFLU
  DO ipatch = 1,region%grid%nPatches
    CALL radi_flimdifffluxpatch( region,region%patches(ipatch) )
  END DO ! ipatch
#endif

! finalize --------------------------------------------------------------------

  CALL deregisterfunction( region%global )

! =============================================================================
!   Flux computation subroutines
! =============================================================================

CONTAINS

  SUBROUTINE computeflux( ijk )

! ... parameters
    INTEGER   :: ijk

! ... local variables
    INTEGER   :: ibeg,iend,jbeg,jend,kbeg,kend, iadd,jadd,kadd
    REAL(RFREAL) :: ac0, ac1

! - Set limits and pointers ---------------------------------------------------

#ifdef RFLO
    IF (ijk==icoord) THEN
      ibeg = ipcbeg+1
      iend = ipcend
      jbeg = jpcbeg
      jend = jpcend
      kbeg = kpcbeg
      kend = kpcend
      iadd = -1
      jadd = 0
      kadd = 0
      qr    => region%levels(ilev)%radi%qri
      grad  => region%levels(ilev)%radi%gradi
      sf    => region%levels(ilev)%grid%si
      avgco => region%levels(ilev)%grid%c2fCoI
    ELSEIF (ijk==jcoord) THEN
      ibeg = ipcbeg
      iend = ipcend
      jbeg = jpcbeg+1
      jend = jpcend
      kbeg = kpcbeg
      kend = kpcend
      iadd = 0
      jadd = -1
      kadd = 0
      qr    => region%levels(ilev)%radi%qrj
      grad  => region%levels(ilev)%radi%gradj
      sf    => region%levels(ilev)%grid%sj
      avgco => region%levels(ilev)%grid%c2fCoJ
    ELSEIF (ijk==kcoord) THEN
      ibeg = ipcbeg
      iend = ipcend
      jbeg = jpcbeg
      jend = jpcend
      kbeg = kpcbeg+1
      kend = kpcend
      iadd = 0
      jadd = 0
      kadd = -1
      qr    => region%levels(ilev)%radi%qrk
      grad  => region%levels(ilev)%radi%gradk
      sf    => region%levels(ilev)%grid%sk
      avgco => region%levels(ilev)%grid%c2fCoK
    ENDIF
#endif
#ifdef RFLU
    ibeg =  1
    iend =  region%grid%nFaces
    f2c  => region%grid%f2c
    qr   => region%radi%qri
    grad => region%radi%gradi
    fn   => region%grid%fn
#endif

! -- flux in ijk-direction (except through boundary) --------------------------

#ifdef RFLO
    DO k=kbeg,kend
      DO j=jbeg,jend
        DO i=ibeg,iend

          ijkc0 = indijk(i     ,j     ,k     ,icoff,ijcoff)
          ijkc1 = indijk(i+iadd,j+jadd,k+kadd,icoff,ijcoff)
          ijkn  = indijk(i     ,j     ,k     ,inoff,ijnoff)
          ac0   = avgco(2,ijkn)
          ac1   = avgco(1,ijkn)
          sface(1)= sf(xcoord,ijkn)
          sface(2)= sf(ycoord,ijkn)
          sface(3)= sf(zcoord,ijkn)
#endif
#ifdef RFLU
    ac0 = 0.5_rfreal
    ac1 = 0.5_rfreal
    DO ijkn = ibeg,iend
          ijkc0 = f2c(1,ijkn)
          ijkc1 = f2c(2,ijkn)
          sface(1)= fn(xcoord,ijkn)*fn(xyzmag,ijkn)
          sface(2)= fn(ycoord,ijkn)*fn(xyzmag,ijkn)
          sface(3)= fn(zcoord,ijkn)*fn(xyzmag,ijkn)           
#endif
          sounda   = ac0*dv(dv_mixt_soun ,ijkc0)+ac1*dv(dv_mixt_soun ,ijkc1)
          flima    = ac0*flim(ijkc0)+ac1*flim(ijkc1)
          coefa    = ac0*coef(ijkc0,radi_coeff_extinct)+ &
                     ac1*coef(ijkc1,radi_coeff_extinct)
          diffcoef = sounda*flima/coefa 

#ifdef RFLO
          radenx = grad(gr_radi_ex,ijkn)
          radeny = grad(gr_radi_ey,ijkn)
          radenz = grad(gr_radi_ez,ijkn)
#endif
#ifdef RFLU
          radenx = grad(xcoord,gr_radi_ex,ijkn)
          radeny = grad(ycoord,gr_radi_ex,ijkn)
          radenz = grad(zcoord,gr_radi_ex,ijkn)
#endif

          fd = diffcoef* &
               (radenx*sface(1)+radeny*sface(2)+radenz*sface(3))

          rdiss(cv_radi_ener,ijkc0) = rdiss(cv_radi_ener,ijkc0) + fd*beta
          rdiss(cv_radi_ener,ijkc1) = rdiss(cv_radi_ener,ijkc1) - fd*beta

! ------- store radiant flux normal to cell face in positive direction in qr
          modsf     = sqrt( sface(1)*sface(1) + &
                            sface(2)*sface(2) + &
                            sface(3)*sface(3))
          qr(ijkn) = fd/modsf
#ifdef RFLO
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k
#else
    ENDDO       ! ijkN
#endif

  END SUBROUTINE computeflux

END SUBROUTINE radi_flimdiffflux

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: RADI_FlimDiffFlux.F90,v $
! Revision 1.3  2008/12/06 08:44:37  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:49  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/09/30 17:49:10  wasistho
! prepared for full FLD radiation model
!
!
!
!******************************************************************************