1 | include(dom.inc)
2 |
3 | SUBROUTINE BAND_INTEG(wkabs, kabs, kscat, Lb, Lo, DWVNB_SI, &
4 | & Gtot, Lbtot, Htot, Q_rtot, Q_ptot)
5 |
6 | USE mod_slave
7 | USE mod_inout
8 | USE mod_pmm
9 |
10 | IMPLICIT NONE
11 |
12 | include 'dom_constants.h'
13 | include 'pmm_constants.h'
14 |
15 | ! IN
16 | DOM_REAL :: DWVNB_SI
17 | DOM_REAL,DIMENSION(is_nkabs) :: wkabs
18 | DOM_REAL,DIMENSION(is_nkabs,is_ncells):: kabs
19 | DOM_REAL,DIMENSION(is_ncells) :: Lb, kscat
20 | DOM_REAL,DIMENSION(is_nbfaces) :: Lo
21 |
22 | ! OUT
23 | DOM_REAL,DIMENSION(3,is_ncells) :: Q_rtot
24 | DOM_REAL,DIMENSION(is_ncells) :: Gtot, Lbtot
25 | DOM_REAL,DIMENSION(is_nbfaces) :: Htot
26 | DOM_REAL,DIMENSION(3,is_nprobes) :: Q_ptot
27 |
28 | ! LOCAL
29 | DOM_INT :: iquad, i
30 | DOM_INT :: n_iter, maxiter
31 | DOM_REAL :: error,error1,error2
32 | DOM_REAL,ALLOCATABLE, DIMENSION(:,:) :: WSGG_W_cell
33 | DOM_REAL,ALLOCATABLE, DIMENSION(:,:) :: WSGG_W_face
34 | DOM_REAL,DIMENSION(is_nkabs,is_ncells) :: Gi
35 | DOM_REAL,DIMENSION(is_ncells) :: G
36 | DOM_REAL,DIMENSION(is_nbfaces) :: H, Hitot
37 | DOM_REAL,DIMENSION(3,is_ncells) :: Qr
38 | DOM_REAL,DIMENSION(3,is_nprobes) :: Qp
39 | DOM_REAL,ALLOCATABLE, DIMENSION(:,:,:) :: Livirt
40 | ! DOM_INT :: mloc(1) , mloc2(3)
41 |
42 | IF (.not.ALLOCATED(WSGG_W_cell)) &
43 | & ALLOCATE(WSGG_W_cell(is_nkabs, is_ncells))
44 | IF (.not.ALLOCATED(WSGG_W_face)) &
45 | & ALLOCATE(WSGG_W_face(is_nkabs, is_nbfaces))
46 |
47 | CALL GATHER(s_WSGG_W, WSGG_W_cell, is_nkabs)
48 | CALL GATHER_FACES(s_WSGG_W, WSGG_W_face, is_nkabs)
49 |
50 | IF((i_dom_npart.gt.1).and.(.not.ALLOCATED(Livirt))) THEN
51 | ALLOCATE(Livirt(is_ntot_vfaces,is_ntotdir, is_nkabs))
52 | ENDIF
53 |
54 | ! -----------------------------------------------------!
55 | ! Initialisation values before reflection subiteration !
56 | ! -----------------------------------------------------!
57 |
58 | error = 1000.
59 | n_iter = 0
60 | maxiter = 50
61 | H = 0.
62 |
63 | DO WHILE ((error.gt.critconv).and.(n_iter.lt.maxiter))
64 |
65 | ! -------------------------------------------!
66 | ! Update values for reflection sub iteration !
67 | ! -------------------------------------------!
68 |
69 | Hitot = Htot
70 |
71 | Gtot = 0.
72 | Lbtot = 0.
73 | Htot = 0.
74 | Q_rtot = 0.
75 | Q_ptot = 0.
76 |
77 | IF(i_dom_npart.gt.1) Livirt = s_Lvirt(:,:,:,1)
78 |
79 | ! print*, " (",pmm_rank,") Doing band integration"
80 |
81 | ! --------------------------!
82 | ! Integrate the narrow band !
83 | ! --------------------------!
84 |
85 | DO iquad=1,is_nkabs
86 |
87 | is_dirbeg = 1
88 | is_dirend = 0
89 |
90 | IF(is_cd.le.is_cf) THEN
91 | IF((iquad.ge.is_cd).and.(iquad.le.is_cf)) THEN
92 | is_dirend = is_ndir
93 | ENDIF
94 | ELSE
95 | IF (iquad.ge.is_cd) THEN
96 | is_dirend = is_ndir -1
97 | ELSEIF (iquad.le.is_cf) THEN
98 | is_dirbeg = 2
99 | is_dirend = is_ndir
100 | ENDIF
101 | ENDIF
102 |
103 | ! PRINT*, pmm_rank, iquad, is_dirbeg,is_dirend
104 |
105 | ! -------------------------------------------!
106 | ! Update values for reflection sub iteration !
107 | ! -------------------------------------------!
108 |
109 | Gi(iquad,:) = G(:)
110 |
111 | ! --------------------------------------------------!
112 | ! SPATIAL SCHEME !
113 | ! WSGG: Le traitement de chaque gaz se fait sur une !
114 | ! fraction du corps noir A*Lb (avec sum(A)=1) !
115 | ! --------------------------------------------------!
116 |
117 | CALL SPATIAL_SCHEME(I_SCHEME, Lb*WSGG_W_cell(iquad,:), &
118 | & Gi(iquad,:), G, kabs(iquad,:), kscat, &
119 | & Qr, H, alpha, Qp, &
120 | & WSGG_W_face(iquad,:)*Lo(:),1,iquad, &
121 | & i_dom_nthread, i_dom_npart)
122 |
123 | ! Special case for FSCK model
124 | ! WSGG_W_face(iquad,:)*Lo(:) -> s_WSGG_Wb(iquad,:)*Lo(:)
125 |
126 | ! -------------------------------------------------------------!
127 | ! Sum over slave quadrature point/directions for global models !
128 | ! -------------------------------------------------------------!
129 |
130 | Gtot(:) = Gtot(:) + G(:) *kabs(iquad,:)*wkabs(iquad) &
131 | & *DWVNB_SI
132 | Lbtot(:) = Lbtot(:) + Lb(:)*4*pi*kabs(iquad,:)*wkabs(iquad) &
133 | & *DWVNB_SI*WSGG_W_cell(iquad,:)
134 | Htot(:) = Htot(:) + H(:) *wkabs(iquad) &
135 | *DWVNB_SI
136 | DO i=1,3
137 | Q_rtot(i,:) = Q_rtot(i,:) + Qr(i,:)*wkabs(iquad)*DWVNB_SI
138 | Q_ptot(i,:) = Q_ptot(i,:) + Qp(i,:)*wkabs(iquad)*DWVNB_SI
139 | ENDDO
140 | ENDDO
141 |
142 | ! ----------------------------------------------!
143 | ! Convergence test for reflection or subdomains !
144 | ! ----------------------------------------------!
145 |
146 | IF(i_dom_npart.eq.1) THEN
147 |
148 | ! mloc = MAXLOC(abs(Htot-Hitot))
149 | ! IF (abs(Hitot(mloc(1))).ne.0.) THEN
150 | ! error1 = (abs(Htot(mloc(1))-Hitot(mloc(1))))/ &
151 | ! & abs(Hitot(mloc(1)))
152 | error2 = 0.
153 | IF(SUM(Hitot).ne.0) THEN
154 | error1 = abs(SUM(Hitot)-SUM(Htot))/SUM(Hitot)
155 | ELSEIF(SUM(ts_boundary%epsil)/is_nbfaces.eq.1) THEN
156 | error1 = 0.
157 | ELSE
158 | error1 = 1000.
159 | ENDIF
160 |
161 | ELSE
162 | ! mloc2 = MAXLOC(abs(s_Lvirt(:,:,:,1)-Livirt))
163 | ! IF (abs(Livirt(mloc2(1),mloc2(2),mloc2(3))).ne.0.) THEN
164 | ! error2 = (abs(s_Lvirt(mloc2(1),mloc2(2),mloc2(3),1)- &
165 | ! & Livirt(mloc2(1),mloc2(2),mloc2(3))) ) &
166 | ! & /abs(Livirt(mloc2(1),mloc2(2),mloc2(3)))
167 | IF (SUM(Livirt).ne.0.) THEN
168 | error2 = abs(SUM(s_Lvirt)-SUM(Livirt))/SUM(Livirt)
169 | ELSE
170 | error2 = 1000.
171 | ENDIF
172 | ENDIF
173 | error = MAX(error1,error2)
174 | IF(SUM(kabs(:,:)).eq.0) error = 0
175 |
176 | n_iter = n_iter + 1
177 |
178 | IF ((pmm_rank.eq.0).and.(error.ne.0)) &
179 | & WRITE(*,'(A14XI2XA7XF12.10)')" *** n_iter=",n_iter, &
180 | & " - err=", error
181 | IF((i_inter.eq.1).and.(pmm_rank.eq.0)) print*
182 |
183 | ! IF (pmm_rank.eq.0) print*, "mloc=",mloc(1),"- Htot=", &
184 | ! & Htot(mloc(1)),"- Hitot=",Hitot(mloc(1))
185 |
186 | ENDDO
187 |
188 | END SUBROUTINE BAND_INTEG
band_integ.F could be called by: