Context Navigation

advect.F90 @ 954

Last change on this file since 954 was 954, checked in by adurocher, 5 years ago

trunk : Added metric terms to kernels parameters to avoid Host/GPU transferts

Metric terms are now subroutine parameters instead of module variables in kernel subroutines. Dummy arguments for metric terms are now defined as fixed-size arrays, and arrays dimensions are well known when entering an 'acc data' region. Array descriptors are no longer transferred form host to device each time the data region is executed.

File size: 19.5 KB

Line
1	MODULE advect_mod
2	USE icosa
3	IMPLICIT NONE
4
5	PRIVATE
6
7	PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz
8
9	CONTAINS
10
11	!==========================================================================
12
13	SUBROUTINE init_advect(normal,tangent,sqrt_leng)
14	IMPLICIT NONE
15	REAL(rstd),INTENT(OUT) :: normal(3iimjjm,3)
16	REAL(rstd),INTENT(OUT) :: tangent(3iimjjm,3)
17	REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm)
18
19	INTEGER :: ij
20
21	!DIR$ SIMD
22	DO ij=ij_begin,ij_end
23
24	CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:))
25	normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)*2)+1e-50)ne(ij,right)
26
27	CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:))
28	normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)*2)+1e-50)ne(ij,lup)
29
30	CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:))
31	normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)*2)+1e-50)ne(ij,ldown)
32
33	tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:)
34	tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50)
35
36	tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:)
37	tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50)
38
39	tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:)
40	tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
41
42	sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))2), &
43	sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))2), &
44	sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))2)) )
45	ENDDO
46
47	END SUBROUTINE init_advect
48
49	!=======================================================================================
50
51	SUBROUTINE compute_gradq3d(qi,sqrt_leng,gradq3d,xyz_i,xyz_v)
52	USE trace
53	USE omp_para
54	IMPLICIT NONE
55	REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm)
56	REAL(rstd),INTENT(IN) :: sqrt_leng(iim*jjm)
57	REAL(rstd),INTENT(IN) :: xyz_i(iim*jjm,3)
58	REAL(rstd),INTENT(IN) :: xyz_v(2iimjjm,3)
59	REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3)
60	REAL(rstd) :: maxq,minq,minq_c,maxq_c
61	REAL(rstd) :: alphamx,alphami,alpha,maggrd
62	REAL(rstd) :: arr(2iimjjm)
63	REAL(rstd) :: ar(iim*jjm)
64	REAL(rstd) :: gradtri(2iimjjm,llm,3)
65	INTEGER :: ij,k,l
66	REAL(rstd) :: detx,dety,detz,det
67	REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33
68	REAL(rstd) :: x1,x2,x3
69	REAL(rstd) :: dq(3)
70
71	CALL trace_start("compute_gradq3d1")
72
73	! TODO : precompute ar, drop arr as output argument of gradq ?
74
75	!========================================================================================== GRADIENT
76	! Compute gradient at triangles solving a linear system
77	! arr = area of triangle joining centroids of hexagons
78	! DO l = ll_begin,ll_end
79	!!DIR$ SIMD
80	! DO ij=ij_begin_ext,ij_end_ext
81	!! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
82	!! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
83	! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
84	! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
85	! END DO
86	! END DO
87
88	!$acc data create(gradtri(:,:,:), arr(:), ar(:)) present(sqrt_leng(:), xyz_i(:,:), xyz_v(:,:), qi(:,:), gradq3d(:,:,:)) async
89
90	!$acc parallel loop collapse(2) async private(A, dq)
91	DO l = ll_begin,ll_end
92	!DIR$ SIMD
93	DO ij=ij_begin_ext,ij_end_ext
94	! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
95
96
97	A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1); A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3)
98	A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1); A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3)
99	A(3,1)=xyz_v(ij+z_up,1); A(3,2)= xyz_v(ij+z_up,2); A(3,3)=xyz_v(ij+z_up,3)
100
101	dq(1) = qi(ij+t_rup,l)-qi(ij,l)
102	dq(2) = qi(ij+t_lup,l)-qi(ij,l)
103	dq(3) = 0.0
104
105
106	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
107
108	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
109	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
110	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
111
112	x1 = a11 * (a22 * a33 - a23 * a32)
113	x2 = a12 * (a21 * a33 - a23 * a31)
114	x3 = a13 * (a21 * a32 - a22 * a31)
115	det = x1 - x2 + x3
116
117	! CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
118
119	a11=dq(1) ; a12=dq(2) ; a13=dq(3)
120	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
121	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
122
123	x1 = a11 * (a22 * a33 - a23 * a32)
124	x2 = a12 * (a21 * a33 - a23 * a31)
125	x3 = a13 * (a21 * a32 - a22 * a31)
126	detx = x1 - x2 + x3
127
128	! CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
129
130	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
131	a21=dq(1) ; a22=dq(2) ; a23=dq(3)
132	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
133
134	x1 = a11 * (a22 * a33 - a23 * a32)
135	x2 = a12 * (a21 * a33 - a23 * a31)
136	x3 = a13 * (a21 * a32 - a22 * a31)
137	dety = x1 - x2 + x3
138
139	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
140
141	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
142	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
143	a31=dq(1) ; a32=dq(2) ; a33=dq(3)
144
145	x1 = a11 * (a22 * a33 - a23 * a32)
146	x2 = a12 * (a21 * a33 - a23 * a31)
147	x3 = a13 * (a21 * a32 - a22 * a31)
148	detz = x1 - x2 + x3
149
150	gradtri(ij+z_up,l,1) = detx
151	gradtri(ij+z_up,l,2) = dety
152	gradtri(ij+z_up,l,3) = detz
153	arr(ij+z_up) = det
154
155	ENDDO
156	ENDDO
157
158	!$acc parallel loop collapse(2) async private(A, dq)
159	DO l = ll_begin,ll_end
160	DO ij=ij_begin_ext,ij_end_ext
161
162
163	! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
164
165	A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1); A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3)
166	A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1); A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3)
167	A(3,1)=xyz_v(ij+z_down,1); A(3,2)= xyz_v(ij+z_down,2); A(3,3)=xyz_v(ij+z_down,3)
168
169	dq(1) = qi(ij+t_ldown,l)-qi(ij,l)
170	dq(2) = qi(ij+t_rdown,l)-qi(ij,l)
171	dq(3) = 0.0
172
173
174	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
175
176	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
177	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
178	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
179
180	x1 = a11 * (a22 * a33 - a23 * a32)
181	x2 = a12 * (a21 * a33 - a23 * a31)
182	x3 = a13 * (a21 * a32 - a22 * a31)
183	det = x1 - x2 + x3
184
185	! CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
186
187	a11=dq(1) ; a12=dq(2) ; a13=dq(3)
188	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
189	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
190
191	x1 = a11 * (a22 * a33 - a23 * a32)
192	x2 = a12 * (a21 * a33 - a23 * a31)
193	x3 = a13 * (a21 * a32 - a22 * a31)
194	detx = x1 - x2 + x3
195
196	! CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
197
198	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
199	a21=dq(1) ; a22=dq(2) ; a23=dq(3)
200	a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
201
202	x1 = a11 * (a22 * a33 - a23 * a32)
203	x2 = a12 * (a21 * a33 - a23 * a31)
204	x3 = a13 * (a21 * a32 - a22 * a31)
205	dety = x1 - x2 + x3
206
207	! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
208
209	a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
210	a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
211	a31=dq(1) ; a32=dq(2) ; a33=dq(3)
212
213	x1 = a11 * (a22 * a33 - a23 * a32)
214	x2 = a12 * (a21 * a33 - a23 * a31)
215	x3 = a13 * (a21 * a32 - a22 * a31)
216	detz = x1 - x2 + x3
217
218	gradtri(ij+z_down,l,1) = detx
219	gradtri(ij+z_down,l,2) = dety
220	gradtri(ij+z_down,l,3) = detz
221	arr(ij+z_down) = det
222
223	END DO
224	END DO
225
226	!DIR$ SIMD
227	!$acc parallel loop async
228	DO ij=ij_begin,ij_end
229	ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
230	ENDDO
231	CALL trace_end("compute_gradq3d1")
232	CALL trace_start2("compute_gradq3d2")
233
234	!$acc parallel loop collapse(3) async
235	DO k=1,3
236	DO l =ll_begin,ll_end
237	!DIR$ SIMD
238	DO ij=ij_begin,ij_end
239	gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k) + &
240	gradtri(ij+z_rup,l,k) + gradtri(ij+z_ldown,l,k) + &
241	gradtri(ij+z_lup,l,k)+ gradtri(ij+z_rdown,l,k) ) / ar(ij)
242	END DO
243	END DO
244	ENDDO
245	CALL trace_end2("compute_gradq3d2")
246	CALL trace_start("compute_gradq3d3")
247
248	!============================================================================================= LIMITING
249	!$acc parallel loop collapse(2) async
250	DO l =ll_begin,ll_end
251	!DIR$ SIMD
252	DO ij=ij_begin,ij_end
253	! maggrd = dot_product_3d(gradq3d(ij,l,:),gradq3d(ij,l,:))
254	maggrd = gradq3d(ij,l,1)gradq3d(ij,l,1) + gradq3d(ij,l,2)gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3)
255	maggrd = sqrt(maggrd)
256	maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij)
257	minq_c = qi(ij,l) - maggrd*sqrt_leng(ij)
258	maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
259	qi(ij+t_rdown,l),qi(ij+t_ldown,l))
260	minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
261	qi(ij+t_rdown,l),qi(ij+t_ldown,l))
262	IF ((maxq_c - qi(ij,l)) /= 0.0) THEN
263	alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
264	alphamx = max(alphamx,0.0)
265	ELSE
266	alphamx = 0.0
267	ENDIF
268	IF ((minq_c - qi(ij,l)) /= 0.0) THEN
269	alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
270	alphami = max(alphami,0.0)
271	ELSE
272	alphami = 0.0
273	ENDIF
274	alpha = min(alphamx,alphami,1.0)
275	! gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
276	gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1)
277	gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2)
278	gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3)
279	END DO
280	END DO
281
282	CALL trace_end("compute_gradq3d3")
283
284	!$acc end data
285
286	CONTAINS
287
288	SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det)
289	IMPLICIT NONE
290	INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
291	REAL(rstd), INTENT(IN) :: q(iim*jjm,llm)
292	! REAL(rstd), INTENT(OUT) :: dq(3), det
293	REAL(rstd), INTENT(OUT) :: dq1,dq2,dq3,det
294	REAL(rstd) :: dq(3)
295
296	REAL(rstd) :: A(3,3)
297
298	! TODO : replace A by A1,A2,A3
299	A(1,1)=xyz_i(n1,1)-xyz_i(n0,1); A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3)
300	A(2,1)=xyz_i(n2,1)-xyz_i(n0,1); A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3)
301	A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3)
302
303	dq(1) = q(n1,l)-q(n0,l)
304	dq(2) = q(n2,l)-q(n0,l)
305	dq(3) = 0.0
306
307	! CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
308	! CALL determinant(A(:,1),A(:,2),A(:,3),det)
309	! CALL determinant(dq,A(:,2),A(:,3),detx)
310	! CALL determinant(A(:,1),dq,A(:,3),dety)
311	! CALL determinant(A(:,1),A(:,2),dq,detz)
312	! dq(1) = detx
313	! dq(2) = dety
314	! dq(3) = detz
315
316	CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
317	CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1)
318	CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2)
319	CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3)
320
321	END SUBROUTINE gradq
322
323	!==========================================================================
324	! PURE SUBROUTINE determinant(a1,a2,a3,det)
325	! IMPLICIT NONE
326	! REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
327	! REAL(rstd), INTENT(OUT) :: det
328	! REAL(rstd) :: x1,x2,x3
329	! x1 = a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
330	! x2 = a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
331	! x3 = a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
332	! det = x1 - x2 + x3
333	! END SUBROUTINE determinant
334
335	SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det)
336	IMPLICIT NONE
337	REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33
338	REAL(rstd), INTENT(OUT) :: det
339	REAL(rstd) :: x1,x2,x3
340	x1 = a11 * (a22 * a33 - a23 * a32)
341	x2 = a12 * (a21 * a33 - a23 * a31)
342	x3 = a13 * (a21 * a32 - a22 * a31)
343	det = x1 - x2 + x3
344	END SUBROUTINE determinant
345
346
347	END SUBROUTINE compute_gradq3d
348
349	! Backward trajectories, for use with Miura approach
350	SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc, &
351	xyz_e, de, wee, le ) ! metrics terms
352	USE trace
353	USE omp_para
354	IMPLICIT NONE
355	REAL(rstd),INTENT(IN) :: normal(3iimjjm,3)
356	REAL(rstd),INTENT(IN) :: tangent(3iimjjm,3)
357	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
358	REAL(rstd),INTENT(OUT) :: cc(3iimjjm,llm,3) ! start of backward trajectory
359	REAL(rstd),INTENT(IN) :: tau
360	! metrics terms
361	REAL(rstd),INTENT(IN) :: xyz_e(iim3jjm,3)
362	REAL(rstd),INTENT(IN) :: de(iim3jjm)
363	REAL(rstd),INTENT(IN) :: wee(iim3jjm,5,2)
364	REAL(rstd),INTENT(IN) :: le(iim3jjm)
365
366	REAL(rstd) :: v_e(3), up_e
367	INTEGER :: ij,l
368
369	CALL trace_start("compute_backward_traj")
370
371	! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
372
373	!$acc data present(ue(:,:), cc(:,:,:), normal(:,:), tangent(:,:), xyz_e(:,:), de(:), wee(:,:,:), le(:)) async
374
375	! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
376	!$acc parallel loop private(up_e, v_e) collapse(2) gang vector async
377	DO l = ll_begin,ll_end
378	!DIR$ SIMD
379	DO ij=ij_begin,ij_end
380	up_e =1/de(ij+u_right)*( &
381	wee(ij+u_right,1,1)le(ij+u_rup)ue(ij+u_rup,l)+ &
382	wee(ij+u_right,2,1)le(ij+u_lup)ue(ij+u_lup,l)+ &
383	wee(ij+u_right,3,1)le(ij+u_left)ue(ij+u_left,l)+ &
384	wee(ij+u_right,4,1)le(ij+u_ldown)ue(ij+u_ldown,l)+ &
385	wee(ij+u_right,5,1)le(ij+u_rdown)ue(ij+u_rdown,l)+ &
386	wee(ij+u_right,1,2)le(ij+t_right+u_ldown)ue(ij+t_right+u_ldown,l)+ &
387	wee(ij+u_right,2,2)le(ij+t_right+u_rdown)ue(ij+t_right+u_rdown,l)+ &
388	wee(ij+u_right,3,2)le(ij+t_right+u_right)ue(ij+t_right+u_right,l)+ &
389	wee(ij+u_right,4,2)le(ij+t_right+u_rup)ue(ij+t_right+u_rup,l)+ &
390	wee(ij+u_right,5,2)le(ij+t_right+u_lup)ue(ij+t_right+u_lup,l) &
391	)
392	v_e = ue(ij+u_right,l)normal(ij+u_right,:) + up_etangent(ij+u_right,:)
393	cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
394
395	up_e=1/de(ij+u_lup)*( &
396	wee(ij+u_lup,1,1)le(ij+u_left)ue(ij+u_left,l)+ &
397	wee(ij+u_lup,2,1)le(ij+u_ldown)ue(ij+u_ldown,l)+ &
398	wee(ij+u_lup,3,1)le(ij+u_rdown)ue(ij+u_rdown,l)+ &
399	wee(ij+u_lup,4,1)le(ij+u_right)ue(ij+u_right,l)+ &
400	wee(ij+u_lup,5,1)le(ij+u_rup)ue(ij+u_rup,l)+ &
401	wee(ij+u_lup,1,2)le(ij+t_lup+u_right)ue(ij+t_lup+u_right,l)+ &
402	wee(ij+u_lup,2,2)le(ij+t_lup+u_rup)ue(ij+t_lup+u_rup,l)+ &
403	wee(ij+u_lup,3,2)le(ij+t_lup+u_lup)ue(ij+t_lup+u_lup,l)+ &
404	wee(ij+u_lup,4,2)le(ij+t_lup+u_left)ue(ij+t_lup+u_left,l)+ &
405	wee(ij+u_lup,5,2)le(ij+t_lup+u_ldown)ue(ij+t_lup+u_ldown,l) &
406	)
407	v_e = ue(ij+u_lup,l)normal(ij+u_lup,:) + up_etangent(ij+u_lup,:)
408	cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
409
410
411	up_e=1/de(ij+u_ldown)*( &
412	wee(ij+u_ldown,1,1)le(ij+u_rdown)ue(ij+u_rdown,l)+ &
413	wee(ij+u_ldown,2,1)le(ij+u_right)ue(ij+u_right,l)+ &
414	wee(ij+u_ldown,3,1)le(ij+u_rup)ue(ij+u_rup,l)+ &
415	wee(ij+u_ldown,4,1)le(ij+u_lup)ue(ij+u_lup,l)+ &
416	wee(ij+u_ldown,5,1)le(ij+u_left)ue(ij+u_left,l)+ &
417	wee(ij+u_ldown,1,2)le(ij+t_ldown+u_lup)ue(ij+t_ldown+u_lup,l)+ &
418	wee(ij+u_ldown,2,2)le(ij+t_ldown+u_left)ue(ij+t_ldown+u_left,l)+ &
419	wee(ij+u_ldown,3,2)le(ij+t_ldown+u_ldown)ue(ij+t_ldown+u_ldown,l)+ &
420	wee(ij+u_ldown,4,2)le(ij+t_ldown+u_rdown)ue(ij+t_ldown+u_rdown,l)+ &
421	wee(ij+u_ldown,5,2)le(ij+t_ldown+u_right)ue(ij+t_ldown+u_right,l) &
422	)
423	v_e = ue(ij+u_ldown,l)normal(ij+u_ldown,:) + up_etangent(ij+u_ldown,:)
424	cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
425	ENDDO
426	END DO
427	!$acc end data
428	CALL trace_end("compute_backward_traj")
429
430	END SUBROUTINE compute_backward_traj
431
432	! Horizontal transport (S. Dubey, T. Dubos)
433	! Slope-limited van Leer approach with hexagons
434	SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass, qi, qfluxt, &
435	Ai, xyz_i) ! metrics terms
436	USE trace
437	USE omp_para
438	USE abort_mod
439	IMPLICIT NONE
440	LOGICAL, INTENT(IN) :: update_mass, diagflux_on
441	REAL(rstd), INTENT(IN) :: gradq3d(iim*jjm,llm,3)
442	REAL(rstd), INTENT(IN) :: hfluxt(3iimjjm,llm) ! mass flux
443	REAL(rstd), INTENT(IN) :: cc(3iimjjm,llm,3) ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
444	REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
445	REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
446	REAL(rstd), INTENT(INOUT) :: qfluxt(3iimjjm,MERGE(llm,1,diagflux_on)) ! time-integrated tracer flux
447	! metrics terms
448	REAL(rstd), INTENT(IN) :: Ai(iim*jjm)
449	REAL(rstd), INTENT(IN) :: xyz_i(iim*jjm,3)
450
451	REAL(rstd) :: dq,dmass,qe,newmass
452	REAL(rstd) :: qflux(3iimjjm,llm)
453	INTEGER :: ij,l,ij_tmp
454
455	IF(diagflux_on) CALL abort_acc("compute_advect_horiz : diagflux_on")
456
457	CALL trace_start("compute_advect_horiz")
458	#include "../kernels/advect_horiz.k90"
459	CALL trace_end("compute_advect_horiz")
460
461	END SUBROUTINE compute_advect_horiz
462
463
464	END MODULE advect_mod

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

Context Navigation

source: codes/icosagcm/trunk/src/transport/advect.F90 @ 954

Download in other formats: