femwind/fortran/w_assembly_test.f90

   1 program w_assembly_test
   2
   3 use module_w_assembly
   4 use module_hexa
   5 use module_utils
   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:: 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(:,:,:),A_m(:,:,:)
  23
  24
  25 real, pointer :: a1(:), a2(:)
  26 integer :: n1(2),lambda_dim(3),u_dim(3), x_dim(3)
  27 real :: Amat(3,3)
  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
  34
  35
  36 integer :: i,j,k,jx
  37 integer :: usize(3)
  38
  39 call read_array(A_m,'A_test')
  40 Amat = reshape(A_m,(/3,3/))
  41
  42 call read_array(lambda, 'lambda_test')
  43 lambda_dim = shape(lambda)
  44
  45 print *, 'lambda array has shape', shape(lambda)
  46
  47 ! read input arrays in ikj index ordering and tight bounds
  48 call read_array(X,'X_test')
  49 call read_array(Y,'Y_test')
  50 call read_array(Z,'Z_test')
  51
  52 call read_array(u0,'u0_test')        !Recovering Inital Wind Arrays
  53 call read_array(v0,'v0_test')
  54 call read_array(w0,'w0_test')
  55
  56 x_dim = shape(X)
  57 u_dim = shape(u0)
  58
  59 !Checking that dimensions of lambda and are X arrays are consistent
  60 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
  61     call crash('Lambda dimensions must equal the dimensions of X')
  62     stop
  63 endif
  64
  65 ifts = 1
  66 ifte = x_dim(1)-1
  67 jfts = 1
  68 jfte = x_dim(3)-1
  69 kfts = 1
  70 kfte = x_dim(2)-1
  71 ifms = ifts-1
  72 ifme = ifte+2
  73 jfms = jfts-1
  74 jfme = jfte+2
  75 kfms = kfts-1
  76 kfme = kfte+2
  77 ifds = ifts
  78 ifde = ifte
  79 jfds = jfts
  80 jfde = jfte
  81 kfds = kfts
  82 kfde = kfte
  83
  84 allocate(lambdamat(ifms:ifme,kfms:kfme,jfms:jfme))
  85
  86
  87
  88 allocate(Xmat(ifms:ifme,kfms:kfme,jfms:jfme))
  89 allocate(Ymat(ifms:ifme,kfms:kfme,jfms:jfme))
  90 allocate(Zmat(ifms:ifme,kfms:kfme,jfms:jfme))
  91
  92 allocate(u0mat(ifms:ifme,kfms:kfme,jfms:jfme))
  93 allocate(v0mat(ifms:ifme,kfms:kfme,jfms:jfme))
  94 allocate(w0mat(ifms:ifme,kfms:kfme,jfms:jfme))
  95
  96 ! copy the input data to tile sized bounds
  97 ! X Y Z are corner based, upper bound larger by one
  98 do j=jfts,jfte+1
  99   do k=kfts,kfte+1
 100     do i=ifts,ifte+1
 101         Xmat(i,k,j) = X(i,k,j)
 102         Ymat(i,k,j) = Y(i,k,j)
 103         Zmat(i,k,j) = Z(i,k,j)
 104         lambdamat(i,k,j) = lambda(i,k,j)
 105     enddo
 106   enddo
 107 enddo
 108
 109 do j=jfts,jfte
 110   do k=kfts,kfte
 111     do i=ifts,ifte
 112         u0mat(i,k,j) = u0(i,k,j)
 113         v0mat(i,k,j) = v0(i,k,j)
 114         w0mat(i,k,j) = w0(i,k,j)
 115     enddo
 116   enddo
 117 enddo
 118
 119 print *, 'u0 is', u0
 120
 121
 122 allocate(U(ifms:ifme,kfms:kfme,jfms:jfme))
 123 allocate(V(ifms:ifme,kfms:kfme,jfms:jfme))
 124 allocate(W(ifms:ifme,kfms:kfme,jfms:jfme))
 125
 126 U = 0.
 127 V = 0.
 128 W = 0.
 129
 130 !write(*,'(a)')'calling w_assembly'
 131 call w_assembly(  &
 132   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 133   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 134   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 135   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 136   lambdamat,u0mat, v0mat, w0mat, Amat,                      & !Input from femwind, u0, v0, w0,Spatial Grid Data
 137   Xmat, Ymat, Zmat,                                         &
 138   U,V,W)
 139
 140 !write(*,'(a,3i8)')'copying the output data to array size ',n2,msize
 141
 142 allocate(Umat(ifts:ifte,kfts:kfte,jfts:jfte))
 143 allocate(Vmat(ifts:ifte,kfts:kfte,jfts:jfte))
 144 allocate(Wmat(ifts:ifte,kfts:kfte,jfts:jfte))
 145 !keep track of indexing
 146 do j=jfts,jfte
 147   do k=kfts,kfte
 148     do i=ifts,ifte
 149             Umat(i,k,j)=U(i,k,j)
 150             Vmat(i,k,j)=V(i,k,j)
 151             Wmat(i,k,j)=W(i,k,j)
 152     enddo
 153   enddo
 154 enddo
 155 !print *, 'Shape of Umat is', shape(Umat)
 156
 157 call write_array(Umat,'U_test')
 158 call write_array(Vmat,'V_test')
 159 call write_array(Wmat,'W_test')
 160
 161
 162 end program w_assembly_test