femwind/fortran/femwind_test.f90

   1 program femwind_test
   2
   3 use module_femwind
   4 use module_utils
   5 use module_common
   6
   7 implicit none
   8
   9 real, pointer:: X_m(:,:,:),Y_m(:,:,:),Z_m(:,:,:), &
  10                 u0_m(:,:,:), v0_m(:,:,:), w0_m(:,:,:),   &
  11                 u_m(:,:,:), v_m(:,:,:), w_m(:,:,:)
  12
  13 integer ::                          &
  14     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  15     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  16     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  17     ifts, ifte, kfts, kfte, jfts, jfte                          ! fire tile bounds
  18
  19 ! A, msize included from module femwind
  20 real, pointer:: X(:,:,:),Y(:,:,:),Z(:,:,:), &
  21                 u0(:,:,:), v0(:,:,:), w0(:,:,:),   &
  22                 u(:,:,:), v(:,:,:), w(:,:,:),Kmat(:,:,:,:),A_m(:,:,:)
  23
  24 integer:: i, j, k, n(3)
  25 real:: rate
  26
  27 call read_array(A_m,'A_input')  ! matrices read from Matlab are _m
  28 call read_array(X_m,'X_input')
  29 call read_array(Y_m,'Y_input')
  30 call read_array(Z_m,'Z_input')
  31 call read_array(u0_m, 'u0_input')
  32 call read_array(v0_m, 'v0_input')
  33 call read_array(w0_m, 'w0_input')
  34
  35 params%A = reshape(A_m,(/3,3/))
  36
  37 n = shape(X_m)
  38 mg(1)%nx = n(1)
  39 mg(1)%ny = n(3)
  40 mg(1)%nz = n(2)
  41
  42 call get_mg_dims(mg(1), &
  43     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
  44     ifms, ifme, kfms,kfme, jfms, jfme,            &
  45     ifps, ifpe, kfps,kfpe, jfps, jfpe,           & ! fire patch bounds
  46     ifts, ifte, kfts,kfte, jfts,jfte)
  47
  48 allocate(mg(1)%X(ifms:ifme,kfms:kfme,jfms:jfme))
  49 allocate(mg(1)%Y(ifms:ifme,kfms:kfme,jfms:jfme))
  50 allocate(mg(1)%Z(ifms:ifme,kfms:kfme,jfms:jfme))
  51 allocate(u0(ifms:ifme,kfms:kfme,jfms:jfme)) ! vector components called u v W
  52 allocate(v0(ifms:ifme,kfms:kfme,jfms:jfme))
  53 allocate(w0(ifms:ifme,kfms:kfme,jfms:jfme))
  54 allocate(u(ifms:ifme,kfms:kfme,jfms:jfme)) ! vector components called u v W
  55 allocate(v(ifms:ifme,kfms:kfme,jfms:jfme))
  56 allocate(w(ifms:ifme,kfms:kfme,jfms:jfme))
  57 allocate(Kmat(ifms:ifme,kfms:kfme,jfms:jfme,msize)) ! stifness matrix
  58
  59 ! copy the input data to tile sized bounds
  60 ! X Y Z are corner based, upper bound larger by one
  61 do j=jfts,jfte+1
  62   do k=kfts,kfte+1
  63     do i=ifts,ifte+1
  64         mg(1)%X(i,k,j) = X_m(i,k,j)
  65         mg(1)%Y(i,k,j) = Y_m(i,k,j)
  66         mg(1)%Z(i,k,j) = Z_m(i,k,j)
  67     enddo
  68   enddo
  69 enddo
  70
  71 mg(1)%dx = mg(1)%X(2,1,1)-mg(1)%X(1,1,1)
  72 mg(1)%dy = mg(1)%Y(1,1,2)-mg(1)%Y(1,1,1)
  73 allocate(mg(1)%dz(mg(1)%nz-1))
  74 do k=kfds,kfte
  75     mg(1)%dz(k)=mg(1)%Z(1,k+1,1)-mg(1)%Z(1,k,1)
  76 enddo
  77
  78 write(*,'(a)')'calling femwind_setup'
  79 call femwind_setup(mg)
  80 write(*,'(a)')'femwind_setup returned OK'
  81
  82 ! u is midpoint based
  83 do j=jfts,jfte
  84   do k=kfts,kfte
  85     do i=ifts,ifte
  86         u0(i,k,j) = u0_m(i,k,j)
  87         v0(i,k,j) = v0_m(i,k,j)
  88         w0(i,k,j) = w0_m(i,k,j)
  89     enddo
  90   enddo
  91 enddo
  92
  93
  94 write(*,'(a)')'calling femwind_solve'
  95 call femwind_solve(  mg,&
  96   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  97   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  98   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  99   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 100   u0, v0, w0,                                  & ! input arrays
 101   u, v, w,                                                  & ! output arrays
 102   rate)
 103 write(*,'(a)')'femwind_solve returned OK'
 104
 105 ! write output as is in 3D but with tight dimensions
 106 call write_array(u(ifts:ifte,kfts:kfte,jfts:jfte),'u')
 107 call write_array(v(ifts:ifte,kfts:kfte,jfts:jfte),'v')
 108 call write_array(w(ifts:ifte,kfts:kfte,jfts:jfte),'w')
 109 call write_scalar(rate,'rate')
 110
 111 end program femwind_test