separated constant part in module_w_assembly.f90
[wrf-fire-matlab.git] / femwind / fortran / module_w_assembly.f90
blobde47fa5534306a33ac61c7b9d59f46aa217ee92a
1 module module_w_assembly
2 use module_hexa
3 use module_lin_alg
5 contains
7 subroutine w_assembly( &
8 ifds, ifde, kfds,kfde, jfds, jfde, & ! fire grid dimensions
9 ifms, ifme, kfms,kfme, jfms, jfme, &
10 ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds
11 ifts, ifte, kfts, kfte, jfts,jfte, &
12 lambda,u0, v0, w0, A, X, Y, Z, & !Input from femwind, u0, v0, w0, Spatial Grid Data
13 U,V,W) !U,V,W
14 !Purpose: Create Arrays of Wind Vector Component Values at Center Points of Spatial Grid
15 !In:
16 !A Coefficient Matrix size 3X3, symmetric positive definite
17 !u0, v0, w0 Initial wind speed values in x,y,z direction at grid cell centers
18 !iflags1 iflags1 = 1 returns Kloc and Jg from hexa, iflags2 = 2 returns Floc and Jg from hexa
19 !iflags2 iflags2 =1 indicates add initial wind to calculated wind
20 !out:
21 ! U,V,W Computed wind values in x,y,z direction
23 implicit none
25 !*** arguments
27 integer, intent(in):: &
28 ifds, ifde, kfds,kfde, jfds, jfde, & ! fire grid dimensions
29 ifms, ifme, kfms,kfme, jfms, jfme, &
30 ifps, ifpe, kfps, kfpe, jfps, jfpe, & ! fire patch bounds
31 ifts, ifte, kfts, kfte, jfts,jfte
33 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: lambda
35 !Input for hexa
36 real, intent(in) :: A(3,3)
37 real, intent(in), dimension(ifms:ifme, kfms:kfme, jfms:jfme)::X,Y,Z,u0, v0, w0
39 !Output
40 real, intent(out), dimension(ifms:ifme, kfms:kfme, jfms:jfme)::U,V,W
42 !*** local
44 integer:: ie1, ie2, ie3, ic1, ic2, ic3, iloc, i, k1, k2, k3
45 real :: Kloc(8,8), Floc(8), Jg(8,3)
46 real :: Xloc(3,8), u0loc(3)
47 real :: grad(3)
48 real :: lambda_loc(8)
49 real :: A_inv(3,3)
50 real :: vol
52 !*** executable
54 lambda_loc = 0.
55 Xloc = 99999.
56 Jg = 0.
57 Kloc = 0.
58 Floc = 0.
59 grad = 0.
60 u0loc =0.
62 !*** u0loc is an input for module_hexa, but is not used to construct K. Do I need to define this?
63 !** executable
65 !print *, 'u0 vector is', u0
66 do ie2=jfts,jfte
67 do ie3=kfts, kfte
68 do ie1=ifts, ifte
69 ! constant part
70 do ic2=0,1
71 do ic3=0,1
72 do ic1=0,1
73 Xloc(1,iloc)=X(ie1 + ic1, ie3 + ic3, ie2 + ic2)
74 Xloc(2,iloc)=Y(ie1 + ic1, ie3 + ic3, ie2 + ic2)
75 Xloc(3,iloc)=Z(ie1 + ic1, ie3 + ic3, ie2 + ic2)
76 enddo
77 enddo
78 enddo
79 u0loc(1) = u0(ie1,ie3,ie2)
80 u0loc(2) = v0(ie1, ie3,ie2)
81 u0loc(3) = w0(ie1, ie3,ie2)
82 !fine print *, 'local lambda is', lambda_loc
83 !print* , 'Xloc is', Xloc
84 !print* , 'u0loc is', u0loc(1)
85 call hexa(A,Xloc,u0loc,Kloc,Floc,Jg,vol,3)
86 !print*, 'Jg is', Jg
87 !print*, shape(jg)
88 !*** end of constant part
89 do ic2=0,1
90 do ic3=0,1
91 do ic1=0,1
92 iloc=1+ic1+2*(ic2+2*ic3); !local index of the node in the element
93 k1 = ie1 + ic1
94 k2 = ie2 + ic2
95 k3 = ie3 + ic3
96 lambda_loc(iloc) = lambda(k1,k3,k2)
97 enddo
98 enddo
99 enddo
100 grad = matmul(transpose(Jg),lambda_loc)
101 !not ok print *,'Grad before multiplication by A_inv is', grad
103 call Inv3(A, A_inv)
104 grad = matmul(transpose(A_inv),grad)
105 ! Not ok print *,'Grad after multiplication by A_inv is', grad
107 U(ie1, ie3, ie2)=grad(1)+ u0(ie1, ie3, ie2)
108 V(ie1, ie3, ie2)=grad(2)+ v0(ie1, ie3, ie2)
109 W(ie1, ie3, ie2)=grad(3)+ w0(ie1, ie3, ie2)
111 enddo
112 !print *, 'lambda array', lambda
113 enddo
114 enddo
115 !print *,'Shape of U', shape(U)
117 end subroutine w_assembly
118 end module module_w_assembly