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