source: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90 @ 658

Last change on this file since 658 was 658, checked in by dubos, 6 years ago

devel/unstructured : updated kernels
File size: 12.0 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- caldyn_slow_NH ----------------------------------
3   !$OMP DO SCHEDULE(STATIC)
4   DO ij = 1, primal_num
5      kdown = 1
6      kup = 1
7      l=1
8      w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij)
9      !DIR$ SIMD
10      DO l = 2, llm
11         kdown = l-1
12         kup = l
13         w_il(l,ij) = 2.*W(l,ij)/(rhodz(kdown,ij)+rhodz(kup,ij))
14      END DO
15      kdown = llm
16      kup = llm
17      l=llm+1
18      w_il(l,ij) = 2.*W(l,ij)/rhodz(kdown,ij)
19   END DO
20   !$OMP END DO
21   !$OMP DO SCHEDULE(STATIC)
22   DO edge = 1, edge_num
23      ij_left = left(edge)
24      ij_right = right(edge)
25      kdown = 1
26      kup = 1
27      l=1
28      ! compute DePhi, v_el, G_el, F_el
29      ! v_el, W2_el and therefore G_el incorporate metric factor le_de
30      ! while DePhil, W_el and F_el do not
31      W_el = .5*( W(l,ij_right)+W(l,ij_left) )
32      DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left))
33      F_el(l,edge) = DePhil(l,edge)*W_el
34      W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) )
35      v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
36      G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
37      !DIR$ SIMD
38      DO l = 2, llm
39         kdown = l-1
40         kup = l
41         ! compute DePhi, v_el, G_el, F_el
42         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
43         ! while DePhil, W_el and F_el do not
44         W_el = .5*( W(l,ij_right)+W(l,ij_left) )
45         DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left))
46         F_el(l,edge) = DePhil(l,edge)*W_el
47         W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) )
48         v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
49         G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
50      END DO
51      kdown = llm
52      kup = llm
53      l=llm+1
54      ! compute DePhi, v_el, G_el, F_el
55      ! v_el, W2_el and therefore G_el incorporate metric factor le_de
56      ! while DePhil, W_el and F_el do not
57      W_el = .5*( W(l,ij_right)+W(l,ij_left) )
58      DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left))
59      F_el(l,edge) = DePhil(l,edge)*W_el
60      W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) )
61      v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
62      G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
63   END DO
64   !$OMP END DO
65   ! compute GradPhi2, dPhi, dW
66   !$OMP DO SCHEDULE(STATIC)
67   DO ij = 1, primal_num
68      ! this VLOOP iterates over primal cell edges
69      SELECT CASE(primal_deg(ij))
70      CASE(4)
71         edge1 = primal_edge(1,ij)
72         edge2 = primal_edge(2,ij)
73         edge3 = primal_edge(3,ij)
74         edge4 = primal_edge(4,ij)
75         le_de1 = le_de(edge1)
76         le_de2 = le_de(edge2)
77         le_de3 = le_de(edge3)
78         le_de4 = le_de(edge4)
79         sign1 = primal_ne(1,ij)
80         sign2 = primal_ne(2,ij)
81         sign3 = primal_ne(3,ij)
82         sign4 = primal_ne(4,ij)
83         !DIR$ SIMD
84         DO l = 1, llm+1
85            gPhi2=0.
86            dP=0.
87            divG=0
88            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2
89            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1)
90            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de
91            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2
92            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2)
93            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de
94            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2
95            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3)
96            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de
97            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2
98            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4)
99            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de
100            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
101            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
102            dW(l,ij) = (-1./Ai(ij))*divG
103         END DO
104      CASE(5)
105         edge1 = primal_edge(1,ij)
106         edge2 = primal_edge(2,ij)
107         edge3 = primal_edge(3,ij)
108         edge4 = primal_edge(4,ij)
109         edge5 = primal_edge(5,ij)
110         le_de1 = le_de(edge1)
111         le_de2 = le_de(edge2)
112         le_de3 = le_de(edge3)
113         le_de4 = le_de(edge4)
114         le_de5 = le_de(edge5)
115         sign1 = primal_ne(1,ij)
116         sign2 = primal_ne(2,ij)
117         sign3 = primal_ne(3,ij)
118         sign4 = primal_ne(4,ij)
119         sign5 = primal_ne(5,ij)
120         !DIR$ SIMD
121         DO l = 1, llm+1
122            gPhi2=0.
123            dP=0.
124            divG=0
125            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2
126            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1)
127            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de
128            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2
129            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2)
130            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de
131            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2
132            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3)
133            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de
134            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2
135            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4)
136            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de
137            gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2
138            dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5)
139            divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de
140            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
141            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
142            dW(l,ij) = (-1./Ai(ij))*divG
143         END DO
144      CASE(6)
145         edge1 = primal_edge(1,ij)
146         edge2 = primal_edge(2,ij)
147         edge3 = primal_edge(3,ij)
148         edge4 = primal_edge(4,ij)
149         edge5 = primal_edge(5,ij)
150         edge6 = primal_edge(6,ij)
151         le_de1 = le_de(edge1)
152         le_de2 = le_de(edge2)
153         le_de3 = le_de(edge3)
154         le_de4 = le_de(edge4)
155         le_de5 = le_de(edge5)
156         le_de6 = le_de(edge6)
157         sign1 = primal_ne(1,ij)
158         sign2 = primal_ne(2,ij)
159         sign3 = primal_ne(3,ij)
160         sign4 = primal_ne(4,ij)
161         sign5 = primal_ne(5,ij)
162         sign6 = primal_ne(6,ij)
163         !DIR$ SIMD
164         DO l = 1, llm+1
165            gPhi2=0.
166            dP=0.
167            divG=0
168            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2
169            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1)
170            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de
171            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2
172            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2)
173            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de
174            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2
175            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3)
176            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de
177            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2
178            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4)
179            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de
180            gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2
181            dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5)
182            divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de
183            gPhi2 = gPhi2 + le_de6*DePhil(l,edge6)**2
184            dP = dP + le_de6*DePhil(l,edge6)*v_el(l,edge6)
185            divG = divG + sign6*G_el(l,edge6) ! -div(G_el), G_el already has le_de
186            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
187            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
188            dW(l,ij) = (-1./Ai(ij))*divG
189         END DO
190      CASE DEFAULT
191         !DIR$ SIMD
192         DO l = 1, llm+1
193            gPhi2=0.
194            dP=0.
195            divG=0
196            DO iedge = 1, primal_deg(ij)
197               edge = primal_edge(iedge,ij)
198               gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2
199               dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge)
200               divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de
201            END DO
202            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
203            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
204            dW(l,ij) = (-1./Ai(ij))*divG
205         END DO
206      END SELECT
207   END DO
208   !$OMP END DO
209   ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
210   !$OMP BARRIER
211   ! Compute berni at scalar points
212   !$OMP DO SCHEDULE(STATIC)
213   DO ij = 1, primal_num
214      ! this VLOOP iterates over primal cell edges
215      SELECT CASE(primal_deg(ij))
216      CASE(4)
217         edge1 = primal_edge(1,ij)
218         edge2 = primal_edge(2,ij)
219         edge3 = primal_edge(3,ij)
220         edge4 = primal_edge(4,ij)
221         le_de1 = le_de(edge1)
222         le_de2 = le_de(edge2)
223         le_de3 = le_de(edge3)
224         le_de4 = le_de(edge4)
225         !DIR$ SIMD
226         DO l = 1, llm
227            u2=0.
228            u2 = u2 + le_de1*u(l,edge1)**2
229            u2 = u2 + le_de2*u(l,edge2)**2
230            u2 = u2 + le_de3*u(l,edge3)**2
231            u2 = u2 + le_de4*u(l,edge4)**2
232            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
233         END DO
234      CASE(5)
235         edge1 = primal_edge(1,ij)
236         edge2 = primal_edge(2,ij)
237         edge3 = primal_edge(3,ij)
238         edge4 = primal_edge(4,ij)
239         edge5 = primal_edge(5,ij)
240         le_de1 = le_de(edge1)
241         le_de2 = le_de(edge2)
242         le_de3 = le_de(edge3)
243         le_de4 = le_de(edge4)
244         le_de5 = le_de(edge5)
245         !DIR$ SIMD
246         DO l = 1, llm
247            u2=0.
248            u2 = u2 + le_de1*u(l,edge1)**2
249            u2 = u2 + le_de2*u(l,edge2)**2
250            u2 = u2 + le_de3*u(l,edge3)**2
251            u2 = u2 + le_de4*u(l,edge4)**2
252            u2 = u2 + le_de5*u(l,edge5)**2
253            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
254         END DO
255      CASE(6)
256         edge1 = primal_edge(1,ij)
257         edge2 = primal_edge(2,ij)
258         edge3 = primal_edge(3,ij)
259         edge4 = primal_edge(4,ij)
260         edge5 = primal_edge(5,ij)
261         edge6 = primal_edge(6,ij)
262         le_de1 = le_de(edge1)
263         le_de2 = le_de(edge2)
264         le_de3 = le_de(edge3)
265         le_de4 = le_de(edge4)
266         le_de5 = le_de(edge5)
267         le_de6 = le_de(edge6)
268         !DIR$ SIMD
269         DO l = 1, llm
270            u2=0.
271            u2 = u2 + le_de1*u(l,edge1)**2
272            u2 = u2 + le_de2*u(l,edge2)**2
273            u2 = u2 + le_de3*u(l,edge3)**2
274            u2 = u2 + le_de4*u(l,edge4)**2
275            u2 = u2 + le_de5*u(l,edge5)**2
276            u2 = u2 + le_de6*u(l,edge6)**2
277            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
278         END DO
279      CASE DEFAULT
280         !DIR$ SIMD
281         DO l = 1, llm
282            u2=0.
283            DO iedge = 1, primal_deg(ij)
284               edge = primal_edge(iedge,ij)
285               u2 = u2 + le_de(edge)*u(l,edge)**2
286            END DO
287            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
288         END DO
289      END SELECT
290   END DO
291   !$OMP END DO
292   !$OMP DO SCHEDULE(STATIC)
293   DO edge = 1, edge_num
294      ij_left = left(edge)
295      ij_right = right(edge)
296      !DIR$ SIMD
297      DO l = 1, llm
298         ! Compute mass flux and grad(berni)
299         uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge) - .5*( F_el(l,edge)+F_el(l+1,edge) )
300         hflux(l,edge) = le_de(edge)*uu
301         du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right))
302      END DO
303   END DO
304   !$OMP END DO
305   !---------------------------- caldyn_slow_NH ----------------------------------
306   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.