femwind/fortran/ndt_boundary_conditions_test.f90

   1 program ndt_boundary_conditions_test
   2
   3 use module_boundary_conditions   ! testing only
   4 use module_io_matlab ! to read and write matrices as text files from matlab
   5
   6 implicit none
   7
   8 real, pointer:: kmat(:,:,:,:), &  ! fortran is not case sensitive
   9                 kmat_m(:,:,:,:),a(:)
  10 integer :: s(4),n(3)
  11
  12 integer :: msize, &
  13     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  14     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  15     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  16     ifts, ifte, kfts, kfte, jfts,jfte                            ! fire tile bounds
  17 integer :: i,j,k,jx
  18
  19 ! read input arrays in ikj index ordering and tight bounds
  20 call read_array_nd(a,s,'kmat')
  21 allocate(kmat_m(s(1),s(2),s(3),s(4)))
  22 kmat_m = reshape(a,s)
  23 n = s(1:3)
  24
  25 ifts = 1
  26 ifte = n(1)
  27 jfts = 1
  28 jfte = n(2)
  29 kfts = 1
  30 kfte = n(3)
  31 msize = s(4)
  32 if(msize.ne.14)call crash('msize must be 14')
  33 ifms = ifts-1
  34 ifme = ifte+1
  35 jfms = jfts-1
  36 jfme = jfte+1
  37 kfms = kfts-1
  38 kfme = kfte+1
  39
  40 ! allocate a little bigger with zeros in extra areas
  41 allocate(kmat(ifms:ifme,kfms:kfme,jfms:jfme,1:msize))
  42 kmat = 0.
  43
  44 ! copy the input data
  45 do j=jfts,jfte
  46   do k=kfts,kfte
  47     do i=ifts,ifte
  48       do jx = 1,msize
  49         kmat(i,k,j,jx) = kmat_m(i,j,k,jx)
  50       enddo
  51     enddo
  52   enddo
  53 enddo
  54
  55 write(*,'(a)')'calling ntd_boundary_conditions'
  56 call ndt_boundary_conditions(  &
  57   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  58   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  59   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  60   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
  61   kmat)
  62
  63 ! copy the output data
  64 do j=jfts,jfte
  65   do k=kfts,kfte
  66     do i=ifts,ifte
  67       do jx = 1,msize
  68         kmat_m(i,k,j,jx) = kmat(i,j,k,jx)
  69       enddo
  70     enddo
  71   enddo
  72 enddo
  73
  74 call write_array_nd(reshape(kmat_m,(/product(s)/)),s,'kmat')
  75
  76 end program ndt_boundary_conditions_test