1 | include(dom.inc)
   2 | 
   3 |        SUBROUTINE BAND_INTEG(neighs, nfcelt, wkabs, all_k_abs, mu,      &
   4 |      &                 eta, ksi, w, Lb, Gtot, Lbtot, Kp, V, k_scat,     &
   5 |      &                 Srtot, Q_rtot, Lo, Lotot, Htot, S, epsil, norm,  &
   6 |      &                 n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, alpha,&
   7 |      &                 critconv,nfacemax,pathway,maxlen,s_s,ALb,ALo,    &
   8 |      &                 mediumtype,nkabs)
   9 | 
  10 |        IMPLICIT NONE
  11 | 
  12 |        include 'dom_constants.h'
  13 | 
  14 |        DOM_INT :: n_k_abs, ncells, I_SCHEME, ndir,nfacemax
  15 |        DOM_INT :: i_gij, n_iter, i, j,nkabs
  16 |        DOM_INT, DIMENSION (ncells,2*nfacemax) :: neighs
  17 |        DOM_INT, DIMENSION (ncells) :: nfcelt
  18 |        CHARACTER*15 :: mediumtype
  19 | 
  20 |        DOM_INT, DIMENSION(ndir,ncells) :: pathway
  21 |        DOM_REAL, DIMENSION(ndir,ncells) :: maxlen
  22 |        DOM_REAL, DIMENSION(ndir,ncells,nfacemax) :: s_s
  23 |        DOM_REAL,DIMENSION (nkabs) :: wkabs
  24 |        DOM_REAL,DIMENSION(ncells,nkabs) :: all_k_abs
  25 |        DOM_REAL,DIMENSION(ndir) ::  mu, eta, ksi, w
  26 |        DOM_REAL,DIMENSION(nfacemax) ::  Dm_ij
  27 |        DOM_REAL,DIMENSION(ncells) :: Lb, Gtot, V, k_scat, Srtot, Lbtot
  28 |        DOM_REAL,DIMENSION(ncells,3) :: Q_rtot
  29 |        DOM_REAL,DIMENSION(ncells,3) :: Qr
  30 |        DOM_REAL, DIMENSION (ncells,nfacemax) :: Lo, Lotot, Htot,S, epsil
  31 |        DOM_REAL, DIMENSION (ncells,nfacemax,3) ::  norm
  32 |        DOM_REAL :: DWVNB_SI, alpha
  33 |        DOM_REAL,DIMENSION(ncells) :: Lpi, Gi, Gj, G, k_abs, Sr, Kp
  34 |        DOM_REAL, DIMENSION (ncells,nfacemax) :: Li, Le, L, H
  35 |        DOM_REAL :: error,critconv
  36 | 
  37 |        DOM_REAL,DIMENSION(ncells,ngg) :: ALb
  38 |        DOM_REAL,DIMENSION(ncells,nfacemax,ngg) :: ALo
  39 | 
  40 |        print*, " Doing narrow band integration"
  41 | 
  42 |        Kp = 0.
  43 | 
  44 |        DO i_gij=1,n_k_abs
  45 | 
  46 |          G=0.
  47 |          Sr=0.
  48 |          Lpi=0.
  49 |          H=0.
  50 |          Le=0.
  51 | 
  52 |          IF (mediumtype.eq.'WSGG') THEN
  53 |            !********************************************************
  54 |            ! Pour le traitement WSGG
  55 |            ! Le traitement de chaque gaz se fait sur une fraction du
  56 |            ! corps noir A*Lb (avec sum(A)=1)
  57 |            !********************************************************
  58 |            Lb(:)=ALb(:,i_gij)
  59 |            Lo(:,:)=epsil(:,:)*ALo(:,:,i_gij)
  60 |          ENDIF
  61 | 
  62 | 
  63 |          !-----------------------------------------------------------!
  64 |          ! POINT DE MISE A JOUR PAR ITERATION (REFLECTION/DIFFUSION) !
  65 |          !-----------------------------------------------------------!
  66 |          error=1.d0
  67 |          n_iter=0
  68 | 
  69 |          DO WHILE(error.gt.critconv)
  70 | 
  71 | !            -------------------------------------------!
  72 | !            Update values for reflection sub iteration !
  73 | !            -------------------------------------------!
  74 |              Gi=G
  75 | 
  76 | !*******************************************************************!
  77 | ! Introduction de la diffusion au niveau du Terme Source S          !
  78 | ! Stockage des luminances au faces mailles et de la luminance       !
  79 | !       calculees pour chaque maille a l'iteration precedente dans  !
  80 | !       les Ij et Gi                                                !
  81 | ! Mise a jour de la luminance aux faces des mailles: Ii=I calculees !
  82 | !       aux parois                                                  !
  83 | ! On remet I=Io et Q=Io pour calculer resp. la luminance emise aux  !
  84 | !       parois avec reflexion et le flux aux parois                 !
  85 | ! Mise a zero des matrices G et Qaxes                               !
  86 | !*******************************************************************!
  87 | 
  88 |              DO i=1,ncells
  89 |                 DO j=1,nfcelt(i)
  90 |                    IF (epsil(i,j)/=-1.) THEN
  91 |                       Li(i,j)=Lo(i,j)+(1-epsil(i,j))/pi*H(i,j)
  92 |                    ELSE
  93 |                       Li(i,j)=Le(i,j)
  94 |                    ENDIF
  95 |                 ENDDO
  96 | !               ----------------!
  97 | !               Testing results !
  98 | !               ----------------!
  99 | !               IF (i.eq.5) THEN
 100 | !                 print*, " Lo = ", Lo(i,:)
 101 | !                 print*, " Li = ", Li(i,:)
 102 | !                 print*, " ===="
 103 | !               ENDIF
 104 | 
 105 |              ENDDO
 106 | 
 107 |              G=0.
 108 |              H=0.
 109 |              Qr=0.
 110 | 
 111 |              !----------------------------------------!
 112 |              ! Lecture des itineraires de resolutions !
 113 |              !----------------------------------------!
 114 | 
 115 |              k_abs=all_k_abs(:,i_gij)
 116 | 
 117 |              SELECT CASE (I_SCHEME)
 118 | 
 119 |                !------------------!
 120 |                ! MEAN FLUX SCHEME !
 121 |                !------------------!
 122 |                CASE(1)
 123 | 
 124 |                  CALL SCHEME_DMFS(neighs,nfcelt,wkabs,mu,eta,ksi,w,     &
 125 |      &               Dm_ij,Lb, Lpi,Gi,G,V,k_scat,k_abs,Qr,              &
 126 |      &              Li,Le,H,S,epsil,norm,alpha,ndir,ncells,i_gij,       &
 127 |      &              nfacemax,pathway,nkabs)
 128 | 
 129 |                !--------------------!
 130 |                ! EXPONENTIAL SCHEME !
 131 |                !--------------------!
 132 |                CASE(2)
 133 | 
 134 |                CALL SCHEME_EXP(neighs,nfcelt, wkabs,mu, eta, ksi, w,    &
 135 |      &                    Dm_ij,Lb, Lpi, Gi, G, V, k_scat, k_abs,       &
 136 |      &                    Li, Le, H,S, epsil, norm, ndir, ncells,i_gij, &
 137 |      &                    nfacemax,pathway,maxlen,s_s,nkabs)
 138 | 
 139 |              END SELECT
 140 | 
 141 | !********************************************************************
 142 | !                   Module d'iteration avec reflexion               *
 143 | !********************************************************************
 144 | !                                                                   *
 145 | ! Tant qu'il existera un element parietal dont la valeur sera       *
 146 | ! modifiee a plus du pourcentage fixe par critconv la convergence   *
 147 | ! ne sera pas satisfaite et la mise a jour des radiosites sera      *
 148 | ! conservee pour une nouvelle estimation des champs de luminance    *
 149 | !                                                                   *
 150 | !********************************************************************
 151 | 
 152 |              DO i=1,ncells
 153 | 
 154 |                 IF ((n_iter/=0).and.(G(i)/=0.)) THEN
 155 |                    Gj(i)=abs((G(i)-Gi(i))/G(i))
 156 |                 ELSE
 157 |                    Gj(i)=1.
 158 |                 ENDIF
 159 | 
 160 |              ENDDO
 161 | 
 162 |              error=maxval(Gj)
 163 |              n_iter=n_iter+1    ! Compteur d'iteration mis a jour
 164 |          ENDDO
 165 | 
 166 |          Gtot=Gtot+G*wkabs(i_gij)*DWVNB_SI !(i_bande)
 167 |          Htot=Htot+H*wkabs(i_gij)*DWVNB_SI
 168 |          Lotot=Lotot+Lo*wkabs(i_gij)*DWVNB_SI
 169 |          Sr=(4*pi*Lb-G)*all_k_abs(:,i_gij)*wkabs(i_gij)*DWVNB_SI
 170 |          Lbtot=Lbtot+4*pi*Lb*all_k_abs(:,i_gij)*wkabs(i_gij)*DWVNB_SI
 171 |          Srtot=Srtot+Sr
 172 | 
 173 |          Q_rtot(:,1) = Q_rtot(:,1)+Qr(:,1)*k_abs*wkabs(i_gij)*DWVNB_SI
 174 |          Q_rtot(:,2) = Q_rtot(:,2)+Qr(:,2)*k_abs*wkabs(i_gij)*DWVNB_SI
 175 |          Q_rtot(:,3) = Q_rtot(:,3)+Qr(:,3)*k_abs*wkabs(i_gij)*DWVNB_SI
 176 | 
 177 |          Kp = Kp + k_abs*wkabs(i_gij)
 178 | 
 179 | !        print*, " ++++++"
 180 | !        print*, " Gmax = ", MAXVAL(Gtot)
 181 | !        print*, " Gmin = ", MINVAL(Gtot)
 182 | !        print*, " Smax = ", MAXVAL(Srtot)
 183 | !        print*, " Smin = ", MINVAL(Srtot)
 184 | !        print*, " ++++++"
 185 | 
 186 | 
 187 |        ENDDO
 188 | 
 189 |        END SUBROUTINE BAND_INTEG