Context Navigation

source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/DYNAMICO/ppsrc/transport/advect.f90 @ 6612

Last change on this file since 6612 was 6612, checked in by acosce, 10 months ago
DYNAMICO used for ICOLMDZORINCA_CO2_Transport_GMD_2023
File size: 22.8 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	!--------------------------------------------------------------------------
459	!---------------------------- advect_horiz ----------------------------------
460
461	!$acc data create(qflux(:,:)) present(qi(:,:), cc(:,:,:), gradq3d(:,:,:), mass(:,:), hfluxt(:,:), Ai(:), xyz_i(:,:)) async
462
463	! evaluate tracer field at point cc using piecewise linear reconstruction
464	! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
465	! sign*hfluxt>0 iff outgoing
466	!$acc parallel loop collapse(2) async
467	DO l = ll_begin, ll_end
468	!DIR$ SIMD
469	DO ij=ij_begin_ext, ij_end_ext
470	IF(ne_right*hfluxt(ij+u_right,l)>0.) THEN
471	ij_tmp = ij
472	ELSE
473	ij_tmp = ij+t_right
474	END IF
475
476	qe = qi(ij_tmp,l)
477	qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
478	qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
479	qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
480
481	qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
482
483	IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
484	ij_tmp = ij
485	ELSE
486	ij_tmp = ij+t_lup
487	END IF
488
489	qe = qi(ij_tmp,l)
490	qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
491	qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
492	qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
493
494	qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
495
496	IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
497	ij_tmp = ij
498	ELSE
499	ij_tmp = ij+t_ldown
500	END IF
501
502	qe = qi(ij_tmp,l)
503	qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
504	qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
505	qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
506
507	qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
508	END DO
509	END DO
510
511	IF(diagflux_on) THEN
512	!$acc parallel loop collapse(2) copy(qfluxt(:,:)) async
513	DO l = ll_begin, ll_end
514	!DIR$ SIMD
515	DO ij=ij_begin_ext, ij_end_ext
516	qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
517	qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
518	qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
519	END DO
520	END DO
521	END IF
522
523	! update q and, if update_mass, update mass
524	!$acc parallel loop collapse(2) async
525	DO l = ll_begin, ll_end
526	!DIR$ SIMD
527	DO ij=ij_begin, ij_end
528	dmass=0.
529	dq=0.
530	dmass = dmass + ne_rup*hfluxt(ij+u_rup,l)
531	dq = dq + ne_rup*qflux(ij+u_rup,l)
532	dmass = dmass + ne_lup*hfluxt(ij+u_lup,l)
533	dq = dq + ne_lup*qflux(ij+u_lup,l)
534	dmass = dmass + ne_left*hfluxt(ij+u_left,l)
535	dq = dq + ne_left*qflux(ij+u_left,l)
536	dmass = dmass + ne_ldown*hfluxt(ij+u_ldown,l)
537	dq = dq + ne_ldown*qflux(ij+u_ldown,l)
538	dmass = dmass + ne_rdown*hfluxt(ij+u_rdown,l)
539	dq = dq + ne_rdown*qflux(ij+u_rdown,l)
540	dmass = dmass + ne_right*hfluxt(ij+u_right,l)
541	dq = dq + ne_right*qflux(ij+u_right,l)
542	newmass = mass(ij,l) - dmass/Ai(ij)
543	qi(ij,l) = (qi(ij,l)*mass(ij,l)-dq/Ai(ij)) / newmass
544	IF(update_mass) mass(ij,l)=newmass
545	END DO
546	END DO
547
548	!$acc end data
549
550	!---------------------------- advect_horiz ----------------------------------
551	!--------------------------------------------------------------------------
552	CALL trace_end("compute_advect_horiz")
553
554	END SUBROUTINE compute_advect_horiz
555
556
557	END MODULE advect_mod

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

Download in other formats: