1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE SCHEME_DMFS(neighs,nfcelt,mu,eta,ksi,w,pathway,Lb,Gi,G,&
   4 |      &                V,k_scat,k_abs,Qr,Li,H,S,norm,alpha,ndir,ncells,  &
   5 |      &                nfacemax,nkabs,dir_d,dir_f,bcell,bface,nbface)
   6 | 
   7 |         USE mod_pmm
   8 | 
   9 |         IMPLICIT NONE
  10 | 
  11 | !       IN
  12 |         DOM_INT           :: ndir, ncells, nfacemax, nkabs
  13 |         DOM_INT           :: dir_d, dir_f, ierr, nbface
  14 |         DOM_INT, DIMENSION (2*nfacemax,ncells):: neighs
  15 |         DOM_INT, DIMENSION (ncells)           :: nfcelt
  16 |         DOM_INT, DIMENSION (ncells,ndir)      :: pathway
  17 | 
  18 |         DOM_REAL          :: alpha
  19 |         DOM_REAL,DIMENSION(ndir)              :: mu, eta, ksi, w
  20 |         DOM_REAL,DIMENSION(ncells)            :: Lb, V, k_scat
  21 |         DOM_REAL,DIMENSION(ncells)            :: k_abs
  22 |         DOM_REAL,DIMENSION(ncells)            :: Gi
  23 |         DOM_REAL,DIMENSION(nfacemax,ncells)   :: S
  24 |         DOM_INT ,DIMENSION(nbface)            :: bcell, bface
  25 |         DOM_REAL,DIMENSION(3,nfacemax,ncells) :: norm
  26 | 
  27 | !       OUT
  28 |         DOM_REAL,DIMENSION(ncells)            :: G
  29 |         DOM_REAL,DIMENSION(3,ncells)          :: Qr
  30 |         DOM_REAL,DIMENSION(nfacemax,ncells)   :: Li
  31 |         DOM_REAL,DIMENSION(nfacemax,ncells)   :: H
  32 | 
  33 | !       LOCAL
  34 |         DOM_INT           :: idir,icell,i,j,n,k
  35 |         DOM_INT           :: neighfaceid, neighcellid
  36 |         DOM_REAL          :: nmin
  37 |         DOM_REAL          :: Lpi
  38 |         DOM_REAL,DIMENSION(nfacemax)          :: Dm_ij
  39 |         DOM_REAL,DIMENSION(ncells)            :: localG
  40 |         DOM_REAL,DIMENSION(nfacemax,ncells)   :: localH
  41 |         DOM_REAL,DIMENSION(nfacemax,ncells)   :: Le
  42 | 
  43 | !       ----------------------------------------------!
  44 | !       Initialise vectors and rounding error catcher !
  45 | !       ----------------------------------------------!
  46 | 
  47 |         H      = 0.
  48 |         G      = 0.
  49 |         Qr     = 0.
  50 |         nmin   = 0.
  51 | 
  52 |         localH   = 0.
  53 |         localG   = 0.
  54 | 
  55 | !       ----------------------------------!
  56 | !       Looping over all slave directions !
  57 | !       ----------------------------------!
  58 | 
  59 |         DO idir=dir_d, dir_f
  60 | 
  61 | !        IF (pmm_rank.eq.0) print*, "--- dir:",                        &
  62 | !    &                      mu(idir),eta(idir),ksi(idir)
  63 | !        IF (pmm_rank.eq.0) print*, "  pathw:", pathway(:,idir)
  64 | 
  65 |           DO icell=1, ncells
  66 | 
  67 |             j = pathway(icell,idir)
  68 | 
  69 | !           IF (pmm_rank.eq.0) print*, "  ++ Cell",icell,":",j
  70 | 
  71 | !           --------------------------------------------!
  72 | !           Calculate the direction.normal at each face !
  73 | !           --------------------------------------------!
  74 | 
  75 |             DO i=1,nfcelt(j)
  76 |               Dm_ij(i)=norm(1,i,j)*mu(idir)+norm(2,i,j)*eta(idir)+      &
  77 |      &                 norm(3,i,j)*ksi(idir)
  78 |             ENDDO
  79 | 
  80 | !           ----------------------------------------------!
  81 | !           Calculate intensity at exit faces of the cell !
  82 | !           ----------------------------------------------!
  83 |                 
  84 |             CALL MFSCHEME(k_abs(j),k_scat(j),V(j),nfcelt(j),S(:,j),     &
  85 |      &                    Dm_ij,Lb(j),Gi(j),Li(:,j),Le(:,j),Lpi,alpha,  &
  86 |      &                    nfacemax)
  87 | 
  88 |             localG(j) = localG(j)    + Lpi*w(idir)
  89 | 
  90 |             Qr(1,j) = Qr(1,j) + Lpi*w(idir)*mu (idir)
  91 |             Qr(2,j) = Qr(2,j) + Lpi*w(idir)*eta(idir)
  92 |             Qr(3,j) = Qr(3,j) + Lpi*w(idir)*ksi(idir)
  93 | 
  94 | !           IF (pmm_rank.eq.0) print*, "     Ip      =",Lpi
  95 | 
  96 | !           -----------------------------------------------------------!
  97 | !           Update incident intensities at the wall and neighbor cells !
  98 | !           -----------------------------------------------------------!
  99 | 
 100 |             DO i=1,nfcelt(j)                  
 101 | 
 102 |               IF (Dm_ij(i).ge.nmin) THEN
 103 | 
 104 | !               ----------------------------------!
 105 | !               Detect if face is in the boundary !
 106 | !               ----------------------------------!
 107 | 
 108 |                 neighcellid = neighs(2*i-1,j)
 109 | 
 110 | !               ----------------------------------------------!
 111 | !               Apply luminance to the neighbor/boundary face !
 112 | !               ----------------------------------------------!
 113 | 
 114 |                 IF (neighcellid.ne.0) THEN
 115 | 
 116 |                   neighfaceid = neighs(2*i,j)
 117 |                   Li(neighfaceid, neighcellid) = Le(i,j)
 118 | 
 119 |                 ELSE
 120 | 
 121 | !                 IF (pmm_rank.eq.0) print*, " face:", i
 122 |                   localH(i,j)= localH(i,j) + Le(i,j)*w(idir)*Dm_ij(i)
 123 | 
 124 |                 ENDIF
 125 | 
 126 |               ENDIF
 127 | 
 128 |             ENDDO
 129 |           
 130 |           ENDDO
 131 | 
 132 |         ENDDO
 133 | 
 134 | !       -------------------------------------!
 135 | !       Making sommation over ALL directions !
 136 | !       -------------------------------------!
 137 | 
 138 |         CALL MPI_ALLREDUCE(localG, G, ncells, MPI_DOUBLE_PRECISION,     &
 139 |      &                     MPI_SUM, MPI_COMM_WORLD, ierr)
 140 | 
 141 |         CALL MPI_ALLREDUCE(localH, H, ncells*nfacemax,                  &
 142 |      &                     MPI_DOUBLE_PRECISION, MPI_SUM,               &
 143 |      &                     MPI_COMM_WORLD, ierr)
 144 |  
 145 |       END SUBROUTINE SCHEME_DMFS