PlasCom2  1.0
XPACC Multi-physics simluation application
FD1.f90
Go to the documentation of this file.
1 SUBROUTINE fd1load
2  WRITE(*,*) 'FD1 Kernel library loaded'
3 
4 END SUBROUTINE fd1load
5 
6 SUBROUTINE fd1interior(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX, U, KN)
7 
8  IMPLICIT NONE
9 
10  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
11  REAL(KIND=8), INTENT(IN) :: DX
12  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
13  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
14 
15  INTEGER(KIND=8) :: I, J, K, LI, LIK, LIJ, IS1, IE1
16  REAL(KIND=8) :: fac
17  INTEGER(KIND=8) :: plane
18 
19  is1 = s1
20  ie1 = e1
21  IF(s1 == 1) is1 = 2
22  IF(e1 == n1) ie1 = n1 - 1
23 
24  fac = 0.5_8 / dx
25  plane = (n1*n2)
26 
27  DO k = s3, e3
28  lik = (k-1)*plane
29  DO j = s2, e2
30  lij = (j-1)*n1 + lik
31  DO i = is1, ie1
32  li = lij + i
33  kn(li) = ( u(li+1) - u(li-1) ) * fac
34  END DO
35  END DO
36  END DO
37 
38 
39 END SUBROUTINE fd1interior
40 
41 SUBROUTINE fd1interiorx(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX, U, KN)
42 
43  IMPLICIT NONE
44  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
45  REAL(KIND=8), INTENT(IN) :: DX
46  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
47  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
48 
49  INTEGER(KIND=8) :: I, J, K, LI, LIK, LIJ, IS1, IE1
50  REAL(KIND=8) :: fac
51  INTEGER(KIND=8) :: plane
52 
53  is1 = s1
54  ie1 = e1
55  IF(s1 == 1) is1 = 2
56  IF(e1 == n1) ie1 = n1 - 1
57 
58  fac = 0.5_8 / dx
59  plane = (n1*n2)
60 
61  DO k = s3, e3
62  lik = (k-1)*plane
63  DO j = s2, e2
64  lij = (j-1)*n1 + lik
65  DO i = is1, ie1
66  li = lij + i
67  kn(li) = ( u(li+1) - u(li-1) ) * fac
68  END DO
69  END DO
70  END DO
71 
72 
73 END SUBROUTINE fd1interiorx
74 
75 SUBROUTINE fd1interiory(N1,N2,N3,S1,E1,S2,E2,S3,E3,DY,U,KN)
76 
77  IMPLICIT NONE
78 
79  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
80  REAL(KIND=8), INTENT(IN) :: DY
81  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
82  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
83 
84  INTEGER(KIND=8) :: I, J, K, LI, LIK, LIJ
85  INTEGER(KIND=8) :: IS1, IS2, IS3, IE1, IE2, IE3
86  REAL(KIND=8) :: fac
87  INTEGER(KIND=8) :: plane
88 
89  is1 = s1
90  ie1 = e1
91  is2 = s2
92  ie2 = e2
93  is3 = s3
94  ie3 = e3
95 
96  IF(s2 == 1) is2 = 2
97  IF(e2 == n2) ie2 = n2 - 1
98 
99  fac = 0.5_8 / dy
100  plane = (n1*n2)
101 
102  DO k = is3, ie3
103  lik = (k-1)*plane
104  DO j = is2, ie2
105  lij = (j-1)*n1 + lik
106  DO i = is1, ie1
107  li = lij + i
108  kn(li) = ( u(li+n1) - u(li-n1) ) * fac
109  END DO
110  END DO
111  END DO
112 
113 
114 END SUBROUTINE fd1interiory
115 
116 SUBROUTINE fd1interiorz(N1,N2,N3,S1,E1,S2,E2,S3,E3,DZ,U,KN)
117 
118  IMPLICIT NONE
119 
120  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
121  REAL(KIND=8), INTENT(IN) :: DZ
122  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
123  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
124 
125  INTEGER(KIND=8) :: I, J, K, LI, LIK, LIJ
126  INTEGER(KIND=8) :: IS1, IS2, IS3, IE1, IE2, IE3
127  REAL(KIND=8) :: fac
128  INTEGER(KIND=8) :: plane
129 
130  is1 = s1
131  ie1 = e1
132  is2 = s2
133  ie2 = e2
134  is3 = s3
135  ie3 = e3
136 
137  IF(s3 == 1) is3 = 2
138  IF(e3 == n3) ie3 = n3 - 1
139 
140  fac = 0.5_8 / dz
141  plane = (n1*n2)
142 
143  DO k = is3, ie3
144  lik = (k-1)*plane
145  DO j = is2, ie2
146  lij = (j-1)*n1 + lik
147  DO i = is1, ie1
148  li = lij + i
149  kn(li) = ( u(li+plane) - u(li-plane) ) * fac
150  END DO
151  END DO
152  END DO
153 
154 
155 END SUBROUTINE fd1interiorz
156 
157 SUBROUTINE fd1interior2(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX,DY,U,KX,KY)
158 
159  IMPLICIT NONE
160 
161  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
162  REAL(KIND=8), INTENT(IN) :: DX,DY
163  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
164  REAL(KIND=8), INTENT(OUT) :: KX(n1*n2*n3)
165  REAL(KIND=8), INTENT(OUT) :: KY(n1*n2*n3)
166 
167  INTEGER(KIND=8) :: I, J, K
168  INTEGER(KIND=8) :: LI, LIK, LIJ, IS1, IE1
169  INTEGER(KIND=8) :: LII, IS2, IE2
170  REAL(KIND=8) :: fac1,fac2
171  INTEGER(KIND=8) :: plane
172 
173  is1 = s1
174  ie1 = e1
175  is2 = s2
176  ie2 = e2
177 
178  IF(s1 == 1) is1 = 2
179  IF(e1 == n1) ie1 = n1 - 1
180  IF(s2 == 1) is2 = 2
181  IF(e2 == n2) ie2 = n2 - 1
182 
183 
184 
185  fac1 = 0.5_8 / dx
186  fac2 = 0.5_8 / dy
187 
188  plane = (n1*n2)
189 
190  DO k = s3, e3
191  lik = (k-1)*plane
192 
193  DO j = s2, e2
194  lij = (j-1)*n1 + lik
195  DO i = is1, ie1
196  li = lij + i
197  kx(li) = ( u(li+1) - u(li-1) ) * fac1
198  END DO
199  END DO
200 
201  DO i = s1, e1
202  lii = lik + i
203  DO j = is2, ie2
204  lij = (j-1)*n1 + lii
205  ky(lij) = ( u(lij+n1) - u(lij-n1) ) * fac2
206  END DO
207  END DO
208 
209  END DO
210 
211 END SUBROUTINE fd1interior2
212 
213 SUBROUTINE fd1interior3(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX,DY,DZ,U,KX,KY,KZ)
214 
215  IMPLICIT NONE
216 
217  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
218  REAL(KIND=8), INTENT(IN) :: DX,DY,DZ
219  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
220  REAL(KIND=8), INTENT(OUT) :: KX(n1*n2*n3)
221  REAL(KIND=8), INTENT(OUT) :: KY(n1*n2*n3)
222  REAL(KIND=8), INTENT(OUT) :: KZ(n1*n2*n3)
223 
224  INTEGER(KIND=8) :: I, J, K
225  INTEGER(KIND=8) :: LI, LIK, LIJ, IS1, IE1
226  INTEGER(KIND=8) :: LII, IS2, IE2, IS3, IE3
227  REAL(KIND=8) :: fac1,fac2,fac3
228  INTEGER(KIND=8) :: plane
229 
230  is1 = s1
231  ie1 = e1
232 
233  is2 = s2
234  ie2 = e2
235 
236  is3 = s3
237  ie3 = e3
238 
239  ! Operate on the interior only
240  IF(s1 == 1) is1 = 2
241  IF(e1 == n1) ie1 = n1 - 1
242 
243  IF(s2 == 1) is2 = 2
244  IF(e2 == n2) ie2 = n2 - 1
245 
246  IF(s3 == 1) is3 = 2
247  IF(e3 == n3) ie3 = n3 - 1
248 
249 
250  fac1 = 0.5_8 / dx
251  fac2 = 0.5_8 / dy
252  fac3 = 0.5_8 / dz
253 
254  plane = (n1*n2)
255 
256  DO k = s3, e3
257 
258  lik = (k-1)*plane
259 
260  DO j = s2, e2
261  lij = (j-1)*n1 + lik
262  DO i = is1, ie1
263  li = lij + i
264  kx(li) = ( u(li+1) - u(li-1) ) * fac1
265  END DO
266  END DO
267 
268  DO i = s1, e1
269  lii = lik + i
270  DO j = is2, ie2
271  lij = (j-1)*n1 + lii
272  ky(lij) = ( u(lij+n1) - u(lij-n1) ) * fac2
273  END DO
274  END DO
275 
276  END DO
277 
278  DO j = s2, e2
279  lij = (j-1)*n1
280  DO i = s1, e1
281  li = lij + i
282  DO k = is3, ie3
283  lik = (k-1)*plane + li
284  kz(lik) = ( u(lik+plane) - u(lik-plane) ) * fac3
285  END DO
286  END DO
287  END DO
288 
289 END SUBROUTINE fd1interior3
290 
291 SUBROUTINE fd1grad3(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX,DY,DZ,U,GRADU)
292 
293  IMPLICIT NONE
294 
295  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
296  REAL(KIND=8), INTENT(IN) :: DX,DY,DZ
297  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3)
298  REAL(KIND=8), INTENT(OUT) :: GRADU(3*n1*n2*n3)
299 ! REAL(KIND=8), INTENT(OUT) :: KY(N1*N2*N3)
300 ! REAL(KIND=8), INTENT(OUT) :: KZ(N1*N2*N3)
301 
302  INTEGER(KIND=8) :: I, J, K
303  INTEGER(KIND=8) :: LI, LIK, LIJ, IS1, IE1
304  INTEGER(KIND=8) :: LII, IS2, IE2, IS3, IE3
305  REAL(KIND=8) :: fac1,fac2,fac3
306  INTEGER(KIND=8) :: plane
307  INTEGER(KIND=8) :: numPoints,numPoints2
308 
309  numpoints = n1*n2*n3
310  numpoints2 = 2*n1*n2*n3
311 
312  is1 = s1
313  ie1 = e1
314 
315  is2 = s2
316  ie2 = e2
317 
318  is3 = s3
319  ie3 = e3
320 
321  ! Operate on the interior only
322  IF(s1 == 1) is1 = 2
323  IF(e1 == n1) ie1 = n1 - 1
324 
325  IF(s2 == 1) is2 = 2
326  IF(e2 == n2) ie2 = n2 - 1
327 
328  IF(s3 == 1) is3 = 2
329  IF(e3 == n3) ie3 = n3 - 1
330 
331 
332  fac1 = 0.5_8 / dx
333  fac2 = 0.5_8 / dy
334  fac3 = 0.5_8 / dz
335 
336  plane = (n1*n2)
337 
338  DO k = s3, e3
339 
340  lik = (k-1)*plane
341 
342  DO j = s2, e2
343  lij = (j-1)*n1 + lik
344  DO i = is1, ie1
345  li = lij + i
346  gradu(li) = ( u(li+1) - u(li-1) ) * fac1
347  END DO
348  END DO
349 
350  DO i = s1, e1
351  lii = lik + i
352  DO j = is2, ie2
353  lij = (j-1)*n1 + lii
354  gradu(lij+numpoints) = ( u(lij+n1) - u(lij-n1) ) * fac2
355  END DO
356  END DO
357 
358  END DO
359 
360  DO j = s2, e2
361  lij = (j-1)*n1
362  DO i = s1, e1
363  li = lij + i
364  DO k = is3, ie3
365  lik = (k-1)*plane + li
366  gradu(lik+numpoints2) = ( u(lik+plane) - u(lik-plane) ) * fac3
367  END DO
368  END DO
369  END DO
370 
371 END SUBROUTINE fd1grad3
372 
373 SUBROUTINE fd1div3(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX,DY,DZ,UVEC,DIVU)
374 
375  IMPLICIT NONE
376 
377  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
378  REAL(KIND=8), INTENT(IN) :: DX,DY,DZ
379  REAL(KIND=8), INTENT(IN) :: UVEC(3*n1*n2*n3)
380  REAL(KIND=8), INTENT(OUT) :: DIVU(n1*n2*n3)
381 ! REAL(KIND=8), INTENT(OUT) :: KY(N1*N2*N3)
382 ! REAL(KIND=8), INTENT(OUT) :: KZ(N1*N2*N3)
383 
384  INTEGER(KIND=8) :: I, J, K
385  INTEGER(KIND=8) :: LI, LIK, LIJ, IS1, IE1
386  INTEGER(KIND=8) :: LII, IS2, IE2, IS3, IE3
387  REAL(KIND=8) :: fac1,fac2,fac3
388  INTEGER(KIND=8) :: plane
389  INTEGER(KIND=8) :: numPoints,numPoints2
390 
391  numpoints = n1*n2*n3
392  numpoints2 = 2*n1*n2*n3
393 
394  is1 = s1
395  ie1 = e1
396 
397  is2 = s2
398  ie2 = e2
399 
400  is3 = s3
401  ie3 = e3
402 
403  ! Operate on the interior only
404  IF(s1 == 1) is1 = 2
405  IF(e1 == n1) ie1 = n1 - 1
406 
407  IF(s2 == 1) is2 = 2
408  IF(e2 == n2) ie2 = n2 - 1
409 
410  IF(s3 == 1) is3 = 2
411  IF(e3 == n3) ie3 = n3 - 1
412 
413 
414  fac1 = 0.5_8 / dx
415  fac2 = 0.5_8 / dy
416  fac3 = 0.5_8 / dz
417 
418  plane = (n1*n2)
419 
420  DO k = s3, e3
421 
422  lik = (k-1)*plane
423 
424  DO j = s2, e2
425  lij = (j-1)*n1 + lik
426  DO i = is1, ie1
427  li = lij + i
428  divu(li) = ( uvec(li+1) - uvec(li-1) ) * fac1
429  END DO
430  END DO
431 
432  DO i = s1, e1
433  lii = lik + i
434  DO j = is2, ie2
435  lij = (j-1)*n1 + lii
436  divu(lij) = divu(lij) + ( uvec((lij+n1)+numpoints) - uvec((lij-n1)+numpoints) ) * fac2
437  END DO
438  END DO
439 
440  END DO
441 
442  DO j = s2, e2
443  lij = (j-1)*n1
444  DO i = s1, e1
445  li = lij + i
446  DO k = is3, ie3
447  lik = (k-1)*plane + li
448  divu(lik) = divu(lik) + ( uvec((lik+plane)+numpoints2) - uvec((lik-plane)+numpoints2) ) * fac3
449  END DO
450  END DO
451  END DO
452 
453 END SUBROUTINE fd1div3
454 
455 ! @ICE block=FD1DIV3ALL
456 SUBROUTINE fd1div3all(N1,N2,N3,OPINTERVAL,DX,DY,DZ,&
457  UVHLEFT,UVHRIGHT,UVHFRONT,UVHBACK,UVHBOTTOM,UVHTOP, &
458  UVEC,DIVU)
459 
460  IMPLICIT NONE
461 
462  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3
463  INTEGER(KIND=8), INTENT(IN) :: OPINTERVAL(6) ! S1, E1, S2, E2, S3, E3
464  REAL(KIND=8), INTENT(IN) :: DX,DY,DZ
465  REAL(KIND=8), INTENT(IN) :: UVHLEFT(3*n2*n3)
466  REAL(KIND=8), INTENT(IN) :: UVHRIGHT(3*n2*n3)
467  REAL(KIND=8), INTENT(IN) :: UVHFRONT(3*n1*n3)
468  REAL(KIND=8), INTENT(IN) :: UVHBACK(3*n1*n3)
469  REAL(KIND=8), INTENT(IN) :: UVHBOTTOM(3*n1*n2)
470  REAL(KIND=8), INTENT(IN) :: UVHTOP(3*n1*n2)
471  REAL(KIND=8), INTENT(IN) :: UVEC(3*n1*n2*n3)
472  REAL(KIND=8), INTENT(OUT) :: DIVU(n1*n2*n3)
473  ! REAL(KIND=8), INTENT(OUT) :: KY(N1*N2*N3)
474  ! REAL(KIND=8), INTENT(OUT) :: KZ(N1*N2*N3)
475 
476  INTEGER(KIND=8) :: I, J, K
477  INTEGER(KIND=8) :: IH,IU,JH,JU,KH,KU
478  INTEGER(KIND=8) :: LI, LIK, LIJ, IS1, IE1
479  INTEGER(KIND=8) :: LII, IS2, IE2, IS3, IE3
480  INTEGER(KIND=8) :: S1,E1,S2,E2,S3,E3
481  REAL(KIND=8) :: fac1,fac2,fac3
482  INTEGER(KIND=8) :: plane,plane2,numHalo
483  INTEGER(KIND=8) :: numPoints,numPoints2
484 
485  divu(:) = 0.0_8
486 
487  numpoints = n1*n2*n3
488  numpoints2 = 2*n1*n2*n3
489 
490  fac1 = 0.5_8 / dx
491  fac2 = 0.5_8 / dy
492  fac3 = 0.5_8 / dz
493 
494  plane = (n1*n2)
495  plane2 = 2*n1*n2
496 
497  s1 = opinterval(1)
498  e1 = opinterval(2)
499 
500  s2 = opinterval(3)
501  e2 = opinterval(4)
502 
503  s3 = opinterval(5)
504  e3 = opinterval(6)
505 
506  is1 = s1
507  ie1 = e1
508 
509  is2 = s2
510  ie2 = e2
511 
512  is3 = s3
513  ie3 = e3
514 ! WRITE(*,*) 'DIV INTERVAL: ',IS1,IE1,IS2,IE2,IS3,IE3
515  ! ============== X DIRECTION =================
516 
517  ! PROCESS LEFT BOUNDARY
518  IF(is1 == 1) THEN ! select only threads with intervals which include the left boundary
519 
520  DO k = s3, e3
521  kh = (k-1)*n2
522  ku = kh*n1 + 1
523  DO j = s2, e2
524  jh = kh + j
525  ju = ku + (j-1)*n1
526  divu(ju) = ( uvec(ju+1) - uvhleft(jh) ) * fac1
527 ! WRITE(*,*) 'LEFT BOUNDARY: DIV(',JU,') = ',DIVU(JU)
528  END DO
529  END DO
530 
531  ENDIF
532 
533  ! PROCESS RIGHT BOUNDARY
534  IF(ie1 == n1) THEN ! select only threads touching right boundary
535 
536  DO k = s3, e3
537  kh = (k-1)*n2
538  ku = kh*n1 + n1
539  DO j = s2, e2
540  jh = kh + j
541  ju = ku + (j-1)*n1
542  divu(ju) = ( uvhright(jh) - uvec(ju-1) ) * fac1
543 ! WRITE(*,*) 'RIGHT BOUNDARY: DIV(',JU,') = ',DIVU(JU)
544  END DO
545  END DO
546 
547  ENDIF
548 
549  ! PROCESS INTERIOR
550  IF((e1 > s1) .OR. ((e1 < n1) .AND. (s1 > 1))) THEN ! Necessary thread selection!!! (awful bug resolved)
551 
552  IF(s1 == 1) is1 = 2
553  IF(e1 == n1) ie1 = n1 - 1
554 
555  DO k = s3, e3
556  ku = (k-1)*plane
557  DO j = s2, e2
558  ju = (j-1)*n1 + ku
559  DO i = is1, ie1
560  iu = ju + i
561  divu(iu) = ( uvec(iu+1) - uvec(iu-1) ) * fac1
562  END DO
563  END DO
564  END DO
565  ENDIF
566 
567  is1 = s1
568  ie1 = e1
569 
570 ! DO K = S3, E3
571 ! KU = (K-1)*plane
572 ! DO J = S2, E2
573 ! JU = (J-1)*N1 + KU
574 ! DO I = IS1, IE1
575 ! IU = JU + I
576 ! WRITE(*,*) 'XDIVU(',IU,') = ', DIVU(IU)
577 ! END DO
578 ! END DO
579 ! END DO
580 
581  ! ========== Y DIRECTION ==============
582 
583  numhalo = n3*n1
584  ! PROCESS FRONT BOUNDARY
585  IF(s2 == 1) THEN
586  DO k = s3, e3
587  kh = (k-1)*n1
588  ku = kh*n2
589  DO i = s1, e1
590  ih = kh + i
591  iu = ku + i
592  divu(iu) = divu(iu) + ( uvec(iu+n1+numpoints) - uvhfront(ih+numhalo) ) * fac2
593  END DO
594  END DO
595  ENDIF
596 ! DO K = S3, E3
597 ! KU = (K-1)*plane
598 ! DO J = S2, E2
599 ! JU = (J-1)*N1 + KU
600 ! DO I = IS1, IE1
601 ! IU = JU + I
602 ! WRITE(*,*) 'after front YDIVU(',IU,') = ', DIVU(IU)
603 ! END DO
604 ! END DO
605 ! END DO
606 
607  ! PROCESS BACK BOUNDARY
608  IF(e2 == n2) THEN
609  j = (n2-1)*n1
610  DO k = s3, e3
611  kh = (k-1)*n1
612  ku = kh*n2 + j
613  DO i = s1, e1
614  ih = kh + i
615  iu = ku + i
616  divu(iu) = divu(iu) + ( uvhback(ih+numhalo) - uvec((iu-n1)+numpoints) ) * fac2
617 ! WRITE(*,*) 'DIVU(',IU,') = uvhyBACK(',IH,',',UVHBACK(IH+numHalo),') - UVECY(',IU-N1,',',UVEC((IU-N1)+numPoints),')'
618  END DO
619  END DO
620  ENDIF
621 
622 ! DO K = S3, E3
623 ! KU = (K-1)*plane
624 ! DO J = S2, E2
625 ! JU = (J-1)*N1 + KU
626 ! DO I = IS1, IE1
627 ! IU = JU + I
628 ! WRITE(*,*) 'after back YDIVU(',IU,') = ', DIVU(IU)
629 ! END DO
630 ! END DO
631 ! END DO
632 
633  ! PROCESS INTERIOR in Y
634  IF((e2 > s2) .OR. ((e2 < n2) .AND. (s2 > 1))) THEN
635 
636  IF(s2 == 1) is2 = 2
637  IF(e2 == n2) ie2 = n2 - 1
638 
639  DO k = s3, e3
640  ku = (k-1)*plane
641  DO j = is2, ie2
642  ju = (j-1)*n1 + ku
643  DO i = s1, e1
644  iu = ju + i
645  divu(iu) = divu(iu) + ( uvec(iu+n1+numpoints) - uvec(iu-n1+numpoints) ) * fac2
646  END DO
647  END DO
648  END DO
649 
650  is2 = s2
651  ie2 = e2
652 
653  ENDIF
654 
655  ! ========== Z DIRECTION ==============
656 
657  ! PROCESS BOTTOM BOUNDARY
658  IF(s3 == 1) THEN
659  k = 1
660  DO j = s2, e2
661  ju = (j-1)*n1
662  DO i = s1, e1
663  iu = ju + i
664  divu(iu) = divu(iu) + ( uvec(iu+plane+numpoints2) - uvhbottom(iu+plane2) ) * fac3
665 ! WRITE(*,*) 'BOTTOM BOUNDARY: DIV(',IU,') = ',DIVU(IU)
666  END DO
667  END DO
668  ENDIF
669 
670  ! PROCESS TOP BOUNDARY
671  IF(e3 == n3) THEN
672  k = (n3-1)*plane
673  DO j = s2, e2
674  jh = (j-1)*n1
675  ju = jh + k
676  DO i = s1, e1
677  ih = jh + i
678  iu = ju + i
679  divu(iu) = divu(iu) + ( uvhtop(ih+plane2) - uvec(iu-plane+numpoints2) ) * fac3
680  ! WRITE(*,*) 'TOP BOUNDARY: DIV(',IU,') = ',DIVU(IU)
681  END DO
682  END DO
683  ENDIF
684 
685  ! PROCESS INTERIOR in Z
686  IF((e3 > s3) .OR. ((e3 < n3) .AND. (s3 > 1))) THEN
687 
688  IF(s3 == 1) is3 = 2
689  IF(e3 == n3) ie3 = n3 - 1
690 
691  DO k = is3, ie3
692  ku = (k-1)*plane
693  DO j = s2, e2
694  ju = (j-1)*n1 + ku
695  DO i = s1, e1
696  iu = ju + i
697  divu(iu) = divu(iu) + ( uvec(iu+plane+numpoints2) - uvec(iu-plane+numpoints2) ) * fac3
698  END DO
699  END DO
700  END DO
701 
702  ENDIF
703 
704 END SUBROUTINE fd1div3all
705 
706 
707 
708 ! Use non-zero tau to specify boundary data
709 SUBROUTINE fd1leftboundary(N1,N2,N3,S1,E1,S2,E2,S3,E3, DX, u0, Tau, U, KN)
710 
711  IMPLICIT NONE
712 
713  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
714  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), DX, Tau, u0
715  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
716 
717  INTEGER(KIND=8) :: I, J, K, LI, plane, KI
718  REAL(KIND=8) :: fac,bcfac,bcterm
719 
720  fac = 1.0_8 / dx
721  plane = (n1*n2)
722  bcfac = tau*fac
723  bcterm = bcfac*u0
724 
725  DO k = s3, e3
726  ki = (k-1)*plane + 1
727  DO j = s2, e2
728  li = ki + (j-1)*n1
729  kn(li) = ( u(li+1) - u(li) ) * fac + bcfac*u(li) - bcterm
730  END DO
731  END DO
732 
733 END SUBROUTINE fd1leftboundary
734 
735 
736 ! Use non-zero tau to specify boundary data
737 SUBROUTINE fd1rightboundary(N1,N2,N3,S1,E1,S2,E2,S3,E3,DX, u0, Tau, U, KN)
738 
739  IMPLICIT NONE
740 
741  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
742  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), DX, Tau, u0
743  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
744  INTEGER(KIND=8) :: J, K, LI, plane, KI
745  REAL(KIND=8) :: fac,bcfac,bcterm
746 
747 
748  fac = 1.0_8 / dx
749  plane = (n1*n2)
750  bcfac = tau*fac
751  bcterm = bcfac*u0
752 
753  DO k = s3, e3
754  ki = (k-1)*plane
755  DO j = s2, e2
756  li = ki + j*n1
757  kn(li) = ( u(li) - u(li-1) ) * fac + bcfac*u(li) - bcterm
758  END DO
759  END DO
760 
761 END SUBROUTINE fd1rightboundary
762 
763 
764 SUBROUTINE fd1haloleft(N1,N2,N3,S1,E1,S2,E2,S3,E3, DX, U, KN, HL)
765 
766 
767  IMPLICIT NONE
768 
769  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
770  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), HL(n2*n3), DX
771  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
772 
773  INTEGER(KIND=8) :: J, K, LI, plane, LII, KI, KII
774  REAL(KIND=8) :: FAC
775 
776  fac = 0.5_8/dx
777  plane = (n1*n2)
778 
779  DO k = s3, e3
780  kii = (k-1)*n2
781  ki = (k-1)*plane + 1
782  DO j = s2, e2
783  lii = kii + j
784  li = ki + (j-1)*n1
785  kn(li) = ( u(li+1) - hl(lii) ) * fac
786  END DO
787  END DO
788 
789 END SUBROUTINE fd1haloleft
790 
791 
792 SUBROUTINE fd1haloright(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, U, KN, HR)
793 
794  IMPLICIT NONE
795 
796  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
797  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), HR(n2*n3), DX
798  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
799 
800  INTEGER(KIND=8) :: J, K, LI, plane, LII, KI, KII
801  REAL(KIND=8) :: fac
802 
803 
804  fac = 0.5_8/dx
805  plane = (n1*n2)
806 
807  DO k = s3, e3
808  kii = (k-1)*n2
809  ki = (k-1)*plane
810  DO j = s2, e2
811  lii = kii + j
812  li = ki + j*n1
813  kn(li) = ( hr(lii) - u(li-1) ) * fac
814  END DO
815  END DO
816 
817 END SUBROUTINE fd1haloright
818 
819 SUBROUTINE fd1halofront(N1,N2,N3,S1,E1,S2,E2,S3,E3, DY, U, KN, HL)
820 
821 
822  IMPLICIT NONE
823 
824  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
825  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), HL(n1*n3), DY
826  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
827 
828  INTEGER(KIND=8) :: I, J, K, LI, plane, LII, KI, KII
829  REAL(KIND=8) :: FAC
830 
831  fac = 0.5_8/dy
832  plane = n1*n2
833 
834  DO k = s3, e3
835  kii = (k-1)*n1
836  ki = (k-1)*plane
837  DO i = s1, e1
838  lii = kii + i
839  li = ki + i
840  kn(li) = ( u(li+n1) - hl(lii) ) * fac
841  END DO
842  END DO
843 
844 END SUBROUTINE fd1halofront
845 
846 SUBROUTINE fd1haloback(N1,N2,N3,S1,E1,S2,E2,S3,E3, DY, U, KN, HL)
847 
848 
849  IMPLICIT NONE
850 
851  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
852  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), HL(n1*n3), DY
853  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
854 
855  INTEGER(KIND=8) :: I, J, K, LI, plane, LII, KI, KII
856  REAL(KIND=8) :: FAC
857 
858  fac = 0.5_8/dy
859  plane = n1*n2
860  j = (n2 - 1)*n1
861 
862  DO k = s3, e3
863  kii = (k-1)*n1
864  ki = (k-1)*plane + j
865  DO i = s1, e1
866  lii = kii + i
867  li = ki + i
868  kn(li) = ( hl(lii) - u(li-n1)) * fac
869  END DO
870  END DO
871 
872 END SUBROUTINE fd1haloback
873 
874 SUBROUTINE fd1halotop(N1,N2,N3,S1,E1,S2,E2,S3,E3, DZ, U, KN, HL)
875 
876 
877  IMPLICIT NONE
878 
879  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
880  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), HL(n1*n2), DZ
881  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
882 
883  INTEGER(KIND=8) :: I, J, K, LI, plane, LII, KI, KII
884  REAL(KIND=8) :: FAC
885 
886  fac = 0.5_8/dz
887  plane = n1*n2
888  k = (n3-1)*plane
889 
890  DO j = s2, e2
891  kii = (j-1)*n1
892  ki = k + kii
893  DO i = s1, e1
894  lii = kii + i
895  li = ki + i
896  kn(li) = ( hl(lii) - u(li-plane) ) * fac
897  END DO
898  END DO
899 
900 END SUBROUTINE fd1halotop
901 
902 SUBROUTINE fd1halobottom(N1,N2,N3,S1,E1,S2,E2,S3,E3, DZ, U, KN, HL)
903 
904 
905  IMPLICIT NONE
906 
907  INTEGER(KIND=8), INTENT(IN) :: N1, N2, N3, S1, E1, S2, E2, S3, E3
908  REAL(KIND=8), INTENT(IN) :: U(n1*n2*n3), HL(n1*n2), DZ
909  REAL(KIND=8), INTENT(OUT) :: KN(n1*n2*n3)
910 
911  INTEGER(KIND=8) :: I, J, K, LI, plane, LII, KI, KII
912  REAL(KIND=8) :: FAC
913 
914  fac = 0.5_8/dz
915  plane = n1*n2
916 
917 
918  DO j = s2, e2
919  kii = (j-1)*n1
920  DO i = s1, e1
921  lii = kii + i
922  kn(lii) = ( u(li+plane) - hl(lii) ) * fac
923  END DO
924  END DO
925 
926 END SUBROUTINE fd1halobottom
927 
fd1halofront
subroutine fd1halofront(N1, N2, N3, S1, E1, S2, E2, S3, E3, DY, U, KN, HL)
Definition: FD1.f90:820
fd1div3all
subroutine fd1div3all(N1, N2, N3, OPINTERVAL, DX, DY, DZ, UVHLEFT, UVHRIGHT, UVHFRONT, UVHBACK, UVHBOTTOM, UVHTOP, UVEC, DIVU)
Definition: FD1.f90:459
fd1interior3
subroutine fd1interior3(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, DY, DZ, U, KX, KY, KZ)
Definition: FD1.f90:214
fd1div3
subroutine fd1div3(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, DY, DZ, UVEC, DIVU)
Definition: FD1.f90:374
fd1grad3
subroutine fd1grad3(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, DY, DZ, U, GRADU)
Definition: FD1.f90:292
fd1interiorz
subroutine fd1interiorz(N1, N2, N3, S1, E1, S2, E2, S3, E3, DZ, U, KN)
Definition: FD1.f90:117
fd1interiory
subroutine fd1interiory(N1, N2, N3, S1, E1, S2, E2, S3, E3, DY, U, KN)
Definition: FD1.f90:76
fd1haloback
subroutine fd1haloback(N1, N2, N3, S1, E1, S2, E2, S3, E3, DY, U, KN, HL)
Definition: FD1.f90:847
fd1halobottom
subroutine fd1halobottom(N1, N2, N3, S1, E1, S2, E2, S3, E3, DZ, U, KN, HL)
Definition: FD1.f90:903
fd1load
subroutine fd1load
Definition: FD1.f90:2
fd1leftboundary
subroutine fd1leftboundary(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, u0, Tau, U, KN)
Definition: FD1.f90:710
fd1haloright
subroutine fd1haloright(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, U, KN, HR)
Definition: FD1.f90:793
fd1rightboundary
subroutine fd1rightboundary(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, u0, Tau, U, KN)
Definition: FD1.f90:738
fd1haloleft
subroutine fd1haloleft(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, U, KN, HL)
Definition: FD1.f90:765
fd1halotop
subroutine fd1halotop(N1, N2, N3, S1, E1, S2, E2, S3, E3, DZ, U, KN, HL)
Definition: FD1.f90:875
fd1interior2
subroutine fd1interior2(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, DY, U, KX, KY)
Definition: FD1.f90:158
fd1interiorx
subroutine fd1interiorx(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, U, KN)
Definition: FD1.f90:42
fd1interior
subroutine fd1interior(N1, N2, N3, S1, E1, S2, E2, S3, E3, DX, U, KN)
Definition: FD1.f90:7