source: codes/icosagcm/trunk/src/transport/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
RevLine 
[22]1MODULE advect_mod
2  USE icosa
3  IMPLICIT NONE
[17]4
[138]5  PRIVATE
6 
7  PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz
8
[17]9CONTAINS
10
[22]11  !==========================================================================
[17]12
[252]13  SUBROUTINE init_advect(normal,tangent,sqrt_leng)
[22]14    IMPLICIT NONE
15    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
16    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
[252]17    REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm)
[17]18
[174]19    INTEGER :: ij
[22]20
[487]21!DIR$ SIMD
[174]22      DO ij=ij_begin,ij_end
[22]23
[174]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)
[22]26
[174]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)
[22]29
[174]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)
[22]32
[174]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)
[22]35
[174]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)
[22]38
[174]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)
[148]41
[252]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)) )
[22]45    ENDDO
46
[17]47  END SUBROUTINE init_advect
48
[22]49  !=======================================================================================
[17]50
[954]51  SUBROUTINE compute_gradq3d(qi,sqrt_leng,gradq3d,xyz_i,xyz_v)
[148]52    USE trace
[151]53    USE omp_para
[22]54    IMPLICIT NONE
[896]55    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
56    REAL(rstd),INTENT(IN)  :: sqrt_leng(iim*jjm)
[186]57    REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
58    REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
[896]59    REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3)
[17]60    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
[899]61    REAL(rstd) :: alphamx,alphami,alpha,maggrd
[17]62    REAL(rstd) :: arr(2*iim*jjm)
[148]63    REAL(rstd) :: ar(iim*jjm)
[17]64    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
[899]65    INTEGER :: ij,k,l
[186]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)
[953]70   
[186]71    CALL trace_start("compute_gradq3d1")
[148]72
[138]73    ! TODO : precompute ar, drop arr as output argument of gradq ?
74
[22]75    !========================================================================================== GRADIENT
[138]76    ! Compute gradient at triangles solving a linear system
77    ! arr = area of triangle joining centroids of hexagons
[186]78!     DO l = ll_begin,ll_end
[487]79!!DIR$ SIMD
[186]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
[953]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
[487]92!DIR$ SIMD
[174]93      DO ij=ij_begin_ext,ij_end_ext
[186]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
[953]156     ENDDO
[186]157     
[953]158     !$acc parallel loop collapse(2) async private(A, dq)
159     DO l = ll_begin,ll_end
[186]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
[22]224    END DO
[17]225
[487]226!DIR$ SIMD
[953]227    !$acc parallel loop async
[186]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
[148]230    ENDDO
[186]231    CALL trace_end("compute_gradq3d1")
232    CALL trace_start2("compute_gradq3d2")
[148]233     
[953]234    !$acc parallel loop collapse(3) async
[148]235    DO k=1,3
[151]236      DO l =ll_begin,ll_end
[487]237!DIR$ SIMD
[174]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) 
[148]242        END DO
243      END DO
244    ENDDO
[186]245    CALL trace_end2("compute_gradq3d2")
246    CALL trace_start("compute_gradq3d3")
[22]247
248    !============================================================================================= LIMITING
[953]249    !$acc parallel loop collapse(2) async
[151]250    DO l =ll_begin,ll_end
[487]251!DIR$ SIMD
[174]252      DO ij=ij_begin,ij_end
[899]253!             maggrd =  dot_product_3d(gradq3d(ij,l,:),gradq3d(ij,l,:))
[186]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) 
[22]255             maggrd = sqrt(maggrd) 
[252]256             maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij)
257             minq_c = qi(ij,l) - maggrd*sqrt_leng(ij)
[174]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))
[953]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
[22]274             alpha   = min(alphamx,alphami,1.0)
[186]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)
[22]279       END DO
[17]280    END DO
[148]281
[186]282  CALL trace_end("compute_gradq3d3")
[953]283
284  !$acc end data
[186]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
[22]347  END SUBROUTINE compute_gradq3d
348
[138]349  ! Backward trajectories, for use with Miura approach
[954]350  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc, &
351                                   xyz_e, de, wee, le         ) ! metrics terms
[148]352    USE trace
[151]353    USE omp_para
[22]354    IMPLICIT NONE
[137]355    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
356    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
357    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
358    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory
359    REAL(rstd),INTENT(IN)    :: tau
[954]360! metrics terms
361    REAL(rstd),INTENT(IN)    :: xyz_e(iim*3*jjm,3)
362    REAL(rstd),INTENT(IN)    :: de(iim*3*jjm)
363    REAL(rstd),INTENT(IN)    :: wee(iim*3*jjm,5,2)
364    REAL(rstd),INTENT(IN)    :: le(iim*3*jjm)
[137]365
[899]366    REAL(rstd) :: v_e(3), up_e   
[174]367    INTEGER :: ij,l
[137]368
[148]369    CALL trace_start("compute_backward_traj")
370
[138]371    ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
[953]372
373    !$acc data present(ue(:,:), cc(:,:,:), normal(:,:), tangent(:,:), xyz_e(:,:), de(:), wee(:,:,:), le(:)) async
374
[137]375    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
[953]376    !$acc parallel loop private(up_e, v_e) collapse(2) gang vector async
[151]377    DO l = ll_begin,ll_end
[487]378!DIR$ SIMD
[174]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)         &
[137]391                  )
[174]392             v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:)
393             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
[137]394
[174]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)           &
[137]406                  )
[174]407             v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:)
408             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
409           
[137]410
[174]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)     &
[137]422                  )
[174]423             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:)
424             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
[137]425       ENDDO
426    END DO
[953]427    !$acc end data
[148]428    CALL trace_end("compute_backward_traj")
429
[137]430  END SUBROUTINE compute_backward_traj
431
432  ! Horizontal transport (S. Dubey, T. Dubos)
433  ! Slope-limited van Leer approach with hexagons
[954]434  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass, qi, qfluxt,  &
435                                  Ai, xyz_i)   ! metrics terms
[148]436    USE trace
[151]437    USE omp_para
[953]438    USE abort_mod
[137]439    IMPLICIT NONE
[599]440    LOGICAL, INTENT(IN)       :: update_mass, diagflux_on
[138]441    REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
442    REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
443    REAL(rstd), INTENT(IN)    :: cc(3*iim*jjm,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)
[894]446    REAL(rstd), INTENT(INOUT) :: qfluxt(3*iim*jjm,MERGE(llm,1,diagflux_on)) ! time-integrated tracer flux
[953]447! metrics terms
[954]448    REAL(rstd), INTENT(IN)    :: Ai(iim*jjm)
449    REAL(rstd), INTENT(IN)    :: xyz_i(iim*jjm,3)
[953]450   
451    REAL(rstd) :: dq,dmass,qe,newmass
[138]452    REAL(rstd) :: qflux(3*iim*jjm,llm)
[953]453    INTEGER :: ij,l,ij_tmp
[17]454
[953]455    IF(diagflux_on) CALL abort_acc("compute_advect_horiz : diagflux_on")
456
[148]457    CALL trace_start("compute_advect_horiz")
[599]458#include "../kernels/advect_horiz.k90"
[148]459    CALL trace_end("compute_advect_horiz")
[953]460
[22]461  END SUBROUTINE compute_advect_horiz
[138]462
[17]463
[22]464END MODULE advect_mod
Note: See TracBrowser for help on using the repository browser.