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