Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

coriolis.k90

Last change on this file was 940, checked in by dubos, 5 years ago
devel : DySL for enstrophy-conserving scheme
File size: 12.4 KB

Line
1	!--------------------------------------------------------------------------
2	!---------------------------- coriolis ----------------------------------
3	!
4	DO iq=1,nqdyn
5	!$OMP DO SCHEDULE(STATIC)
6	DO edge = 1, edge_num
7	ij_left = left(edge)
8	ij_right = right(edge)
9	!DIR$ SIMD
10	DO l = 1, llm
11	Ftheta(l,edge) = .5(theta(l,ij_left,iq)+theta(l,ij_right,iq))hflux(l,edge)
12	END DO
13	END DO
14	!$OMP END DO
15	!$OMP DO SCHEDULE(STATIC)
16	DO ij = 1, primal_num
17	! this VLOOP iterates over primal cell edges
18	SELECT CASE(primal_deg(ij))
19	CASE(4)
20	edge1 = primal_edge(1,ij)
21	edge2 = primal_edge(2,ij)
22	edge3 = primal_edge(3,ij)
23	edge4 = primal_edge(4,ij)
24	sign1 = primal_ne(1,ij)
25	sign2 = primal_ne(2,ij)
26	sign3 = primal_ne(3,ij)
27	sign4 = primal_ne(4,ij)
28	!DIR$ SIMD
29	DO l = 1, llm
30	divF=0.
31	divF = divF + Ftheta(l,edge1)*sign1
32	divF = divF + Ftheta(l,edge2)*sign2
33	divF = divF + Ftheta(l,edge3)*sign3
34	divF = divF + Ftheta(l,edge4)*sign4
35	dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
36	END DO
37	CASE(5)
38	edge1 = primal_edge(1,ij)
39	edge2 = primal_edge(2,ij)
40	edge3 = primal_edge(3,ij)
41	edge4 = primal_edge(4,ij)
42	edge5 = primal_edge(5,ij)
43	sign1 = primal_ne(1,ij)
44	sign2 = primal_ne(2,ij)
45	sign3 = primal_ne(3,ij)
46	sign4 = primal_ne(4,ij)
47	sign5 = primal_ne(5,ij)
48	!DIR$ SIMD
49	DO l = 1, llm
50	divF=0.
51	divF = divF + Ftheta(l,edge1)*sign1
52	divF = divF + Ftheta(l,edge2)*sign2
53	divF = divF + Ftheta(l,edge3)*sign3
54	divF = divF + Ftheta(l,edge4)*sign4
55	divF = divF + Ftheta(l,edge5)*sign5
56	dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
57	END DO
58	CASE(6)
59	edge1 = primal_edge(1,ij)
60	edge2 = primal_edge(2,ij)
61	edge3 = primal_edge(3,ij)
62	edge4 = primal_edge(4,ij)
63	edge5 = primal_edge(5,ij)
64	edge6 = primal_edge(6,ij)
65	sign1 = primal_ne(1,ij)
66	sign2 = primal_ne(2,ij)
67	sign3 = primal_ne(3,ij)
68	sign4 = primal_ne(4,ij)
69	sign5 = primal_ne(5,ij)
70	sign6 = primal_ne(6,ij)
71	!DIR$ SIMD
72	DO l = 1, llm
73	divF=0.
74	divF = divF + Ftheta(l,edge1)*sign1
75	divF = divF + Ftheta(l,edge2)*sign2
76	divF = divF + Ftheta(l,edge3)*sign3
77	divF = divF + Ftheta(l,edge4)*sign4
78	divF = divF + Ftheta(l,edge5)*sign5
79	divF = divF + Ftheta(l,edge6)*sign6
80	dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
81	END DO
82	CASE DEFAULT
83	!DIR$ SIMD
84	DO l = 1, llm
85	divF=0.
86	DO iedge = 1, primal_deg(ij)
87	edge = primal_edge(iedge,ij)
88	divF = divF + Ftheta(l,edge)*primal_ne(iedge,ij)
89	END DO
90	dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
91	END DO
92	END SELECT
93	END DO
94	!$OMP END DO
95	END DO ! iq
96	!
97	!$OMP DO SCHEDULE(STATIC)
98	DO ij = 1, primal_num
99	! this VLOOP iterates over primal cell edges
100	SELECT CASE(primal_deg(ij))
101	CASE(4)
102	edge1 = primal_edge(1,ij)
103	edge2 = primal_edge(2,ij)
104	edge3 = primal_edge(3,ij)
105	edge4 = primal_edge(4,ij)
106	sign1 = primal_ne(1,ij)
107	sign2 = primal_ne(2,ij)
108	sign3 = primal_ne(3,ij)
109	sign4 = primal_ne(4,ij)
110	!DIR$ SIMD
111	DO l = 1, llm
112	divF=0.
113	divF = divF + hflux(l,edge1)*sign1
114	divF = divF + hflux(l,edge2)*sign2
115	divF = divF + hflux(l,edge3)*sign3
116	divF = divF + hflux(l,edge4)*sign4
117	convm(l,ij) = -divF / Ai(ij)
118	END DO
119	CASE(5)
120	edge1 = primal_edge(1,ij)
121	edge2 = primal_edge(2,ij)
122	edge3 = primal_edge(3,ij)
123	edge4 = primal_edge(4,ij)
124	edge5 = primal_edge(5,ij)
125	sign1 = primal_ne(1,ij)
126	sign2 = primal_ne(2,ij)
127	sign3 = primal_ne(3,ij)
128	sign4 = primal_ne(4,ij)
129	sign5 = primal_ne(5,ij)
130	!DIR$ SIMD
131	DO l = 1, llm
132	divF=0.
133	divF = divF + hflux(l,edge1)*sign1
134	divF = divF + hflux(l,edge2)*sign2
135	divF = divF + hflux(l,edge3)*sign3
136	divF = divF + hflux(l,edge4)*sign4
137	divF = divF + hflux(l,edge5)*sign5
138	convm(l,ij) = -divF / Ai(ij)
139	END DO
140	CASE(6)
141	edge1 = primal_edge(1,ij)
142	edge2 = primal_edge(2,ij)
143	edge3 = primal_edge(3,ij)
144	edge4 = primal_edge(4,ij)
145	edge5 = primal_edge(5,ij)
146	edge6 = primal_edge(6,ij)
147	sign1 = primal_ne(1,ij)
148	sign2 = primal_ne(2,ij)
149	sign3 = primal_ne(3,ij)
150	sign4 = primal_ne(4,ij)
151	sign5 = primal_ne(5,ij)
152	sign6 = primal_ne(6,ij)
153	!DIR$ SIMD
154	DO l = 1, llm
155	divF=0.
156	divF = divF + hflux(l,edge1)*sign1
157	divF = divF + hflux(l,edge2)*sign2
158	divF = divF + hflux(l,edge3)*sign3
159	divF = divF + hflux(l,edge4)*sign4
160	divF = divF + hflux(l,edge5)*sign5
161	divF = divF + hflux(l,edge6)*sign6
162	convm(l,ij) = -divF / Ai(ij)
163	END DO
164	CASE DEFAULT
165	!DIR$ SIMD
166	DO l = 1, llm
167	divF=0.
168	DO iedge = 1, primal_deg(ij)
169	edge = primal_edge(iedge,ij)
170	divF = divF + hflux(l,edge)*primal_ne(iedge,ij)
171	END DO
172	convm(l,ij) = -divF / Ai(ij)
173	END DO
174	END SELECT
175	END DO
176	!$OMP END DO
177	!
178	SELECT CASE(caldyn_conserv)
179	CASE(conserv_energy) ! energy-conserving TRiSK
180	!$OMP DO SCHEDULE(STATIC)
181	DO edge = 1, edge_num
182	! this VLOOP iterates over the TRISK stencil
183	SELECT CASE(trisk_deg(edge))
184	CASE(4)
185	!DIR$ SIMD
186	DO l = 1, llm
187	du_trisk=0.
188	itrisk = 1
189	edge_trisk = trisk(1,edge)
190	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
191	itrisk = 2
192	edge_trisk = trisk(2,edge)
193	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
194	itrisk = 3
195	edge_trisk = trisk(3,edge)
196	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
197	itrisk = 4
198	edge_trisk = trisk(4,edge)
199	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
200	du(l,edge) = du(l,edge) + .5*du_trisk
201	END DO
202	CASE(10)
203	!DIR$ SIMD
204	DO l = 1, llm
205	du_trisk=0.
206	itrisk = 1
207	edge_trisk = trisk(1,edge)
208	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
209	itrisk = 2
210	edge_trisk = trisk(2,edge)
211	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
212	itrisk = 3
213	edge_trisk = trisk(3,edge)
214	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
215	itrisk = 4
216	edge_trisk = trisk(4,edge)
217	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
218	itrisk = 5
219	edge_trisk = trisk(5,edge)
220	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
221	itrisk = 6
222	edge_trisk = trisk(6,edge)
223	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
224	itrisk = 7
225	edge_trisk = trisk(7,edge)
226	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
227	itrisk = 8
228	edge_trisk = trisk(8,edge)
229	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
230	itrisk = 9
231	edge_trisk = trisk(9,edge)
232	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
233	itrisk = 10
234	edge_trisk = trisk(10,edge)
235	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
236	du(l,edge) = du(l,edge) + .5*du_trisk
237	END DO
238	CASE DEFAULT
239	!DIR$ SIMD
240	DO l = 1, llm
241	du_trisk=0.
242	DO itrisk = 1, trisk_deg(edge)
243	edge_trisk = trisk(itrisk,edge)
244	du_trisk = du_trisk + wee(itrisk,edge,1)hflux(l,edge_trisk)(qu(l,edge)+qu(l,edge_trisk))
245	END DO
246	du(l,edge) = du(l,edge) + .5*du_trisk
247	END DO
248	END SELECT
249	END DO
250	!$OMP END DO
251	CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
252	!$OMP DO SCHEDULE(STATIC)
253	DO edge = 1, edge_num
254	! this VLOOP iterates over the TRISK stencil
255	SELECT CASE(trisk_deg(edge))
256	CASE(4)
257	!DIR$ SIMD
258	DO l = 1, llm
259	du_trisk=0.
260	itrisk = 1
261	edge_trisk = trisk(1,edge)
262	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
263	itrisk = 2
264	edge_trisk = trisk(2,edge)
265	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
266	itrisk = 3
267	edge_trisk = trisk(3,edge)
268	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
269	itrisk = 4
270	edge_trisk = trisk(4,edge)
271	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
272	du(l,edge) = du(l,edge) + du_trisk*qu(l,edge)
273	END DO
274	CASE(10)
275	!DIR$ SIMD
276	DO l = 1, llm
277	du_trisk=0.
278	itrisk = 1
279	edge_trisk = trisk(1,edge)
280	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
281	itrisk = 2
282	edge_trisk = trisk(2,edge)
283	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
284	itrisk = 3
285	edge_trisk = trisk(3,edge)
286	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
287	itrisk = 4
288	edge_trisk = trisk(4,edge)
289	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
290	itrisk = 5
291	edge_trisk = trisk(5,edge)
292	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
293	itrisk = 6
294	edge_trisk = trisk(6,edge)
295	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
296	itrisk = 7
297	edge_trisk = trisk(7,edge)
298	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
299	itrisk = 8
300	edge_trisk = trisk(8,edge)
301	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
302	itrisk = 9
303	edge_trisk = trisk(9,edge)
304	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
305	itrisk = 10
306	edge_trisk = trisk(10,edge)
307	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
308	du(l,edge) = du(l,edge) + du_trisk*qu(l,edge)
309	END DO
310	CASE DEFAULT
311	!DIR$ SIMD
312	DO l = 1, llm
313	du_trisk=0.
314	DO itrisk = 1, trisk_deg(edge)
315	edge_trisk = trisk(itrisk,edge)
316	du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
317	END DO
318	du(l,edge) = du(l,edge) + du_trisk*qu(l,edge)
319	END DO
320	END SELECT
321	END DO
322	!$OMP END DO
323	END SELECT
324	!---------------------------- coriolis ----------------------------------
325	!--------------------------------------------------------------------------

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format