source: codes/icosagcm/devel/src/kernels_unst/wind_centered.k90 @ 913

Last change on this file since 913 was 796, checked in by jisesh, 5 years ago

unstructured : reconstruction of velocity components (bugged)

File size: 8.1 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- wind_centered ----------------------------------
3   ! Perot reconstruction based on Gauss theorem
4   ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
5   !$OMP DO SCHEDULE(STATIC)
6   DO ij = 1, primal_num
7      ! this VLOOP iterates over primal cell edges
8      SELECT CASE(primal_deg(ij))
9      CASE(4)
10         edge1 = primal_edge(1,ij)
11         edge2 = primal_edge(2,ij)
12         edge3 = primal_edge(3,ij)
13         edge4 = primal_edge(4,ij)
14         sign1 = primal_ne(1,ij)
15         sign2 = primal_ne(2,ij)
16         sign3 = primal_ne(3,ij)
17         sign4 = primal_ne(4,ij)
18         ij_up1 = up(edge1)
19         ij_up2 = up(edge2)
20         ij_up3 = up(edge3)
21         ij_up4 = up(edge4)
22         ij_down1 = down(edge1)
23         ij_down2 = down(edge2)
24         ij_down3 = down(edge3)
25         ij_down4 = down(edge4)
26         !DIR$ SIMD
27         DO l = 1, llm
28            cx=centroid(ij,1)
29            cy=centroid(ij,2)
30            cz=centroid(ij,3)
31            ux=0. ; uy=0. ; uz=0.
32            ue_le = sign1*ue(l,edge1)*le(edge1)
33            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
34            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
35            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
36            ue_le = sign2*ue(l,edge2)*le(edge2)
37            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
38            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
39            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
40            ue_le = sign3*ue(l,edge3)*le(edge3)
41            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
42            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
43            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
44            ue_le = sign4*ue(l,edge4)*le(edge4)
45            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
46            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
47            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
48            fac = scale*(1./Ai(ij))
49            ucenter(l,ij,1)=ux*fac
50            ucenter(l,ij,2)=uy*fac
51            ucenter(l,ij,3)=uz*fac
52         END DO
53      CASE(5)
54         edge1 = primal_edge(1,ij)
55         edge2 = primal_edge(2,ij)
56         edge3 = primal_edge(3,ij)
57         edge4 = primal_edge(4,ij)
58         edge5 = primal_edge(5,ij)
59         sign1 = primal_ne(1,ij)
60         sign2 = primal_ne(2,ij)
61         sign3 = primal_ne(3,ij)
62         sign4 = primal_ne(4,ij)
63         sign5 = primal_ne(5,ij)
64         ij_up1 = up(edge1)
65         ij_up2 = up(edge2)
66         ij_up3 = up(edge3)
67         ij_up4 = up(edge4)
68         ij_up5 = up(edge5)
69         ij_down1 = down(edge1)
70         ij_down2 = down(edge2)
71         ij_down3 = down(edge3)
72         ij_down4 = down(edge4)
73         ij_down5 = down(edge5)
74         !DIR$ SIMD
75         DO l = 1, llm
76            cx=centroid(ij,1)
77            cy=centroid(ij,2)
78            cz=centroid(ij,3)
79            ux=0. ; uy=0. ; uz=0.
80            ue_le = sign1*ue(l,edge1)*le(edge1)
81            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
82            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
83            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
84            ue_le = sign2*ue(l,edge2)*le(edge2)
85            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
86            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
87            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
88            ue_le = sign3*ue(l,edge3)*le(edge3)
89            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
90            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
91            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
92            ue_le = sign4*ue(l,edge4)*le(edge4)
93            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
94            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
95            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
96            ue_le = sign5*ue(l,edge5)*le(edge5)
97            ux = ux + ue_le*(.5*(xyz_v(ij_up5,1)+xyz_v(ij_down5,1))-cx)
98            uy = uy + ue_le*(.5*(xyz_v(ij_up5,2)+xyz_v(ij_down5,2))-cy)
99            uz = uz + ue_le*(.5*(xyz_v(ij_up5,3)+xyz_v(ij_down5,3))-cz)
100            fac = scale*(1./Ai(ij))
101            ucenter(l,ij,1)=ux*fac
102            ucenter(l,ij,2)=uy*fac
103            ucenter(l,ij,3)=uz*fac
104         END DO
105      CASE(6)
106         edge1 = primal_edge(1,ij)
107         edge2 = primal_edge(2,ij)
108         edge3 = primal_edge(3,ij)
109         edge4 = primal_edge(4,ij)
110         edge5 = primal_edge(5,ij)
111         edge6 = primal_edge(6,ij)
112         sign1 = primal_ne(1,ij)
113         sign2 = primal_ne(2,ij)
114         sign3 = primal_ne(3,ij)
115         sign4 = primal_ne(4,ij)
116         sign5 = primal_ne(5,ij)
117         sign6 = primal_ne(6,ij)
118         ij_up1 = up(edge1)
119         ij_up2 = up(edge2)
120         ij_up3 = up(edge3)
121         ij_up4 = up(edge4)
122         ij_up5 = up(edge5)
123         ij_up6 = up(edge6)
124         ij_down1 = down(edge1)
125         ij_down2 = down(edge2)
126         ij_down3 = down(edge3)
127         ij_down4 = down(edge4)
128         ij_down5 = down(edge5)
129         ij_down6 = down(edge6)
130         !DIR$ SIMD
131         DO l = 1, llm
132            cx=centroid(ij,1)
133            cy=centroid(ij,2)
134            cz=centroid(ij,3)
135            ux=0. ; uy=0. ; uz=0.
136            ue_le = sign1*ue(l,edge1)*le(edge1)
137            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
138            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
139            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
140            ue_le = sign2*ue(l,edge2)*le(edge2)
141            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
142            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
143            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
144            ue_le = sign3*ue(l,edge3)*le(edge3)
145            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
146            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
147            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
148            ue_le = sign4*ue(l,edge4)*le(edge4)
149            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
150            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
151            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
152            ue_le = sign5*ue(l,edge5)*le(edge5)
153            ux = ux + ue_le*(.5*(xyz_v(ij_up5,1)+xyz_v(ij_down5,1))-cx)
154            uy = uy + ue_le*(.5*(xyz_v(ij_up5,2)+xyz_v(ij_down5,2))-cy)
155            uz = uz + ue_le*(.5*(xyz_v(ij_up5,3)+xyz_v(ij_down5,3))-cz)
156            ue_le = sign6*ue(l,edge6)*le(edge6)
157            ux = ux + ue_le*(.5*(xyz_v(ij_up6,1)+xyz_v(ij_down6,1))-cx)
158            uy = uy + ue_le*(.5*(xyz_v(ij_up6,2)+xyz_v(ij_down6,2))-cy)
159            uz = uz + ue_le*(.5*(xyz_v(ij_up6,3)+xyz_v(ij_down6,3))-cz)
160            fac = scale*(1./Ai(ij))
161            ucenter(l,ij,1)=ux*fac
162            ucenter(l,ij,2)=uy*fac
163            ucenter(l,ij,3)=uz*fac
164         END DO
165      CASE DEFAULT
166         !DIR$ SIMD
167         DO l = 1, llm
168            cx=centroid(ij,1)
169            cy=centroid(ij,2)
170            cz=centroid(ij,3)
171            ux=0. ; uy=0. ; uz=0.
172            DO iedge = 1, primal_deg(ij)
173               edge = primal_edge(iedge,ij)
174               ij_up = up(edge)
175               ij_down = down(edge)
176               ue_le = primal_ne(iedge,ij)*ue(l,edge)*le(edge)
177               ux = ux + ue_le*(.5*(xyz_v(ij_up,1)+xyz_v(ij_down,1))-cx)
178               uy = uy + ue_le*(.5*(xyz_v(ij_up,2)+xyz_v(ij_down,2))-cy)
179               uz = uz + ue_le*(.5*(xyz_v(ij_up,3)+xyz_v(ij_down,3))-cz)
180            END DO
181            fac = scale*(1./Ai(ij))
182            ucenter(l,ij,1)=ux*fac
183            ucenter(l,ij,2)=uy*fac
184            ucenter(l,ij,3)=uz*fac
185         END DO
186      END SELECT
187   END DO
188   !$OMP END DO
189   !---------------------------- wind_centered ----------------------------------
190   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.