1subroutine formFunction_C(nx, ny, nz, h, t, x, xdot, f) & 2 bind(C, name="formFunction") 3 use ISO_C_BINDING, only: C_INT, C_DOUBLE 4 implicit none 5 integer(kind=C_INT), intent(in) :: nx, ny, nz 6 real(kind=C_DOUBLE), intent(in) :: h(3), t 7 real(kind=C_DOUBLE), intent(in) :: x(nx,ny,nz), xdot(nx,ny,nz) 8 real(kind=C_DOUBLE), intent(inout) :: f(nx,ny,nz) 9 call formfunction_f(nx, ny, nz, h, t, x, xdot, f) 10end subroutine formFunction_C 11 12subroutine formInitial_C(nx, ny, nz, h, t, x) & 13 bind(C, name="formInitial") 14 use ISO_C_BINDING, only: C_INT, C_DOUBLE 15 implicit none 16 integer(kind=C_INT), intent(in) :: nx, ny, nz 17 real(kind=C_DOUBLE), intent(in) :: h(3), t 18 real(kind=C_DOUBLE), intent(inout) :: x(nx,ny,nz) 19 call forminitial_f(nx, ny, nz, h, t, x) 20end subroutine formInitial_C 21 22subroutine evalK (P, K) 23 real(kind=8), intent(in) :: P 24 real(kind=8), intent(out) :: K 25 if (P >= 0.0) then 26 K = 1.0 27 else 28 K = 1.0/(1+P**2) 29 end if 30end subroutine evalK 31 32subroutine fillK (P, K) 33 real(kind=8), intent(in) :: P(-1:1) 34 real(kind=8), intent(out) :: K(-1:1) 35 real(kind=8) :: Ka, Kb 36 call evalK((P(-1)+P( 0))/2.0, Ka) 37 call evalK((P( 0)+P( 1))/2.0, Kb) 38 K(-1) = -Ka 39 K( 0) = Ka + Kb 40 K(+1) = -Kb 41end subroutine fillK 42 43subroutine forminitial_f(nx, ny, nz, h, t, x) 44 implicit none 45 integer, intent(in) :: nx, ny, nz 46 real(kind=8), intent(in) :: h(3), t 47 real(kind=8), intent(inout) :: x(nx,ny,nz) 48 ! 49 x(:,:,:) = 0.0 50end subroutine forminitial_f 51 52subroutine formfunction_f(nx, ny, nz, h, t, x, xdot, f) 53 implicit none 54 integer, intent(in) :: nx, ny, nz 55 real(kind=8), intent(in) :: h(3), t 56 real(kind=8), intent(in) :: x(nx,ny,nz), xdot(nx,ny,nz) 57 real(kind=8), intent(inout) :: f(nx,ny,nz) 58 ! 59 integer :: i,j,k,ii(-1:1),jj(-1:1),kk(-1:1) 60 real(kind=8) :: K1(-1:1), K2(-1:1), K3(-1:1) 61 real(kind=8) :: P1(-1:1), P2(-1:1), P3(-1:1) 62 ! 63 do k=1,nz 64 do j=1,ny 65 do i=1,nx 66 ! 67 ii = [i-1, i, i+1] 68 if (i == 1) then 69 ii(-1) = i 70 else if (i == nx) then 71 ii(+1) = i 72 end if 73 ! 74 jj = [j-1, j, j+1] 75 if (j == 1) then 76 jj(-1) = j 77 else if (j == ny) then 78 jj(+1) = j 79 end if 80 ! 81 kk = [k-1, k, k+1] 82 if (k == 1) then 83 kk(-1) = k 84 else if (k == nz) then 85 kk(+1) = k 86 end if 87 ! 88 P1 = x(ii,j,k) 89 P2 = x(i,jj,k) 90 P3 = x(i,j,kk) 91 call fillK(P1,K1) 92 call fillK(P2,K2) 93 call fillK(P3,K3) 94 f(i,j,k) = & 95 xdot(i,j,k) + & 96 sum(K1*P1)/h(1)**2 + & 97 sum(K2*P2)/h(2)**2 + & 98 sum(K3*P3)/h(3)**2 99 ! 100 end do !i 101 end do !j 102 end do !k 103 ! 104 i = nx/4+1 105 j = ny/4+1 106 k = nz/2+1 107 f(i,j,k:nz) = f(i,j,k:nz) + 300.0 108 ! 109end subroutine formfunction_f 110