Line data Source code
1 : module fvm_overlap_mod
2 : use shr_kind_mod, only: r8=>shr_kind_r8
3 :
4 : real (kind=r8),parameter, private :: bignum = 1.0e20_r8
5 : real (kind=r8),parameter, private :: tiny = 1.0e-12_r8
6 : real (kind=r8),parameter, private :: fuzzy_width = 10.0_r8*tiny
7 :
8 : public:: compute_weights_cell
9 :
10 : private
11 : integer, parameter :: max_cross = 10
12 : contains
13 194400 : subroutine compute_weights_cell(nvertex,lexact_horizontal_line_integrals,&
14 194400 : xcell_in,ycell_in,jx,jy,nreconstruction,xgno,ygno,igno_min,igno_max,&
15 : jx_min, jx_max, jy_min, jy_max,&
16 194400 : ngauss,gauss_weights,abscissae,weights,weights_eul_index,jcollect,jmax_segments)
17 :
18 : implicit none
19 : integer , intent(in) :: nvertex
20 : logical, intent(in) :: lexact_horizontal_line_integrals
21 : integer , intent(in):: nreconstruction, jx,jy,ngauss,jmax_segments
22 : !
23 : ! dimension(nvertex)
24 : !
25 : real (kind=r8) , dimension(4), intent(in):: xcell_in,ycell_in
26 : !
27 : integer , intent(in) :: jx_min, jy_min, jx_max, jy_max,igno_min,igno_max
28 : !
29 : ! dimension(-ihalo:nc+2+ihalo)
30 : !
31 : real (kind=r8), dimension(igno_min:igno_max), intent(in) :: xgno, ygno
32 : !
33 : ! for Gaussian quadrature
34 : !
35 : real (kind=r8), dimension(:), intent(in) :: gauss_weights, abscissae !dimension(ngauss)
36 : !
37 : ! Number of Eulerian sub-cell integrals for the cell in question
38 : !
39 : integer , intent(out) :: jcollect
40 : !
41 : ! local workspace
42 : !
43 : !
44 : ! max number of line segments is:
45 : !
46 : ! (number of longitudes)*(max average number of crossings per line segment = 3)*ncube*2
47 : !
48 : real (kind=r8) , &
49 : dimension(jmax_segments,nreconstruction), intent(out) :: weights
50 : integer , &
51 : dimension(jmax_segments,2), intent(out) :: weights_eul_index
52 :
53 : integer :: jsegment
54 : !
55 : ! variables for registering crossings with Eulerian latitudes and longitudes
56 : !
57 : integer :: jcross_lat
58 : !
59 : ! max. crossings per side is 2*ihalo
60 : !
61 : real (kind=r8), &
62 : dimension(max_cross,2) :: r_cross_lat
63 : integer , &
64 : dimension(max_cross,2) :: cross_lat_eul_index
65 388800 : real (kind=r8) , dimension(nvertex) :: xcell,ycell
66 :
67 1166400 : xcell = xcell_in(1:nvertex)
68 972000 : ycell = ycell_in(1:nvertex)
69 :
70 194400 : jsegment = 0
71 59680800 : weights = 0.0_r8
72 194400 : jcross_lat = 0
73 :
74 : call side_integral(lexact_horizontal_line_integrals,xcell,ycell,nvertex,jsegment,jmax_segments,&
75 : weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,igno_min,igno_max,jx_min, jx_max, jy_min, jy_max,&
76 : ngauss,gauss_weights,abscissae,&
77 194400 : jcross_lat,r_cross_lat,cross_lat_eul_index)
78 : !
79 : !**********************
80 : !
81 : ! Do inner integrals
82 : !
83 : !**********************
84 : !
85 : call compute_inner_line_integrals_lat(lexact_horizontal_line_integrals,&
86 : r_cross_lat,cross_lat_eul_index,&
87 : jcross_lat,jsegment,xgno,igno_min,igno_max,jx_min, jx_max, jy_min, jy_max,&
88 : weights,weights_eul_index,&
89 194400 : nreconstruction,ngauss,gauss_weights,abscissae)
90 :
91 194400 : IF (ABS((jcross_lat/2)-DBLE(jcross_lat)/2.0_r8)>tiny) then
92 0 : WRITE(*,*) "number of latitude crossings are not even: ABORT",jcross_lat,jx,jy
93 0 : STOP
94 : END IF
95 :
96 : !
97 : ! collect line-segment that reside in the same Eulerian cell
98 : !
99 194400 : if (jsegment>0) then
100 194400 : call collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments)
101 : else
102 0 : jcollect = 0
103 : end if
104 194400 : end subroutine compute_weights_cell
105 : !
106 : !****************************************************************************
107 : !
108 : ! organize data and store it
109 : !
110 : !****************************************************************************
111 : !
112 194400 : subroutine collect(weights,weights_eul_index,nreconstruction,jcollect,jsegment,jmax_segments)
113 : implicit none
114 : integer , INTENT(IN ) :: jsegment,jmax_segments
115 : integer , intent(in) :: nreconstruction
116 : !
117 : real (kind=r8) , dimension(:,:), intent(inout) :: weights !dimension(jmax_segments,nreconstruction)
118 : integer , dimension(:,:), intent(inout) :: weights_eul_index !dimension(jmax_segments,2)
119 : integer , INTENT(OUT ) :: jcollect
120 : !
121 : ! local workspace
122 : !
123 : integer :: imin, imax, jmin, jmax, i,j,k,h
124 : logical :: ltmp
125 :
126 388800 : real (kind=r8) , dimension(jmax_segments,nreconstruction) :: weights_out
127 388800 : integer , dimension(jmax_segments,2 ) :: weights_eul_index_out
128 :
129 59680800 : weights_out = 0.0_r8
130 20023200 : weights_eul_index_out = -100
131 :
132 1166400 : imin = MINVAL(weights_eul_index(1:jsegment,1))
133 972000 : imax = MAXVAL(weights_eul_index(1:jsegment,1))
134 972000 : jmin = MINVAL(weights_eul_index(1:jsegment,2))
135 972000 : jmax = MAXVAL(weights_eul_index(1:jsegment,2))
136 :
137 194400 : ltmp = .FALSE.
138 :
139 194400 : jcollect = 1
140 :
141 388800 : do j=jmin,jmax
142 583200 : do i=imin,imax
143 972000 : do k=1,jsegment
144 972000 : if (weights_eul_index(k,1)==i.AND.weights_eul_index(k,2)==j) then
145 : weights_out(jcollect,1:nreconstruction) = &
146 5443200 : weights_out(jcollect,1:nreconstruction) + weights(k,1:nreconstruction)
147 : ltmp = .TRUE.
148 : h = k
149 : endif
150 : enddo
151 194400 : if (ltmp) then
152 583200 : weights_eul_index_out(jcollect,:) = weights_eul_index(h,:)
153 194400 : jcollect = jcollect+1
154 : endif
155 388800 : ltmp = .FALSE.
156 : enddo
157 : enddo
158 194400 : jcollect = jcollect-1
159 59875200 : weights = weights_out
160 20217600 : weights_eul_index = weights_eul_index_out
161 194400 : end subroutine collect
162 : !
163 : !*****************************************************************************************
164 : !
165 : ! compute crossings with Eulerian latitudes and longitudes
166 : !
167 : !*****************************************************************************************
168 : !
169 583200 : subroutine compute_inner_line_integrals_lat(lexact_horizontal_line_integrals,r_cross_lat,&
170 194400 : cross_lat_eul_index,&
171 194400 : jcross_lat,jsegment,xgno,igno_min,igno_max,jx_min,jx_max,jy_min, jy_max,weights,weights_eul_index,&
172 194400 : nreconstruction,ngauss,gauss_weights,abscissae)
173 : implicit none
174 : logical, intent(in) :: lexact_horizontal_line_integrals
175 : !
176 : ! variables for registering crossings with Eulerian latitudes and longitudes
177 : !
178 : integer , intent(in):: jcross_lat, nreconstruction,ngauss,igno_min,igno_max
179 : integer , intent(inout):: jsegment
180 : !
181 : ! for Gaussian quadrature
182 : !
183 : real (kind=r8), dimension(ngauss), intent(in) :: gauss_weights, abscissae
184 : !
185 : ! max. crossings per side is 2*ihalo
186 : !
187 :
188 : real (kind=r8) , dimension(:,:), intent(in):: r_cross_lat ! dimension(8*ihalo,2)
189 : integer , dimension(:,:), intent(in):: cross_lat_eul_index ! ! dimension(8*ihalo,2)
190 : integer , intent(in):: jx_min, jx_max, jy_min, jy_max
191 :
192 : real (kind=r8), dimension(igno_min:igno_max), intent(in) :: xgno !dimension(-ihalo:nc+2+ihalo)
193 : !
194 : ! dimension(jmax_segments,nreconstruction)
195 : !
196 : real (kind=r8), dimension(:,:), intent(inout) :: weights
197 : !
198 : ! dimension(jmax_segments,2)
199 : !
200 : integer , dimension(:,:), intent(inout) :: weights_eul_index
201 :
202 388800 : real (kind=r8) , dimension(nreconstruction) :: weights_tmp
203 : integer :: imin,imax,i,j,k,h
204 : real (kind=r8), dimension(2) :: rstart,rend,rend_tmp
205 : real (kind=r8), dimension(2) :: xseg, yseg
206 : 5 FORMAT(10e14.6)
207 194400 : if (jcross_lat>0) then
208 0 : do i=MINVAL(cross_lat_eul_index(1:jcross_lat,2)),MAXVAL(cross_lat_eul_index(1:jcross_lat,2))
209 : !
210 : ! find "first" crossing with Eulerian cell i
211 : !
212 0 : do k=1,jcross_lat
213 0 : if (cross_lat_eul_index(k,2)==i) exit
214 : enddo
215 0 : do j=k+1,jcross_lat
216 : !
217 : ! find "second" crossing with Eulerian cell i
218 : !
219 0 : if (cross_lat_eul_index(j,2)==i) then
220 0 : if (r_cross_lat(k,1)<r_cross_lat(j,1)) then
221 0 : rstart = r_cross_lat(k,1:2)
222 0 : rend = r_cross_lat(j,1:2)
223 0 : imin = cross_lat_eul_index(k,1)
224 0 : imax = cross_lat_eul_index(j,1)
225 : else
226 0 : rstart = r_cross_lat(j,1:2)
227 0 : rend = r_cross_lat(k,1:2)
228 0 : imin = cross_lat_eul_index(j,1)
229 0 : imax = cross_lat_eul_index(k,1)
230 : endif
231 0 : do h=imin,imax
232 0 : if (h==imax) then
233 0 : rend_tmp = rend
234 : else
235 0 : rend_tmp(1) = xgno(h+1)
236 0 : rend_tmp(2) = r_cross_lat(k,2)
237 : endif
238 0 : xseg(1) = rstart(1)
239 0 : xseg(2) = rend_tmp(1)
240 0 : yseg(1) = rstart(2)
241 0 : yseg(2) = rend_tmp(2)
242 : call get_weights_exact(lexact_horizontal_line_integrals, weights_tmp,xseg,yseg,&
243 0 : nreconstruction,ngauss,gauss_weights,abscissae)
244 :
245 0 : if (i.LE.jy_max-1.AND.i.GE.jy_min.AND.h.LE.jx_max-1.AND.h.GE.jx_min) then
246 0 : jsegment=jsegment+1
247 0 : weights_eul_index(jsegment,1) = h
248 0 : weights_eul_index(jsegment,2) = i
249 0 : weights(jsegment,1:nreconstruction) = -weights_tmp
250 : endif
251 : !
252 : ! subtract the same weights on the west side of the line
253 : !
254 0 : if (i.LE.jy_max.AND.i.GE.jy_min+1.AND.h.LE.jx_max-1.AND.h.GE.jx_min) then
255 0 : jsegment = jsegment+1
256 0 : weights_eul_index(jsegment,1) = h
257 0 : weights_eul_index(jsegment,2) = i-1
258 0 : weights(jsegment,1:nreconstruction) = weights_tmp
259 : endif
260 : !
261 : ! prepare for next iteration
262 : !
263 0 : rstart = rend_tmp
264 : enddo
265 : endif
266 : enddo
267 : enddo
268 : endif
269 194400 : end subroutine compute_inner_line_integrals_lat
270 :
271 : !
272 : ! line integral from (a1_in,a2_in) to (b1_in,b2_in)
273 : ! If line is coniciding with an Eulerian longitude or latitude the routine
274 : ! needs to know where an adjacent side is located to determine which
275 : ! reconstruction must be used. therefore (c1,c2) is passed to the routine
276 : !
277 : !
278 :
279 194400 : subroutine side_integral(lexact_horizontal_line_integrals,&
280 194400 : x_in,y_in,nvertex,jsegment,jmax_segments,&
281 194400 : weights,weights_eul_index,nreconstruction,jx,jy,xgno,ygno,igno_min,igno_max,&
282 : jx_min,jx_max,jy_min,jy_max,&
283 194400 : ngauss,gauss_weights,abscissae,&!)!phl add jx_min etc.
284 : jcross_lat,r_cross_lat,cross_lat_eul_index)
285 : implicit none
286 :
287 :
288 : logical, intent(in) :: lexact_horizontal_line_integrals
289 : integer , intent(in) :: nreconstruction,jx,jy,jmax_segments,ngauss
290 : integer , intent(in) :: nvertex,igno_min,igno_max
291 : !
292 : ! for Gaussian quadrature
293 : !
294 : real (kind=r8), dimension(:), intent(in) :: gauss_weights, abscissae !dimension(ngauss)
295 : real (kind=r8), dimension(:) , intent(in) :: x_in,y_in !dimension(1:nvertex)
296 :
297 : integer , intent(in) :: jx_min, jy_min, jx_max, jy_max
298 : real (kind=r8), dimension(igno_min:igno_max), intent(in) :: xgno, ygno !dimension(-ihalo:nc+2+ihalo)
299 : integer , intent(inout) :: jsegment
300 : ! integer ,dimension(0:2),intent(in) :: jx_eul_in, jy_eul_in
301 : real (kind=r8) , dimension(:,:), intent(out) :: weights !dimension(jmax_segments,nreconstruction)
302 : integer , &
303 : dimension(jmax_segments,2), intent(out) :: weights_eul_index
304 :
305 : !
306 : ! variables for registering crossings with Eulerian latitudes and longitudes
307 : !
308 : integer , intent(inout):: jcross_lat
309 : !
310 : ! max. crossings per side is 2*ihalo
311 : !
312 : real (kind=r8), &
313 : dimension(max_cross,2), intent(inout):: r_cross_lat
314 : integer , &
315 : dimension(max_cross,2), intent(inout):: cross_lat_eul_index
316 : !
317 : ! local variables
318 : !
319 : real (kind=r8), dimension(2) :: xseg,yseg
320 : real (kind=r8), dimension(0:3) :: x,y
321 : real (kind=r8) :: xeul,yeul,xcross,ycross,slope
322 : integer :: jx_eul_tmp,jy_eul_tmp
323 : integer :: xsgn1,ysgn1,xsgn2,ysgn2
324 : integer :: iter
325 : logical :: lcontinue, lsame_cell_x, lsame_cell_y
326 :
327 : integer :: jx_eul, jy_eul, side_count
328 388800 : real (kind=r8), dimension(0:nvertex+2) :: xcell,ycell
329 :
330 :
331 : 5 FORMAT(10e14.6)
332 : !
333 : !***********************************************
334 : !
335 : ! find jx_eul and jy_eul for (x(1),y(1))
336 : !
337 : !***********************************************
338 : !
339 194400 : jx_eul = jx; jy_eul = jy
340 1749600 : xcell(1:nvertex)=x_in; ycell(1:nvertex)=y_in
341 972000 : DO iter=1,nvertex
342 777600 : CALL truncate_vertex(xcell(iter),jx_eul,xgno,igno_min,igno_max)
343 972000 : CALL truncate_vertex(ycell(iter),jy_eul,ygno,igno_min,igno_max)
344 : END DO
345 194400 : xcell(0) = xcell(nvertex); xcell(nvertex+1)=xcell(1); xcell(nvertex+2)=xcell(2);
346 194400 : ycell(0) = ycell(nvertex); ycell(nvertex+1)=ycell(1); ycell(nvertex+2)=ycell(2);
347 :
348 :
349 7192800 : IF ((&
350 388800 : MAXVAL(xcell).LE.xgno(jx_min).OR.MINVAL(xcell).GE.xgno(jx_max).OR.&
351 388800 : MAXVAL(ycell).LE.ygno(jy_min).OR.MINVAL(ycell).GE.ygno(jy_max))) THEN
352 : !
353 : ! entire cell off panel
354 : !
355 : ELSE
356 194400 : jx_eul = MIN(MAX(jx,jx_min),jx_max-1)
357 194400 : jy_eul = MIN(MAX(jy,jy_min),jy_max-1)
358 194400 : CALL which_eul_cell(xcell(1:3),jx_eul,xgno,igno_min,igno_max)
359 194400 : CALL which_eul_cell(ycell(1:3),jy_eul,ygno,igno_min,igno_max)
360 :
361 194400 : side_count = 1
362 972000 : DO WHILE (side_count<nvertex+1)
363 777600 : iter = 0
364 777600 : lcontinue = .TRUE.
365 6998400 : x(0:3) = xcell(side_count-1:side_count+2); y(0:3) = ycell(side_count-1:side_count+2);
366 1555200 : DO while (lcontinue)
367 777600 : iter = iter+1
368 777600 : IF (iter>10) THEN
369 0 : WRITE(*,*) "search not converging",iter
370 0 : STOP
371 : END IF
372 777600 : lsame_cell_x = (x(2).GE.xgno(jx_eul).AND.x(2).LE.xgno(jx_eul+1))
373 777600 : lsame_cell_y = (y(2).GE.ygno(jy_eul).AND.y(2).LE.ygno(jy_eul+1))
374 777600 : IF (lsame_cell_x.AND.lsame_cell_y) THEN
375 : !
376 : !****************************
377 : !
378 : ! same cell integral
379 : !
380 : !****************************
381 : !
382 777600 : xseg(1) = x(1); yseg(1) = y(1); xseg(2) = x(2); yseg(2) = y(2)
383 777600 : jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
384 777600 : lcontinue = .FALSE.
385 : !
386 : ! prepare for next side if (x(2),y(2)) is on a grid line
387 : !
388 777600 : IF (x(2).EQ.xgno(jx_eul+1).AND.x(3)>xgno(jx_eul+1)) THEN
389 : !
390 : ! cross longitude jx_eul+1
391 : !
392 0 : jx_eul=jx_eul+1
393 777600 : ELSE IF (x(2).EQ.xgno(jx_eul ).AND.x(3)<xgno(jx_eul)) THEN
394 : !
395 : ! cross longitude jx_eul
396 : !
397 0 : jx_eul=jx_eul-1
398 : END IF
399 777600 : IF (y(2).EQ.ygno(jy_eul+1).AND.y(3)>ygno(jy_eul+1)) THEN
400 : !
401 : ! register crossing with latitude: line-segments point Northward
402 : !
403 0 : jcross_lat = jcross_lat + 1
404 0 : jy_eul = jy_eul + 1
405 0 : cross_lat_eul_index(jcross_lat,1) = jx_eul
406 0 : cross_lat_eul_index(jcross_lat,2) = jy_eul
407 0 : r_cross_lat(jcross_lat,1) = x(2)
408 0 : r_cross_lat(jcross_lat,2) = y(2)
409 : ! write(*,*) "A register crossing with latitude",x(2),y(2),jx_eul,jy_eul
410 777600 : ELSE IF (y(2).EQ.ygno(jy_eul ).AND.y(3)<ygno(jy_eul)) THEN
411 : !
412 : ! register crossing with latitude: line-segments point Southward
413 : !
414 0 : jcross_lat = jcross_lat+1
415 0 : cross_lat_eul_index(jcross_lat,1) = jx_eul
416 0 : cross_lat_eul_index(jcross_lat,2) = jy_eul
417 0 : r_cross_lat(jcross_lat,1) = x(2)
418 0 : r_cross_lat(jcross_lat,2) = y(2)
419 : ! write(*,*) "B register crossing with latitude",x(2),y(2),jx_eul,jy_eul
420 : !
421 0 : jy_eul=jy_eul-1
422 : END IF
423 : lcontinue=.FALSE.
424 : ELSE
425 : !
426 : !****************************
427 : !
428 : ! not same cell integral
429 : !
430 : !****************************
431 : !
432 0 : IF (lsame_cell_x) THEN
433 0 : ysgn1 = (1+INT(SIGN(1.0_r8,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0"
434 0 : ysgn2 = INT(SIGN(1.0_r8,y(2)-y(1))) !"1" if y(2)>y(1) else "-1"
435 : !
436 : !*******************************************************************************
437 : !
438 : ! there is at least one crossing with latitudes but no crossing with longitudes
439 : !
440 : !*******************************************************************************
441 : !
442 0 : yeul = ygno(jy_eul+ysgn1)
443 0 : IF (x(1).EQ.x(2)) THEN
444 : !
445 : ! line segment is parallel to longitude (infinite slope)
446 : !
447 : xcross = x(1)
448 : ELSE
449 0 : slope = (y(2)-y(1))/(x(2)-x(1))
450 0 : xcross = x_cross_eul_lat(x(1),y(1),yeul,slope)
451 : !
452 : ! constrain crossing to be "physically" possible
453 : !
454 0 : xcross = MIN(MAX(xcross,xgno(jx_eul)),xgno(jx_eul+1))
455 : !
456 : ! debugging
457 : !
458 0 : IF (xcross.GT.xgno(jx_eul+1).OR.xcross.LT.xgno(jx_eul)) THEN
459 0 : WRITE(*,*) "xcross is out of range",jx,jy
460 0 : WRITE(*,*) "xcross-xgno(jx_eul+1), xcross-xgno(jx_eul))",&
461 0 : xcross-xgno(jx_eul+1), xcross-ygno(jx_eul)
462 0 : STOP
463 : END IF
464 : END IF
465 0 : xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul
466 0 : jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
467 : !
468 : ! prepare for next iteration
469 : !
470 0 : x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2
471 : !
472 : ! register crossing with latitude
473 : !
474 0 : jcross_lat = jcross_lat+1
475 0 : cross_lat_eul_index(jcross_lat,1) = jx_eul
476 0 : if (ysgn2>0) then
477 0 : cross_lat_eul_index(jcross_lat,2) = jy_eul
478 : else
479 0 : cross_lat_eul_index(jcross_lat,2) = jy_eul+1
480 : end if
481 0 : r_cross_lat(jcross_lat,1) = xcross
482 0 : r_cross_lat(jcross_lat,2) = yeul
483 0 : ELSE IF (lsame_cell_y) THEN
484 : !
485 : !*******************************************************************************
486 : !
487 : ! there is at least one crossing with longitudes but no crossing with latitudes
488 : !
489 : !*******************************************************************************
490 : !
491 0 : xsgn1 = (1+INT(SIGN(1.0_r8,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0"
492 0 : xsgn2 = INT(SIGN(1.0_r8,x(2)-x(1))) !"1" if x(2)>x(1) else "-1"
493 0 : xeul = xgno(jx_eul+xsgn1)
494 0 : IF (ABS(x(2)-x(1))<fuzzy_width) THEN
495 : ! fuzzy crossing
496 0 : ycross = 0.5_r8*(y(2)-y(1))
497 : ELSE
498 0 : slope = (y(2)-y(1))/(x(2)-x(1))
499 0 : ycross = y_cross_eul_lon(x(1),y(1),xeul,slope)
500 : END IF
501 : !
502 : ! constrain crossing to be "physically" possible
503 : !
504 0 : ycross = MIN(MAX(ycross,ygno(jy_eul)),ygno(jy_eul+1))
505 : !
506 : ! debugging
507 : !
508 0 : IF (ycross.GT.ygno(jy_eul+1).OR.ycross.LT.ygno(jy_eul)) THEN
509 0 : WRITE(*,*) "ycross is out of range"
510 0 : WRITE(*,*) "jx,jy,jx_eul,jy_eul",jx,jy,jx_eul,jy_eul
511 0 : WRITE(*,*) "ycross-ygno(jy_eul+1), ycross-ygno(jy_eul))",&
512 0 : ycross-ygno(jy_eul+1), ycross-ygno(jy_eul)
513 0 : STOP
514 : END IF
515 0 : xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross
516 0 : jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
517 : !
518 : ! prepare for next iteration
519 : !
520 0 : x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2
521 : ELSE
522 : !
523 : !*******************************************************************************
524 : !
525 : ! there are crossings with longitude(s) and latitude(s)
526 : !
527 : !*******************************************************************************
528 : !
529 0 : xsgn1 = (1+INT(SIGN(1.0_r8,x(2)-x(1))))/2 !"1" if x(2)>x(1) else "0"
530 0 : xsgn2 = (INT(SIGN(1.0_r8,x(2)-x(1)))) !"1" if x(2)>x(1) else "0"
531 0 : xeul = xgno(jx_eul+xsgn1)
532 0 : ysgn1 = (1+INT(SIGN(1.0_r8,y(2)-y(1))))/2 !"1" if y(2)>y(1) else "0"
533 0 : ysgn2 = INT(SIGN(1.0_r8,y(2)-y(1))) !"1" if y(2)>y(1) else "-1"
534 0 : yeul = ygno(jy_eul+ysgn1)
535 :
536 0 : slope = (y(2)-y(1))/(x(2)-x(1))
537 0 : IF (ABS(x(2)-x(1))<fuzzy_width) THEN
538 0 : ycross = 0.5_r8*(y(2)-y(1))
539 : ELSE
540 0 : ycross = y_cross_eul_lon(x(1),y(1),xeul,slope)
541 : END IF
542 0 : xcross = x_cross_eul_lat(x(1),y(1),yeul,slope)
543 :
544 :
545 0 : IF ((xsgn2>0.AND.xcross.LE.xeul).OR.(xsgn2<0.AND.xcross.GE.xeul)) THEN
546 : !
547 : ! cross latitude
548 : !
549 0 : xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xcross; yseg(2) = yeul
550 0 : jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
551 : !
552 : ! prepare for next iteration
553 : !
554 0 : x(0) = x(1); y(0) = y(1); x(1) = xcross; y(1) = yeul; jy_eul = jy_eul+ysgn2
555 : !
556 : ! register crossing with latitude
557 : !
558 0 : jcross_lat = jcross_lat+1
559 0 : cross_lat_eul_index(jcross_lat,1) = jx_eul
560 0 : if (ysgn2>0) then
561 0 : cross_lat_eul_index(jcross_lat,2) = jy_eul
562 : else
563 0 : cross_lat_eul_index(jcross_lat,2) = jy_eul+1
564 : end if
565 0 : r_cross_lat(jcross_lat,1) = xcross
566 0 : r_cross_lat(jcross_lat,2) = yeul
567 : ! write(*,*) "D register crossing with latitude",xcross,yeul,jx_eul,cross_lat_eul_index(jcross_lat,2)
568 : ELSE
569 : !
570 : ! cross longitude
571 : !
572 0 : xseg(1) = x(1); yseg(1) = y(1); xseg(2) = xeul; yseg(2) = ycross
573 0 : jx_eul_tmp = jx_eul; jy_eul_tmp = jy_eul;
574 : !
575 : ! prepare for next iteration
576 : !
577 0 : x(0) = x(1); y(0) = y(1); x(1) = xeul; y(1) = ycross; jx_eul = jx_eul+xsgn2
578 : END IF
579 :
580 : END IF
581 : END IF
582 : !
583 : ! register line-segment (don't register line-segment if outside of panel)
584 : !
585 : if (jx_eul_tmp>=jx_min.AND.jy_eul_tmp>=jy_min.AND.&
586 1555200 : jx_eul_tmp<=jx_max-1.AND.jy_eul_tmp<=jy_max-1) then
587 777600 : jsegment=jsegment+1
588 777600 : weights_eul_index(jsegment,1) = jx_eul_tmp
589 777600 : weights_eul_index(jsegment,2) = jy_eul_tmp
590 :
591 : call get_weights_exact(lexact_horizontal_line_integrals.AND.ABS(yseg(2)-yseg(1))<tiny,&
592 : weights(jsegment,:),&
593 777600 : xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
594 : !old call get_weights_gauss(weights(jsegment,1:nreconstruction),&
595 : !old xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
596 : ELSE
597 : !
598 : ! segment outside of panel
599 : !
600 : END IF
601 :
602 : END DO
603 777600 : side_count = side_count+1
604 : END DO
605 : END IF
606 194400 : end subroutine side_integral
607 :
608 :
609 0 : real (kind=r8) function y_cross_eul_lon(x,y,xeul,slope)
610 : implicit none
611 : real (kind=r8), intent(in) :: x,y
612 : real (kind=r8) , intent(in) :: xeul,slope
613 : !
614 : ! line: y=a*x+b
615 : !
616 : real (kind=r8) :: b
617 :
618 0 : b = y-slope*x
619 0 : y_cross_eul_lon = slope*xeul+b
620 0 : end function y_cross_eul_lon
621 :
622 0 : real (kind=r8) function x_cross_eul_lat(x,y,yeul,slope)
623 : implicit none
624 : real (kind=r8), intent(in) :: x,y
625 : real (kind=r8) , intent(in) :: yeul,slope
626 :
627 0 : if (fuzzy(ABS(slope),fuzzy_width)>0) THEN
628 0 : x_cross_eul_lat = x+(yeul-y)/slope
629 : ELSE
630 : x_cross_eul_lat = bignum
631 : END IF
632 0 : end function x_cross_eul_lat
633 :
634 777600 : subroutine get_weights_exact(lexact_horizontal_line_integrals,weights,xseg,yseg,nreconstruction,&
635 777600 : ngauss,gauss_weights,abscissae)
636 : use fvm_analytic_mod, only: I_00, I_10, I_01, I_20, I_02, I_11
637 : implicit none
638 : logical, intent(in) :: lexact_horizontal_line_integrals
639 : integer , intent(in) :: nreconstruction, ngauss
640 : real (kind=r8), intent(out) :: weights(:)
641 : real (kind=r8), dimension(:), intent(in) :: gauss_weights, abscissae !dimension(ngauss)
642 :
643 :
644 : real (kind=r8), dimension(:), intent(in) :: xseg,yseg !dimension(2)
645 : !
646 : ! compute weights
647 : !
648 777600 : if(lexact_horizontal_line_integrals) then
649 388800 : weights(1) = ((I_00(xseg(2),yseg(2))-I_00(xseg(1),yseg(1))))
650 388800 : if (ABS(weights(1))>1.0_r8) THEN
651 0 : WRITE(*,*) "1 exact weights(jsegment)",weights(1),xseg,yseg
652 0 : stop
653 : end if
654 388800 : if (nreconstruction>1) then
655 388800 : weights(2) = ((I_10(xseg(2),yseg(2))-I_10(xseg(1),yseg(1))))
656 388800 : weights(3) = ((I_01(xseg(2),yseg(2))-I_01(xseg(1),yseg(1))))
657 : endif
658 388800 : if (nreconstruction>3) then
659 388800 : weights(4) = ((I_20(xseg(2),yseg(2))-I_20(xseg(1),yseg(1))))
660 388800 : weights(5) = ((I_02(xseg(2),yseg(2))-I_02(xseg(1),yseg(1))))
661 388800 : weights(6) = ((I_11(xseg(2),yseg(2))-I_11(xseg(1),yseg(1))))
662 : endif
663 : else
664 388800 : call get_weights_gauss(weights,xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
665 : endif
666 777600 : end subroutine get_weights_exact
667 :
668 :
669 :
670 388800 : subroutine get_weights_gauss(weights,xseg,yseg,nreconstruction,ngauss,gauss_weights,abscissae)
671 : use fvm_analytic_mod, only: F_00, F_10, F_01, F_20, F_02, F_11
672 : implicit none
673 : integer , intent(in) :: nreconstruction,ngauss
674 : real (kind=r8), intent(out) :: weights(:)
675 : real (kind=r8), dimension(2 ), intent(in) :: xseg,yseg
676 : real (kind=r8) :: slope
677 : !
678 : ! compute weights
679 : !
680 : !
681 : ! for Gaussian quadrature
682 : !
683 : real (kind=r8), dimension(ngauss), intent(in) :: gauss_weights, abscissae
684 :
685 : ! if line-segment parallel to x or y use exact formulaes else use qudrature
686 : !
687 : real (kind=r8) :: b,integral,dx2,xc,x,y
688 : integer :: i
689 :
690 : ! if (fuzzy(abs(xseg(1) -xseg(2)),fuzzy_width)==0)then
691 388800 : if (xseg(1).EQ.xseg(2))then
692 2721600 : weights = 0.0_r8
693 : else
694 0 : slope = (yseg(2)-yseg(1))/(xseg(2)-xseg(1))
695 0 : b = yseg(1)-slope*xseg(1)
696 0 : dx2 = 0.5_r8*(xseg(2)-xseg(1))
697 0 : xc = 0.5_r8*(xseg(1)+xseg(2))
698 0 : integral = 0.0_r8
699 0 : do i=1,ngauss
700 0 : x = xc+abscissae(i)*dx2
701 0 : y = slope*x+b
702 0 : integral = integral+gauss_weights(i)*F_00(x,y)
703 : enddo
704 0 : weights(1) = integral*dx2
705 0 : if (nreconstruction>1) then
706 : integral = 0.0_r8
707 0 : do i=1,ngauss
708 0 : x = xc+abscissae(i)*dx2
709 0 : y = slope*x+b
710 0 : integral = integral+gauss_weights(i)*F_10(x,y)
711 : enddo
712 0 : weights(2) = integral*dx2
713 0 : integral = 0.0_r8
714 0 : do i=1,ngauss
715 0 : x = xc+abscissae(i)*dx2
716 0 : y = slope*x+b
717 0 : integral = integral+gauss_weights(i)*F_01(x,y)
718 : enddo
719 0 : weights(3) = integral*dx2
720 : endif
721 0 : if (nreconstruction>3) then
722 : integral = 0.0_r8
723 0 : do i=1,ngauss
724 0 : x = xc+abscissae(i)*dx2
725 0 : y = slope*x+b
726 0 : integral = integral+gauss_weights(i)*F_20(x,y)
727 : enddo
728 0 : weights(4) = integral*dx2
729 0 : integral = 0.0_r8
730 0 : do i=1,ngauss
731 0 : x = xc+abscissae(i)*dx2
732 0 : y = slope*x+b
733 0 : integral = integral+gauss_weights(i)*F_02(x,y)
734 : enddo
735 0 : weights(5) = integral*dx2
736 0 : integral = 0.0_r8
737 0 : do i=1,ngauss
738 0 : x = xc+abscissae(i)*dx2
739 0 : y = slope*x+b
740 0 : integral = integral+gauss_weights(i)*F_11(x,y)
741 : enddo
742 0 : weights(6) = integral*dx2
743 : endif
744 : end if
745 388800 : end subroutine get_weights_gauss
746 :
747 1555200 : subroutine truncate_vertex(x,j_eul,gno,igno_min,igno_max)
748 : implicit none
749 : integer , intent(inout) :: j_eul
750 : integer , intent(in) :: igno_min,igno_max
751 :
752 : real (kind=r8) , intent(inout) :: x
753 : real (kind=r8), dimension(igno_min:igno_max), intent(in) :: gno !dimension(-ihalo:nc+2+ihalo)
754 : ! real (kind=r8), intent(in) :: eps
755 :
756 : logical :: lcontinue
757 : integer :: iter, xsgn
758 : real (kind=r8) :: dist,dist_new,tmp
759 :
760 1555200 : lcontinue = .TRUE.
761 1555200 : iter = 0
762 1555200 : dist = bignum
763 :
764 1555200 : xsgn = INT(SIGN(1.0_r8,x-gno(j_eul)))
765 :
766 3693600 : DO WHILE (lcontinue)
767 2138400 : if ((j_eul<igno_min) .or. (j_eul>igno_max)) then
768 0 : write(*,*) 'something is wrong', j_eul,igno_min,igno_max, iter
769 0 : stop
770 : endif
771 2138400 : iter = iter+1
772 2138400 : tmp = x-gno(j_eul)
773 2138400 : dist_new = ABS(tmp)
774 2138400 : IF (dist_new>dist) THEN
775 : lcontinue = .FALSE.
776 2138400 : ELSE IF (ABS(tmp)<1.0E-9_r8) THEN
777 1555200 : x = gno(j_eul)
778 1555200 : lcontinue = .FALSE.
779 : ELSE
780 583200 : j_eul = j_eul+xsgn
781 583200 : dist = dist_new
782 : END IF
783 3693600 : IF (iter>100) THEN
784 0 : WRITE(*,*) "truncate vertex not converging"
785 0 : STOP
786 : END IF
787 : END DO
788 1555200 : END subroutine truncate_vertex
789 :
790 388800 : subroutine which_eul_cell(x,j_eul,gno,igno_min,igno_max)
791 : implicit none
792 : integer , intent(inout) :: j_eul
793 : integer , intent(in) :: igno_min,igno_max
794 : real (kind=r8), dimension(:) , intent(in):: x !dimension(3)
795 : real (kind=r8), dimension(igno_min:igno_max), intent(in) :: gno ! dimension(-ihalo:nc+2+ihalo)
796 :
797 : logical :: lcontinue
798 : integer :: iter
799 :
800 388800 : lcontinue = .TRUE.
801 388800 : iter = 0
802 :
803 777600 : DO WHILE (lcontinue)
804 388800 : iter = iter+1
805 388800 : IF (x(1).GE.gno(j_eul).AND.x(1).LT.gno(j_eul+1)) THEN
806 388800 : lcontinue = .FALSE.
807 : !
808 : ! special case when x(1) is on top of grid line
809 : !
810 388800 : IF (x(1).EQ.gno(j_eul)) THEN
811 : !
812 : ! x(1) is on top of gno(J_eul)
813 : !
814 388800 : IF (x(2).GT.gno(j_eul)) THEN
815 : j_eul = j_eul
816 194400 : ELSE IF (x(2).LT.gno(j_eul)) THEN
817 0 : j_eul = j_eul-1
818 : ELSE
819 : !
820 : ! x(2) is on gno(j_eul) grid line; need x(3) to determine Eulerian cell
821 : !
822 194400 : IF (x(3).GT.gno(j_eul)) THEN
823 : !
824 : ! x(3) to the right
825 : !
826 : j_eul = j_eul
827 0 : ELSE IF (x(3).LT.gno(j_eul)) THEN
828 : !
829 : ! x(3) to the left
830 : !
831 0 : j_eul = j_eul-1
832 : ELSE
833 0 : WRITE(*,*) "inconsistent cell: x(1)=x(2)=x(3)",x(1),x(2),x(3)
834 0 : STOP
835 : END IF
836 : END IF
837 : END IF
838 : ELSE
839 : !
840 : ! searching - prepare for next iteration
841 : !
842 0 : IF (x(1).GE.gno(j_eul+1)) THEN
843 0 : j_eul = j_eul + 1
844 : ELSE
845 : !
846 : ! x(1).LT.gno(j_eul)
847 : !
848 0 : j_eul = j_eul - 1
849 : END IF
850 : END IF
851 777600 : IF (iter>1000.OR.j_eul<igno_min.OR.j_eul>igno_max) THEN
852 0 : WRITE(*,*) "search is which_eul_cell not converging!", iter, j_eul,igno_min,igno_max
853 0 : WRITE(*,*) "gno", gno(igno_min), gno(igno_max)
854 0 : write(*,*) gno
855 0 : STOP
856 : END IF
857 : END DO
858 388800 : END subroutine which_eul_cell
859 :
860 :
861 0 : function fuzzy(x,epsilon)
862 : implicit none
863 :
864 : integer :: fuzzy
865 : real (kind=r8), intent(in) :: epsilon
866 : real (kind=r8) :: x
867 :
868 0 : IF (ABS(x)<epsilon) THEN
869 : fuzzy = 0
870 0 : ELSE IF (x >epsilon) THEN
871 : fuzzy = 1
872 : ELSE !IF (x < fuzzy_width) THEN
873 0 : fuzzy = -1
874 : ENDIF
875 0 : end function
876 :
877 : end module fvm_overlap_mod
|
