1 | include(dom.inc)
2 |
3 | program PRISSMA
4 |
5 | !-----------!
6 | ! Variables !
7 | !----------------------------------------------------------------------!
8 | ! Lpi = Luminance emise par la maille dans la direction iquad
9 | ! Lb = Luminance noire du gaz emise au centre de chaque maille
10 | ! Lo = Luminances init aux faces des mailles pour une bande etroite
11 | ! Lotot = Luminances initiales totale aux faces des mailles
12 | ! Li = Luminances aux faces des mailles : en court
13 | ! L_ = Luminances aux faces d mail:futures aux parois (si itration)
14 | ! H = Eclairement aux parois (W.m-2 pour une bande troite
15 | ! Htot = Eclairement totale aux parois (W.m-2)
16 | ! Spi = Terme source (emission et diffusion) au centre de maille ds
17 | ! la direction iquad (W.m-2.str-1)
18 | ! G = Raynnmt incident centre de maille actuel pr une bnd etroite
19 | ! Gtot = Rayonnement incident totale au centre de maille actuel
20 | ! Gi = Rayonnement incident (W.m-2): a l'iteration precedente
21 | ! Gj = Raynnmnt incid(W.m-2): correction apportee par l'iteration
22 | ! Qin = Composantes vecteur Qr sur chaque axe
23 | ! Qintot = Composantes vecteur Qr sur chaque axe
24 | ! Sr = Terme source de chaleur par rayonnement (W.m-3.str-1)
25 | ! Srtot = Composantes vecteur Qr sur chaque axe
26 | ! mu, eta, ksi, w =coordonnees sr X, Y et Z e poids associe du vecteur
27 | ! unitaire lie a iquad
28 | ! epsil = emissivite aux parois
29 | ! a , V = aire des faces et volume d'une maille
30 | ! norm = vecteurs normaux aux faces de chaques mailles
31 | ! celldata = Donne caract. d'une maille physq : T,P,Xh20,Xco2,Xco,Xsuies
32 | ! KCO,DCO = Donnees SNB de l'EM2C pour le CO
33 | ! KC,DC = Donnees SNB de l'EM2C pour le CO2
34 | ! KH,DH = Donnees SNB de l'EM2C pour le H2O
35 | ! k_abs = Absorpt coeffs discretized over the narrow band rearranged
36 | ! wkabs =
37 | ! ka_g_i =
38 | !-----------------------------------------------------------------------!
39 |
40 | IMPLICIT NONE
41 |
42 | include 'dom_constants.h'
43 |
44 | CHARACTER*15 :: mediumtype, datastatus
45 | CHARACTER*3 :: homosyst
46 | CHARACTER*80 :: path, gfile
47 | CHARACTER*80 :: pathspec
48 | CHARACTER*80 :: outpath
49 | LOGICAL :: LICO, LICO2, LIH2O
50 | DOM_INT :: nkabs
51 | DOM_INT :: finprog, ncells, i_scheme, ndir, nbounds, n_k_abs
52 | DOM_INT :: i_bande, ICO2, ICO, IH2O, end_prog,nfacemax,ngbands
53 | DOM_INT :: i_nnodes,i_nfaces,i_nfacesmax
54 | DOM_REAL :: kbar, phi_f, DWVNB
55 | DOM_REAL :: critconv, alpha, DWVNB_SI, WVNB, WVNB_SI
56 | DOM_REAL,DIMENSION(2) :: SNBDATA
57 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: wkabs, ka_g_i, k_absmel
58 | DOM_REAL,DIMENSION(n_SNBmax) :: WV
59 | DOM_REAL,DIMENSION(n_SNBmax) :: all_WVNB,all_DWVNB
60 | DOM_REAL,DIMENSION(14,n_SNBmax) :: KCO,DCO,KC,DC,KH,DH
61 | DOM_REAL,DIMENSION(ngg) :: k_wsgg
62 | DOM_REAL,DIMENSION(ngg,6) :: alpha_wsgg
63 |
64 | DOM_INT, ALLOCATABLE,DIMENSION(:,:) :: neighs
65 | DOM_INT, ALLOCATABLE,DIMENSION(:) :: nfcelt,npcelt
66 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: WVSOOT,xc,yc,zc,V
67 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: Lb,Lpi,Spi,Gi,Gj,G,Gtot
68 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: Kp, Kptot
69 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: Lbtot,k_scat,k_abs,Sr
70 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: kabs_gray
71 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: T,Srtot,Qin,Qintot,wk
72 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: Gtotpar,Srtotpar
73 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:) :: celldata,all_k_abs
74 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:) :: Q_rtot,H,Htot,S
75 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:) :: Lo,Lotot,Li,Le,L
76 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:) :: epsil,Tf,Htotpar
77 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:) :: Lototpar,ALb
78 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:) :: WSGG_W,xfc,yfc,zfc
79 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:,:) :: norm,ALo
80 | DOM_REAL, ALLOCATABLE,DIMENSION(:) :: mu, eta, ksi, w
81 |
82 | DOM_INT , DIMENSION(:,:), ALLOCATABLE :: pathway
83 | DOM_REAL, DIMENSION(:,:), ALLOCATABLE :: maxlen
84 | DOM_REAL, DIMENSION(:,:,:),ALLOCATABLE :: s_s
85 |
86 |
87 | ! -------------------------------------!
88 | ! Get input files path and global data !
89 | ! -------------------------------------!
90 |
91 | OPEN(UNIT=1,FILE='prissma.choices')
92 | READ(1,*) path
93 | READ(1,*)
94 | READ(1,*)
95 | READ(1,*)
96 | READ(1,*)
97 | READ(1,*)
98 | READ(1,*) nkabs
99 | CLOSE(1)
100 |
101 | gfile = path(1:len_trim(path))//'/Global.in'
102 |
103 | OPEN(UNIT=10,FILE=gfile,FORM='unformatted')
104 | READ(10) ncells,i_nnodes,i_nfaces,nfacemax,ndir
105 | CLOSE(10)
106 |
107 | ! ------------------!
108 | ! Allocate vecrtors !
109 | ! ------------------!
110 |
111 | IF (ALLOCATED(pathway)) DEALLOCATE(pathway)
112 | IF (ALLOCATED(maxlen)) DEALLOCATE(maxlen)
113 | IF (ALLOCATED(s_s)) DEALLOCATE(s_s)
114 |
115 | ALLOCATE(pathway(ndir,ncells))
116 | ALLOCATE(maxlen(ndir,ncells))
117 | ALLOCATE(s_s(ndir,ncells,nfacemax))
118 |
119 | ALLOCATE(neighs(ncells,2*nfacemax))
120 | ALLOCATE(nfcelt(ncells))
121 | ALLOCATE(npcelt(ncells))
122 | ALLOCATE(WVSOOT(ncells))
123 | ALLOCATE(Lb(ncells))
124 | ALLOCATE(Lpi(ncells))
125 | ALLOCATE(Spi(ncells))
126 | ALLOCATE(Gi(ncells))
127 | ALLOCATE(Gj(ncells))
128 | ALLOCATE(G(ncells))
129 | ALLOCATE(Gtot(ncells))
130 | ALLOCATE(Kp(ncells))
131 | ALLOCATE(Kptot(ncells))
132 | ALLOCATE(T(ncells))
133 | ALLOCATE(Lbtot(ncells))
134 | ALLOCATE(xc(ncells))
135 | ALLOCATE(yc(ncells))
136 | ALLOCATE(zc(ncells))
137 | ALLOCATE(V(ncells))
138 | ALLOCATE(k_scat(ncells))
139 | ALLOCATE(k_abs(ncells))
140 | ALLOCATE(Sr(ncells))
141 | ALLOCATE(Srtot(ncells))
142 | ALLOCATE(Qin(ncells))
143 | ALLOCATE(Qintot(ncells))
144 | ALLOCATE(wk(ncells))
145 | ALLOCATE(Gtotpar(ncells))
146 | ALLOCATE(Srtotpar(ncells))
147 | ALLOCATE(kabs_gray(ncells))
148 | ALLOCATE(Q_rtot(ncells,3))
149 | ALLOCATE(celldata(ncells,8))
150 | ALLOCATE(all_k_abs(ncells,nkabs))
151 | ALLOCATE(Lo(ncells,nfacemax))
152 | ALLOCATE(Lotot(ncells,nfacemax))
153 | ALLOCATE(Li(ncells,nfacemax))
154 | ALLOCATE(Le(ncells,nfacemax))
155 | ALLOCATE(L(ncells,nfacemax))
156 | ALLOCATE(H(ncells,nfacemax))
157 | ALLOCATE(Htot(ncells,nfacemax))
158 | ALLOCATE(S(ncells,nfacemax))
159 | ALLOCATE(epsil(ncells,nfacemax))
160 | ALLOCATE(Tf(ncells,nfacemax))
161 | ALLOCATE(xfc(ncells,nfacemax))
162 | ALLOCATE(yfc(ncells,nfacemax))
163 | ALLOCATE(zfc(ncells,nfacemax))
164 | ALLOCATE(Htotpar(ncells,nfacemax))
165 | ALLOCATE(Lototpar(ncells,nfacemax))
166 | ALLOCATE(norm(ncells,nfacemax,3))
167 | ALLOCATE(ALb(ncells,ngg))
168 | ALLOCATE(ALo(ncells,nfacemax,ngg))
169 | ALLOCATE(WSGG_W(ncells,ngg))
170 | ALLOCATE(mu(ndir))
171 | ALLOCATE(eta(ndir))
172 | ALLOCATE(ksi(ndir))
173 | ALLOCATE(w(ndir))
174 | ALLOCATE(wkabs(nkabs))
175 | ALLOCATE(ka_g_i(nkabs))
176 | ALLOCATE(k_absmel(nkabs))
177 |
178 | ! --------------------------------------------------!
179 | ! Read data from input files and from the data base !
180 | ! --------------------------------------------------!
181 |
182 | CALL READ_DATA(mediumtype,datastatus,neighs,nfcelt,KCO,DCO,KC, &
183 | & DC,KH,DH,xc, yc, zc, V, k_scat, kabs_gray,S, &
184 | & epsil, xfc, yfc, zfc, Tf,norm,celldata,mu, eta, &
185 | & ksi, w,critconv,alpha,ncells,I_SCHEME,ndir, &
186 | & nbounds,end_prog, all_WVNB,all_DWVNB,alpha_wsgg,&
187 | & k_wsgg,homosyst,nfacemax,pathway,maxlen,s_s, &
188 | & path,pathspec,outpath,n_SNBmax)
189 |
190 | finprog=0
191 |
192 | ! ------------------------------!
193 | ! Results vectos initialisation !
194 | ! ------------------------------!
195 |
196 | Kptot=0.
197 | Gtot=0.
198 | Lbtot=0.
199 | Htot=0.
200 | Lotot=0.
201 | Srtot=0.
202 | Q_rtot=0.
203 |
204 | ! ---------------------!
205 | ! GRAY CASE TREATEMENT !
206 | ! ---------------------!
207 |
208 | IF (mediumtype.eq.'GRAY') THEN
209 | n_k_abs=1
210 | wkabs=1.0
211 | DWVNB_SI=1.0
212 |
213 | CALL EMISSIV(nfcelt, celldata, Lb, Lo, epsil, Tf, ncells, &
214 | & nfacemax)
215 |
216 | CALL GRAY_CASE(all_k_abs,kabs_gray,ncells,nkabs,homosyst)
217 |
218 | CALL BAND_INTEG(neighs, nfcelt, wkabs, all_k_abs, mu, eta, &
219 | & ksi,w,Lb,Gtot,Lbtot,Kp,V,k_scat,Srtot, &
220 | & Q_rtot, Lo, Lotot, Htot, S, epsil, norm, &
221 | & n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, &
222 | & alpha,critconv,nfacemax,pathway,maxlen,s_s, &
223 | & ALb,ALo,mediumtype,nkabs)
224 |
225 | Kptot = Kp
226 |
227 | ! ---------------------!
228 | ! WSGG CASE TREATEMENT !
229 | ! ---------------------!
230 |
231 | ELSEIF (mediumtype.eq.'WSGG') THEN
232 | n_k_abs=ngg
233 | wkabs=1.0
234 | DWVNB_SI=1.0
235 |
236 | CALL EMISSIV(nfcelt, celldata, Lb, Lo, epsil, Tf, ncells, &
237 | & nfacemax)
238 |
239 | CALL WSGG_CASE(all_k_abs,wkabs,Lb,celldata,ncells,nfcelt,Tf, &
240 | & ALb,ALo,WSGG_W,alpha_wsgg,k_wsgg,all_WVNB, &
241 | & all_DWVNB,nfacemax,nkabs)
242 |
243 | CALL BAND_INTEG(neighs, nfcelt, wkabs, all_k_abs, mu, eta, &
244 | & ksi,w,Lb,Gtot,Lbtot,Kp,V,k_scat,Srtot, &
245 | & Q_rtot, Lo, Lotot, Htot, S, epsil, norm, &
246 | & n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, &
247 | & alpha,critconv,nfacemax,pathway,maxlen,s_s,
248 | & ALb,ALo,mediumtype,nkabs)
249 |
250 | Kptot = Kp
251 |
252 | ! ------------------------!
253 | ! SNB-FSK CASE TREATEMENT !
254 | ! ------------------------!
255 |
256 | ELSEIF (mediumtype.eq.'FSK') THEN
257 | n_k_abs=nkabs
258 | DWVNB_SI=1.0
259 | ngbands=n_SNBmax
260 |
261 | CALL EMISSIV(nfcelt, celldata, Lb, Lo, epsil, Tf, ncells, &
262 | & nfacemax)
263 |
264 | CALL FSK_CASE(all_k_abs,wkabs,Lb,celldata,ncells,all_WVNB, &
265 | & all_DWVNB,KCO,DCO,KC,DC,KH,DH,homosyst, &
266 | & ngbands,n_SNBmax,nkabs)
267 |
268 | CALL BAND_INTEG(neighs, nfcelt, wkabs, all_k_abs, mu, eta, &
269 | & ksi,w,Lb, Gtot,Lbtot,Kp,V,k_scat,Srtot, &
270 | & Q_rtot, Lo, Lotot, Htot, S, epsil, norm, &
271 | & n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, &
272 | & alpha,critconv,nfacemax,pathway,maxlen,s_s, &
273 | & ALb,ALo,mediumtype,nkabs)
274 |
275 | Kptot = Kp
276 |
277 | ! -------------------!
278 | ! CK CASE TREATEMENT !
279 | ! -------------------!
280 |
281 | ELSEIF (mediumtype.eq.'CK') THEN
282 |
283 | ! ------------------------------!
284 | ! Process each one of the bands !
285 | ! ------------------------------!
286 |
287 | DO i_bande=1,n_SNBmax
288 |
289 | print*, 'Bande: ',i_bande,'/',n_SNBmax
290 |
291 | n_k_abs = nkabs
292 | ngbands = 1
293 | WVNB = all_WVNB(i_bande)
294 | DWVNB = all_DWVNB(i_bande)
295 | WVNB_SI = 100.*WVNB
296 | DWVNB_SI = 100.*DWVNB
297 |
298 | CALL EMISSIV_SNB(nfcelt,celldata,Lb,Lo,epsil,Tf,WVNB_SI, &
299 | & ncells,nfacemax)
300 |
301 | CALL CK_CASE(all_k_abs,wkabs,Lb,celldata,ncells,all_WVNB, &
302 | & all_DWVNB,KCO,DCO,KC,DC,KH,DH,homosyst, &
303 | & ngbands,i_bande,n_SNBmax,nkabs)
304 |
305 | ! ----------------------------------------------!
306 | ! Band integration using nkabs point quadrature !
307 | ! ----------------------------------------------!
308 |
309 | print*, " Doing band integration..."
310 | print*, mediumtype
311 | CALL BAND_INTEG(neighs,nfcelt,wkabs,all_k_abs,mu,eta,ksi, &
312 | & w,Lb,Gtot,Lbtot,Kp,V,k_scat,Srtot, &
313 | & Q_rtot, Lo, Lotot, Htot, S, epsil, norm, &
314 | & n_k_abs, ncells, I_SCHEME,ndir,DWVNB_SI, &
315 | & alpha,critconv,nfacemax,pathway,maxlen, &
316 | & s_s,ALb,ALo,mediumtype,nkabs)
317 |
318 | Kptot = Kptot + Kp*(Lb*pi*DWVNB_SI) &
319 | & *pi/(sigma*celldata(:,1)**4.d0)
320 |
321 | ENDDO
322 |
323 | ELSE
324 | finprog=1
325 | PRINT*, 'Error no computation done'
326 | PRINT*, 'fill the choice file correctly and run again'
327 | ENDIF
328 |
329 | ! ----------------------------!
330 | ! Out-processing and finalize !
331 | ! ----------------------------!
332 |
333 | IF (finprog .ne. 1) THEN
334 |
335 | CALL OUTPROCESSING(ncells,nbounds,xc,yc,zc,Gtot,Srtot,Htot, &
336 | & Lotot,Q_rtot,Lb,Kptot,epsil,nfacemax,path,&
337 | & outpath)
338 |
339 | ENDIF
340 |
341 | END program PRISSMA
prissma.F could be called by: