1 | include(dom.inc)
   2 | 
   3 |        SUBROUTINE SCHEME_DMFS(neighs,nfcelt,wkabs,mu,eta,ksi,w,Dm_ij,   &
   4 |      &                     Lb, Lpi,Gi,G,V,k_scat,k_abs,Qr,Li,Le,H,S,    &
   5 |      &                     epsil,norm,alpha,ndir,ncells,i_gij,nfacemax, &
   6 |      &                     pathway,nkabs)
   7 | 
   8 |        IMPLICIT NONE
   9 | 
  10 |        include 'dom_constants.h'
  11 | 
  12 |        DOM_INT  :: ndir,ncells,i_gij,nfacemax
  13 |        DOM_INT  :: iquad,icell,i,j,n,nkabs
  14 |        DOM_INT, DIMENSION (ncells,2*nfacemax) :: neighs
  15 |        DOM_INT, DIMENSION (ncells) :: nfcelt
  16 | 
  17 |        DOM_INT, DIMENSION(ndir,ncells) :: pathway
  18 |        DOM_REAL, DIMENSION (nkabs) :: wkabs
  19 |        DOM_REAL, DIMENSION(ndir) ::  mu, eta, ksi, w
  20 |        DOM_REAL, DIMENSION(nfacemax) ::  Dm_ij
  21 |        DOM_REAL, DIMENSION(ncells) :: Lb, Lpi, Gi, G, V, k_scat, k_abs
  22 |        DOM_REAL, DIMENSION(ncells,3) :: Qr
  23 |        DOM_REAL, DIMENSION (ncells,nfacemax) :: Li, Le, H, S, epsil
  24 |        DOM_REAL, DIMENSION (ncells,nfacemax,3) ::  norm
  25 |        DOM_REAL  :: alpha
  26 | 
  27 |        DO iquad=1,ndir
  28 | 
  29 | !         print*, " dir ", iquad
  30 | 
  31 |           DO icell=1,ncells
  32 | 
  33 |              j=pathway(iquad,icell)
  34 | 
  35 |              DO i=1,nfcelt(j)
  36 |                 Dm_ij(i)=norm(j,i,1)*mu(iquad)+norm(j,i,2)*             &
  37 |      &                   eta(iquad)+norm(j,i,3)*ksi(iquad)
  38 |              ENDDO
  39 | 
  40 |              CALL MFSCHEME(k_abs(j),wkabs(i_gij),k_scat(j),             &
  41 |      &        V(j),nfcelt(j),S(j,:),Dm_ij,Lb(j),Gi(j),Li(j,:),          &
  42 |      &        epsil(j,:),Le(j,:),Lpi(j),alpha,nfacemax)
  43 | 
  44 |              G(j)=G(j)+Lpi(j)*w(iquad)
  45 |              Qr(j,1)=Qr(j,1)+Lpi(j)*w(iquad)*mu(iquad)
  46 |              Qr(j,2)=Qr(j,2)+Lpi(j)*w(iquad)*eta(iquad)
  47 |              Qr(j,3)=Qr(j,3)+Lpi(j)*w(iquad)*ksi(iquad)
  48 | 
  49 | !           ----------------!
  50 | !           Testing results !
  51 | !           ----------------!
  52 | !
  53 | !           IF (j.eq.5) THEN
  54 | !             print*, " ----"
  55 | !             print*, " w   : ", w(iquad)
  56 | !             print*, " G   : ", G(j)
  57 | !             print*, " Lb  : ", Lb(j)
  58 | !             print*, " Lpi : ", Lpi(j)
  59 | !             print*, " S   : ", S(j,:)
  60 | !             print*, " Dij : ", Dm_ij(:)
  61 | !             print*, " Li  : ", Li(j,:)
  62 | !             print*, " kabs: ", k_abs(j)
  63 | !             print*, " ksca: ", k_scat(j)
  64 | !             print*, " alph: ", alpha
  65 | !             print*, " V   : ", V(j)
  66 | !             print*, " ----"
  67 | !           ENDIF
  68 | 
  69 |              !------------------------------------------------!
  70 |              ! Mise a jour des faces de sorties de la cellule !
  71 |              ! Calcul du flux incident Hw a la parois         !
  72 |              !------------------------------------------------!
  73 | 
  74 |              DO i=1,nfcelt(j)
  75 |                 IF (Dm_ij(i)>=0.) THEN
  76 |                    IF (epsil(j,i)==-1.) THEN
  77 |                       Li(neighs(j,(2*i-1)),neighs(j,(2*i)))=Le(j,i)
  78 | !                     IF (neighs(j,(2*i-1)).eq.5) THEN
  79 | !                       print*, " ******"
  80 | !                       print*, " neigh: ", j
  81 | !                       print*, " Le   : ", Le(j,i)
  82 | !                       print*, " ******"
  83 | !                     ENDIF
  84 |                    ELSE
  85 |                       H(j,i)=H(j,i)+Le(j,i)*w(iquad)*Dm_ij(i)
  86 |                    ENDIF
  87 |                 ENDIF
  88 |              ENDDO
  89 | 
  90 |           ENDDO
  91 |        ENDDO
  92 | 
  93 |        END SUBROUTINE SCHEME_DMFS