1 | include(dom.inc)
2 |
3 | SUBROUTINE BAND_INTEG(neighs, nfcelt, wkabs, all_k_abs, mu, eta, &
4 | & ksi, w, Lb, WSGG_W, V, k_scat, maxlen, Lo, S, &
5 | & epsil, s_s, norm, n_k_abs, ncells, nfacemax, &
6 | & I_SCHEME, ndir, DWVNB_SI, alpha, nkabs, &
7 | & pathway, Gtot, Htot, Lotot, Lbtot, Q_rtot, dir_d, &
8 | & dir_f, ngg, mediumtype, WVNB_SI, critconv, bcell, &
9 | & bface, nbface)
10 |
11 | ! ONLY for print test
12 | USE mod_pmm
13 |
14 | IMPLICIT NONE
15 |
16 | include 'dom_constants.h'
17 | include 'pmm_constants.h'
18 |
19 | ! IN
20 | DOM_INT :: nkabs, n_k_abs, ncells, nfacemax, ngg
21 | DOM_INT :: I_SCHEME, ndir, dir_d, dir_f, nbface
22 | DOM_INT, DIMENSION (2*nfacemax,ncells):: neighs
23 | DOM_INT, DIMENSION (ncells) :: nfcelt
24 | DOM_INT, DIMENSION (ncells,ndir) :: pathway
25 |
26 | DOM_REAL :: DWVNB_SI, alpha, critconv
27 | DOM_REAL :: Gmax, Gmin, Smax, Smin
28 | DOM_REAL :: WVNB_SI
29 | DOM_REAL,DIMENSION(nkabs) :: wkabs
30 | DOM_REAL,DIMENSION(nkabs,ncells) :: all_k_abs
31 | DOM_REAL,DIMENSION(ndir) :: mu, eta, ksi, w
32 | DOM_REAL,DIMENSION(ncells) :: Lb, V, k_scat
33 | DOM_REAL,DIMENSION(ncells,ndir) :: maxlen
34 | DOM_REAL,DIMENSION(nfacemax,ncells) :: Lo
35 | DOM_REAL,DIMENSION(nfacemax,ncells) :: S
36 | DOM_REAL,DIMENSION(nbface) :: epsil
37 | DOM_INT ,DIMENSION(nbface) :: bcell, bface
38 | DOM_REAL,DIMENSION(3,nfacemax,ncells) :: norm
39 | DOM_REAL,DIMENSION(nfacemax,ncells,ndir) ::s_s
40 | DOM_REAL,DIMENSION(ngg,ncells) :: WSGG_W
41 |
42 | CHARACTER*80 :: mediumtype
43 |
44 | ! OUT
45 | DOM_REAL,DIMENSION(3,ncells) :: Q_rtot
46 | DOM_REAL,DIMENSION(ncells) :: Gtot, Lbtot
47 | DOM_REAL,DIMENSION(nfacemax,ncells) :: Lotot, Htot
48 |
49 | ! LOCAL
50 | DOM_INT :: i_gij, n_iter, i, j, maxiter, ibnd
51 | DOM_INT :: mloc(1)
52 |
53 | DOM_REAL :: aiminan, planck, blae, error
54 | DOM_REAL,DIMENSION(ncells) :: Gi,Gj,G
55 | DOM_REAL,DIMENSION(ncells) :: k_abs
56 | DOM_REAL,DIMENSION(nfacemax,ncells) :: Li, H
57 | DOM_REAL,DIMENSION(3,ncells) :: Qr
58 | DOM_REAL,DIMENSION(nfacemax,ncells) :: Lo_org
59 | DOM_REAL,DIMENSION(ncells) :: Lb_org
60 |
61 | ! print*, " (",pmm_rank,") Doing band integration"
62 |
63 | ! --------------------------!
64 | ! Integrate the narrow band !
65 | ! --------------------------!
66 |
67 | Lo_org = Lo
68 | Lb_org = Lb
69 |
70 | DO i_gij=1,n_k_abs
71 |
72 | ! -----------------------------------------------------!
73 | ! Initialisation values before reflection subiteration !
74 | ! -----------------------------------------------------!
75 |
76 | G = 0.
77 | H = 0.
78 | Li = 0.
79 |
80 | error = 1.
81 | n_iter = 0
82 | maxiter = 500
83 |
84 | DO WHILE ((error.gt.critconv).and.(n_iter.lt.maxiter))
85 |
86 | ! -------------------------------------------!
87 | ! Update values for reflection sub iteration !
88 | ! -------------------------------------------!
89 |
90 | Gi = G
91 |
92 | ! --------------------------------------------------!
93 | ! WSGG: Le traitement de chaque gaz se fait sur une !
94 | ! fraction du corps noir A*Lb (avec sum(A)=1) !
95 | ! --------------------------------------------------!
96 | IF (trim(mediumtype).eq.'WSGG') THEN
97 | DO j = 1, ncells
98 | Lb(j) = WSGG_W(i_gij,j) * Lb_org(j)
99 | DO i = 1, nfcelt(j)
100 | Lo(i,j) = WSGG_W(i_gij,j) * Lo_org(i,j)
101 | ENDDO
102 | ENDDO
103 | ENDIF
104 |
105 | ! --------------------------------------------!
106 | ! Update faces "in" lumminance for reflection !
107 | ! --------------------------------------------!
108 |
109 | ibnd = 1
110 |
111 | DO WHILE (ibnd.le.nbface)
112 |
113 | j = bcell(ibnd)
114 | i = bface(ibnd)
115 |
116 | Li(i,j) = Lo(i,j) + (1.-epsil(ibnd))*H(i,j)/pi
117 | ibnd = ibnd + 1
118 |
119 | ENDDO
120 |
121 | ! ---------------------------------------!
122 | ! Lecture des itineraires de resolutions !
123 | ! ---------------------------------------!
124 |
125 | k_abs=all_k_abs(i_gij,:)
126 |
127 | SELECT CASE (I_SCHEME)
128 |
129 | ! -----------------!
130 | ! MEAN FLUX SCHEME !
131 | ! -----------------!
132 | CASE(1)
133 |
134 | ! print*," (",pmm_rank,") Doing DMFS spatial discretization"
135 | CALL SCHEME_DMFS(neighs,nfcelt,mu,eta,ksi,w,pathway,Lb, &
136 | & Gi,G,V,k_scat,k_abs,Qr,Li,H,S,norm, &
137 | & alpha,ndir,ncells,nfacemax,nkabs,dir_d,dir_f, &
138 | & bcell, bface, nbface)
139 |
140 | ! -------------------!
141 | ! EXPONENTIAL SCHEME !
142 | ! -------------------!
143 | CASE(2)
144 |
145 | CALL SCHEME_EXPO(neighs, nfcelt, mu, eta, ksi, w, pathway,&
146 | & Lb, Gi, G, V, k_scat, k_abs, maxlen, Qr, Li, H, S, &
147 | & s_s,norm, ndir, nfacemax,ncells, dir_d, dir_f, &
148 | & bcell, bface, nbface)
149 |
150 | END SELECT
151 |
152 | ! ---------------------------------------!
153 | ! Convergence test for reflexion problem !
154 | ! ---------------------------------------!
155 |
156 | mloc = MAXLOC(abs(G-Gi))
157 | IF ((n_iter.ne.0).and.((abs(Gi(mloc(1)))).ne.0.)) THEN
158 | error = (abs(G(mloc(1))-Gi(mloc(1))))/abs(Gi(mloc(1)))
159 | ELSE
160 | error = 1.
161 | ENDIF
162 |
163 | n_iter = n_iter + 1
164 |
165 | ! IF (pmm_rank.eq.0) print*, "Nq=",i_gij,"/",n_k_abs,"-", &
166 | ! & "n_iter=",n_iter,"- err=", error
167 |
168 | ! IF (pmm_rank.eq.0) print*, "mloc=",mloc(1),"- G=", &
169 | ! & G(mloc(1)),"- Gi=",Gi(mloc(1))
170 |
171 | ENDDO
172 |
173 | ! ----------------------------------------!
174 | ! Sum over all slave directions and bands !
175 | ! ----------------------------------------!
176 |
177 | G = G / pmm_n_p
178 | H = H / pmm_n_p
179 |
180 | ! Qr = Qr / pmm_n_p
181 |
182 | Gtot = Gtot + G*k_abs*wkabs(i_gij)*DWVNB_SI
183 | Htot = Htot + H*wkabs(i_gij)*DWVNB_SI
184 | Lotot = Lotot + Lo*wkabs(i_gij)*DWVNB_SI
185 | Lbtot = Lbtot + Lb*4*pi*k_abs*wkabs(i_gij)*DWVNB_SI
186 |
187 | Q_rtot(1,:) = Q_rtot(1,:)+Qr(1,:)*wkabs(i_gij)*DWVNB_SI
188 | Q_rtot(2,:) = Q_rtot(2,:)+Qr(2,:)*wkabs(i_gij)*DWVNB_SI
189 | Q_rtot(3,:) = Q_rtot(3,:)+Qr(3,:)*wkabs(i_gij)*DWVNB_SI
190 |
191 | ENDDO
192 |
193 | END SUBROUTINE BAND_INTEG
band_integ.F could be called by: