femwind/fortran/ndt_w_test2

   1 program w_assembly_test
   2
   3 use module_w_assembly
   4 use module_hexa
   5 use module_io_matlab
   6
   7 !Purpose: Create Arrays of Wind Vector Component Values at Center Points of Spatial Grid
   8 !In:
   9 !A Coefficient Matrix size 3X3, symmetric positive definite
  10 !u0, v0, w0  Initial wind speed values in x,y,z direction at grid cell centers
  11 !X,Y,Z       3-D Physical Location in the Mesh Grid
  12 !iflags1 iflags1 = 1 returns Kloc and Jg from hexa, iflags2 = 2 returns Floc and Jg from hexa
  13 !iflags2 iflags2 =1  indicates add initial wind to calculated wind
  14 !out:
  15 !U,V,W Computed wind values in x,y,z direction
  16
  17 implicit none
  18
  19 real, pointer:: Amat(:,:), u0mat(:,:,:),v0mat(:,:,:), w0mat(:,:,:), Umat(:,:,:),              &
  20                 Vmat(:,:,:), Wmat(:,:,:), u0(:,:,:), v0(:,:,:), w0(:,:,:),         &
  21                 U(:,:,:), V(:,:,:),W(:,:,:),lambda(:,:,:), lambdamat(:,:,:),       &  ! Calculated final windFinal
  22                 Xmat(:,:,:),Ymat(:,:,:),Zmat(:,:,:), X(:,:,:),Y(:,:,:),Z(:,:,:)
  23
  24
  25 real, pointer :: a1(:), a2(:)
  26 integer :: n1(2),lambda_dim(3),u_dim(3), x_dim(3)
  27
  28 integer :: &
  29     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  30     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  31     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  32     ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
  33     iuts, iute, kuts, kute, juts, jute,                       &
  34     iums, iume, kums, kume, jums, jume
  35
  36
  37
  38 integer :: i,j,k,jx
  39 integer :: aflags(2) = (/3,1/)                 !Set iflags=1 to construct K in hexa module, iflags = 3 to construct Jg
  40 !integer :: iflags2 = 1
  41 integer :: usize(3)
  42
  43 call read_array_nd(a1,n1,'A') !Recovering X-Matrix and dimension of X matrix
  44 if (n1(1).ne.3.or.n1(2).ne.3)then
  45     call crash('A must be 3 by 3')
  46     stop
  47 endif
  48
  49 Amat = reshape(a1,n1)
  50
  51 !call read_array_nd(a,n,'u')
  52 !allocate(u_m(n(1),n(2),n(3)))
  53 !u_m = reshape(a,n)
  54
  55 call read_array(lambda, 'lambda')
  56 lambda_dim = shape(lambda)
  57
  58
  59 ! read input arrays in ikj index ordering and tight bounds
  60 call read_array(X,'X')
  61 call read_array(Y,'Y')
  62 call read_array(Z,'Z')
  63
  64 call read_array(u0,'u0')        !Recovering Inital Wind Arrays
  65 call read_array(v0,'v0')
  66 call read_array(w0,'w0')
  67
  68 x_dim = shape(X)
  69 u_dim = shape(u0)
  70
  71 !Checking that dimensions of lambda and are X arrays are consistent
  72 if (x_dim(1).ne.lambda_dim(1).or.x_dim(2).ne.lambda_dim(2).or.x_dim(3).ne.lambda_dim(3)) then
  73     call crash('Lambda dimensions must equal the dimensions of X')
  74     stop
  75 endif
  76
  77 ifts = 1
  78 ifte = x_dim(1)
  79 jfts = 1
  80 jfte = x_dim(3)
  81 kfts = 1
  82 kfte = x_dim(2)
  83 ifms = ifts - 1
  84 ifme = ifte + 1
  85 jfms = jfts - 1
  86 jfme = jfte + 1
  87 kfms = kfts - 1
  88 kfme = kfte + 1
  89
  90 iuts = 1
  91 iute = u_dim(1)
  92 juts = 1
  93 jute = u_dim(3)
  94 kuts = 1
  95 kute = u_dim(2)
  96 iums = iuts - 1
  97 iume = iute + 1
  98 jums = juts - 1
  99 jume = jute + 1
 100 kums = kuts - 1
 101 kume = kute + 1
 102
 103
 104
 105 allocate(lambdamat(ifms:ifme, kfms:kfme, jfms:jfme))
 106
 107
 108 allocate(Xmat(ifms:ifme,kfms:kfme,jfms:jfme))
 109 allocate(Ymat(ifms:ifme,kfms:kfme,jfms:jfme))
 110 allocate(Zmat(ifms:ifme,kfms:kfme,jfms:jfme))
 111
 112 allocate(u0mat(iums:iume,kums:kume,jums:jume))
 113 allocate(v0mat(iums:iume,kums:kume,jums:jume))
 114 allocate(w0mat(iums:iume,kums:kume,jums:jume))
 115
 116
 117
 118
 119 do j=jfts,jfte
 120   do k=kfts,kfte
 121     do i=ifts,ifte
 122         Xmat(i,k,j) = X(i,k,j)
 123         Ymat(i,k,j) = Y(i,k,j)
 124         Zmat(i,k,j) = Z(i,k,j)
 125         lambdamat(i,k,j) = lambda(i,k,j)
 126     enddo
 127   enddo
 128 enddo
 129
 130 ! copy the input data to tile sized bounds
 131 do j=juts,jute
 132   do k=kuts,kute
 133     do i=iuts,iute
 134     u0mat(i,k,j) = u0(i,k,j)
 135         v0mat(i,k,j) = v0(i,k,j)
 136         w0mat(i,k,j) = w0(i,k,j)
 137     enddo
 138   enddo
 139 enddo
 140
 141
 142
 143
 144 !write(*,'(a)')'calling w_assembly'
 145 call w_assembly(  &
 146   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 147   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 148   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 149   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 150   lambda,u0, v0, w0,                            & !Input from femwind, u0, v0, w0
 151   Amat, X, Y, Z,                                      & !Spatial Grid Data
 152   U,V,W)
 153
 154 !write(*,'(a,3i8)')'copying the output data to array size ',n2,msize
 155
 156 allocate(Umat(1:u_dim(1),1:u_dim(3),1:u_dim(2)))
 157 allocate(Vmat(1:u_dim(1),1:u_dim(3),1:u_dim(2)))
 158 allocate(Wmat(1:u_dim(1),1:u_dim(3),1:u_dim(2)))
 159 !keep track of indexing
 160 do j=juts,jute
 161   do k=kuts,kute
 162     do i=iuts,iute
 163             Umat(i,k,j)=U(i,k,j)
 164             Vmat(i,k,j)=V(i,k,j)
 165             Wmat(i,k,j)=W(i,k,j)
 166     enddo
 167   enddo
 168 enddo
 169
 170
 171 usize = (/ifte-ifts+1,kfte-kfts+1,jfte-jfts+1/)
 172 call write_array_nd(reshape(Umat,(/product(usize)/)),usize,'U')
 173 call write_array_nd(reshape(Vmat,(/product(usize)/)),usize,'V')
 174 call write_array_nd(reshape(Wmat,(/product(usize)/)),usize,'W')
 175
 176
 177 end program w_assembly_test