source: codes/icosagcm/trunk/src/advect.f90 @ 238

Last change on this file since 238 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
File size: 21.9 KB
Line 
1MODULE 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
9CONTAINS
10
11  !==========================================================================
12
13  SUBROUTINE init_advect(normal,tangent,one_over_sqrt_leng)
14    IMPLICIT NONE
15    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
16    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
17    REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)
18
19    INTEGER :: ij
20
21!$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          one_over_sqrt_leng(ij) = 1./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_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
52    USE trace
53    USE omp_para
54    IMPLICIT NONE
55    REAL(rstd),INTENT(IN)  :: qi_in(iim*jjm,llm)
56    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng_in(iim*jjm)
57    REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
58    REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
59    REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3) 
60    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
61    REAL(rstd) :: alphamx,alphami,alpha ,maggrd
62    REAL(rstd) :: leng1,leng2
63    REAL(rstd) :: arr(2*iim*jjm)
64    REAL(rstd) :: ar(iim*jjm)
65    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
66    INTEGER :: ij,k,ind,l
67    REAL(rstd)  :: qi(iim*jjm,llm)
68    REAL(rstd)  :: one_over_sqrt_leng(iim*jjm)
69    REAL(rstd) :: gradq3d(iim*jjm,llm,3) 
70    REAL(rstd) :: detx,dety,detz,det
71    REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33
72    REAL(rstd) :: x1,x2,x3
73    REAL(rstd) :: dq(3)
74
75    qi=qi_in
76    one_over_sqrt_leng=one_over_sqrt_leng_in
77   
78    CALL trace_start("compute_gradq3d1")
79
80    ! TODO : precompute ar, drop arr as output argument of gradq ?
81
82    !========================================================================================== GRADIENT
83    ! Compute gradient at triangles solving a linear system
84    ! arr = area of triangle joining centroids of hexagons
85!     DO l = ll_begin,ll_end
86!!$SIMD
87!      DO ij=ij_begin_ext,ij_end_ext
88!!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
89!!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
90!             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))
91!             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))
92!       END DO
93!    END DO
94
95     DO l = ll_begin,ll_end 
96!$SIMD
97      DO ij=ij_begin_ext,ij_end_ext
98!       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))
99
100       
101        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) 
102        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) 
103        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)
104   
105        dq(1) = qi(ij+t_rup,l)-qi(ij,l)
106        dq(2) = qi(ij+t_lup,l)-qi(ij,l)
107        dq(3) = 0.0
108
109
110!        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)
111
112         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
113         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
114         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
115
116         x1 =  a11 * (a22 * a33 - a23 * a32)
117         x2 =  a12 * (a21 * a33 - a23 * a31)
118         x3 =  a13 * (a21 * a32 - a22 * a31)
119         det =  x1 - x2 + x3
120                 
121!        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)
122
123         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
124         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
125         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
126
127         x1 =  a11 * (a22 * a33 - a23 * a32)
128         x2 =  a12 * (a21 * a33 - a23 * a31)
129         x3 =  a13 * (a21 * a32 - a22 * a31)
130         detx =  x1 - x2 + x3
131       
132!        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)
133
134         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
135         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
136         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
137
138         x1 =  a11 * (a22 * a33 - a23 * a32)
139         x2 =  a12 * (a21 * a33 - a23 * a31)
140         x3 =  a13 * (a21 * a32 - a22 * a31)
141         dety =  x1 - x2 + x3
142
143!        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)
144
145         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
146         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
147         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
148
149         x1 =  a11 * (a22 * a33 - a23 * a32)
150         x2 =  a12 * (a21 * a33 - a23 * a31)
151         x3 =  a13 * (a21 * a32 - a22 * a31)
152         detz =  x1 - x2 + x3
153
154        gradtri(ij+z_up,l,1) = detx
155        gradtri(ij+z_up,l,2) = dety
156        gradtri(ij+z_up,l,3) = detz
157        arr(ij+z_up) = det
158       
159      ENDDO
160     
161      DO ij=ij_begin_ext,ij_end_ext
162
163
164!        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))
165
166        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) 
167        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) 
168        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)
169 
170        dq(1) = qi(ij+t_ldown,l)-qi(ij,l)
171        dq(2) = qi(ij+t_rdown,l)-qi(ij,l)
172        dq(3) = 0.0
173
174
175!        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)
176
177         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
178         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
179         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
180
181         x1 =  a11 * (a22 * a33 - a23 * a32)
182         x2 =  a12 * (a21 * a33 - a23 * a31)
183         x3 =  a13 * (a21 * a32 - a22 * a31)
184         det =  x1 - x2 + x3
185                 
186!        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)
187
188         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
189         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
190         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
191
192         x1 =  a11 * (a22 * a33 - a23 * a32)
193         x2 =  a12 * (a21 * a33 - a23 * a31)
194         x3 =  a13 * (a21 * a32 - a22 * a31)
195         detx =  x1 - x2 + x3
196       
197!        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)
198
199         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
200         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
201         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
202
203         x1 =  a11 * (a22 * a33 - a23 * a32)
204         x2 =  a12 * (a21 * a33 - a23 * a31)
205         x3 =  a13 * (a21 * a32 - a22 * a31)
206         dety =  x1 - x2 + x3
207
208!        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)
209
210         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
211         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
212         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
213
214         x1 =  a11 * (a22 * a33 - a23 * a32)
215         x2 =  a12 * (a21 * a33 - a23 * a31)
216         x3 =  a13 * (a21 * a32 - a22 * a31)
217         detz =  x1 - x2 + x3
218
219         gradtri(ij+z_down,l,1) = detx
220         gradtri(ij+z_down,l,2) = dety
221         gradtri(ij+z_down,l,3) = detz
222         arr(ij+z_down) = det
223
224      END DO
225    END DO
226
227!$SIMD
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    DO k=1,3
235      DO l =ll_begin,ll_end
236!$SIMD
237      DO ij=ij_begin,ij_end
238              gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k)   +   &
239                                 gradtri(ij+z_rup,l,k) +  gradtri(ij+z_ldown,l,k)   +   & 
240                                 gradtri(ij+z_lup,l,k)+    gradtri(ij+z_rdown,l,k) ) / ar(ij) 
241        END DO
242      END DO
243    ENDDO
244    CALL trace_end2("compute_gradq3d2")
245    CALL trace_start("compute_gradq3d3")
246
247    !============================================================================================= LIMITING
248    DO l =ll_begin,ll_end
249!$SIMD
250      DO ij=ij_begin,ij_end
251!             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
252             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) 
253             maggrd = sqrt(maggrd) 
254             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
255             minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij)
256             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), &
257                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
258             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), &
259                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
260             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
261             alphamx = max(alphamx,0.0)
262             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
263             alphami = max(alphami,0.0) 
264             alpha   = min(alphamx,alphami,1.0)
265!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
266             gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1)
267             gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2)
268             gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3)
269       END DO
270    END DO
271
272  CALL trace_end("compute_gradq3d3")
273 
274  gradq3d_out=gradq3d
275 
276  CONTAINS
277
278    SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det)
279    IMPLICIT NONE
280    INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
281    REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm)
282!    REAL(rstd), INTENT(OUT)    :: dq(3), det
283    REAL(rstd), INTENT(OUT)    :: dq1,dq2,dq3,det
284    REAL(rstd)    :: dq(3)
285   
286    REAL(rstd)  ::detx,dety,detz
287    INTEGER    :: info
288    INTEGER    :: IPIV(3)
289
290    REAL(rstd) :: A(3,3)
291    REAL(rstd) :: B(3)
292   
293    ! TODO : replace A by A1,A2,A3
294    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) 
295    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) 
296    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3)
297
298    dq(1) = q(n1,l)-q(n0,l)
299    dq(2) = q(n2,l)-q(n0,l)
300    dq(3) = 0.0
301
302    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
303!    CALL determinant(A(:,1),A(:,2),A(:,3),det)
304!    CALL determinant(dq,A(:,2),A(:,3),detx)
305!    CALL determinant(A(:,1),dq,A(:,3),dety)
306!    CALL determinant(A(:,1),A(:,2),dq,detz)
307!    dq(1) = detx
308!    dq(2) = dety
309!    dq(3) = detz
310
311    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)
312    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)
313    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)
314    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)
315
316  END SUBROUTINE gradq
317
318  !==========================================================================
319!  PURE SUBROUTINE determinant(a1,a2,a3,det)
320!    IMPLICIT NONE
321!    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
322!    REAL(rstd), INTENT(OUT) :: det
323!    REAL(rstd) ::  x1,x2,x3
324!    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
325!    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
326!    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
327!    det =  x1 - x2 + x3
328!  END SUBROUTINE determinant
329
330   SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det)
331    IMPLICIT NONE
332    REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33
333    REAL(rstd), INTENT(OUT) :: det
334    REAL(rstd) ::  x1,x2,x3
335    x1 =  a11 * (a22 * a33 - a23 * a32)
336    x2 =  a12 * (a21 * a33 - a23 * a31)
337    x3 =  a13 * (a21 * a32 - a22 * a31)
338    det =  x1 - x2 + x3
339  END SUBROUTINE determinant
340
341
342  END SUBROUTINE compute_gradq3d
343
344  ! Backward trajectories, for use with Miura approach
345  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
346    USE trace
347    USE omp_para
348    IMPLICIT NONE
349    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
350    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
351    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
352    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory
353    REAL(rstd),INTENT(IN)    :: tau
354
355    REAL(rstd) :: v_e(3), up_e, qe, ed(3)   
356    INTEGER :: ij,l
357
358    CALL trace_start("compute_backward_traj")
359
360    ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
361   
362    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
363    DO l = ll_begin,ll_end
364!$SIMD
365      DO ij=ij_begin,ij_end
366             up_e =1/de(ij+u_right)*(                                                   &
367                  wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
368                  wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
369                  wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+                      &
370                  wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                    &
371                  wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    & 
372                  wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+    &
373                  wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+    &
374                  wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+    &
375                  wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+        &
376                  wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+t_right+u_lup,l)         &
377                  )
378             v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:)
379             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
380
381             up_e=1/de(ij+u_lup)*(                                                      &
382                  wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+                        &
383                  wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                      &
384                  wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                      &
385                  wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+                      &
386                  wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                          & 
387                  wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+          &
388                  wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+              &
389                  wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+              &
390                  wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+            &
391                  wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+t_lup+u_ldown,l)           &
392                  )
393             v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:)
394             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
395           
396
397             up_e=1/de(ij+u_ldown)*(                                                    &
398                  wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    &
399                  wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+                    &
400                  wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
401                  wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
402                  wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+                      & 
403                  wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+        &
404                  wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+      &
405                  wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+    &
406                  wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+    &
407                  wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+t_ldown+u_right,l)     &
408                  )
409             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:)
410             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
411       ENDDO
412    END DO
413
414    CALL trace_end("compute_backward_traj")
415
416  END SUBROUTINE compute_backward_traj
417
418  ! Horizontal transport (S. Dubey, T. Dubos)
419  ! Slope-limited van Leer approach with hexagons
420  SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
421    USE trace
422    USE omp_para
423    IMPLICIT NONE
424    LOGICAL, INTENT(IN)       :: update_mass
425    REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
426    REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
427    REAL(rstd), INTENT(IN)    :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
428    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
429    REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
430
431    REAL(rstd) :: dq,dmass,qe,ed(3), newmass
432    REAL(rstd) :: qflux(3*iim*jjm,llm)
433    INTEGER :: ij,k,l
434
435    CALL trace_start("compute_advect_horiz")
436
437    ! evaluate tracer field at point cc using piecewise linear reconstruction
438    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
439    ! ne*hfluxt>0 iff outgoing
440    DO l = ll_begin,ll_end
441
442!$SIMD
443      DO ij=ij_begin_ext,ij_end_ext
444
445             IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN
446                ed = cc(ij+u_right,l,:) - xyz_i(ij,:)
447!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
448                 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
449             ELSE
450                ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:)
451!                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed)
452                qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3) 
453             ENDIF
454             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
455             
456             IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN
457                ed = cc(ij+u_lup,l,:) - xyz_i(ij,:)
458!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
459                qe = qi(ij,l)  + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
460             ELSE
461                ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:)
462!                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed)
463                qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3) 
464             ENDIF
465             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 
466
467             IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN
468                ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:)
469!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
470                qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
471             ELSE
472                ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:)
473!                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed)
474                qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3) 
475             ENDIF
476             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
477       END DO
478    END DO
479
480    ! update q and, if update_mass, update mass
481    DO l = ll_begin,ll_end
482!$SIMD
483      DO ij=ij_begin,ij_end
484             ! sign convention as Ringler et al. (2010) eq. 21
485             dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   &
486                     + hfluxt(ij+u_lup,l)  *ne(ij,lup)     &
487                     + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   &
488                     + hfluxt(ij+u_rup,l)  *ne(ij,rup)     &
489                     + hfluxt(ij+u_left,l) *ne(ij,left)    &
490                     + hfluxt(ij+u_rdown,l)*ne(ij,rdown)
491
492             dq    =   qflux(ij+u_right,l) *ne(ij,right)   &
493                     + qflux(ij+u_lup,l)   *ne(ij,lup)     &
494                     + qflux(ij+u_ldown,l) *ne(ij,ldown)   &
495                     + qflux(ij+u_rup,l)   *ne(ij,rup)     &
496                     + qflux(ij+u_left,l)  *ne(ij,left)    &
497                     + qflux(ij+u_rdown,l) *ne(ij,rdown)
498
499             
500             newmass = mass(ij,l) - dmass/Ai(ij)
501             qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass
502             IF(update_mass) mass(ij,l) = newmass
503
504       END DO
505    END DO
506
507    CALL trace_end("compute_advect_horiz")
508  CONTAINS
509
510    !====================================================================================
511    PURE REAL(rstd) FUNCTION sum2(a1,a2)
512      IMPLICIT NONE
513      REAL(rstd),INTENT(IN):: a1(3), a2(3)
514      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
515      ! sum2 = 0. ! Godunov scheme
516    END FUNCTION sum2
517
518  END SUBROUTINE compute_advect_horiz
519
520
521END MODULE advect_mod
Note: See TracBrowser for help on using the repository browser.