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

Last change on this file since 176 was 174, checked in by ymipsl, 11 years ago

Transform 2 loops on i and j in one loop ij for efficiency, vectorization and future GPU programing

YM
File size: 14.5 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,one_over_sqrt_leng,gradq3d)
52    USE trace
53    USE omp_para
54    IMPLICIT NONE
55    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
56    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng(iim*jjm)
57    REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3) 
58    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
59    REAL(rstd) :: alphamx,alphami,alpha ,maggrd
60    REAL(rstd) :: leng1,leng2
61    REAL(rstd) :: arr(2*iim*jjm)
62    REAL(rstd) :: ar(iim*jjm)
63    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
64    INTEGER :: ij,k,ind,l
65
66    CALL trace_start("compute_gradq3d")
67
68    ! TODO : precompute ar, drop arr as output argument of gradq ?
69
70    !========================================================================================== GRADIENT
71    ! Compute gradient at triangles solving a linear system
72    ! arr = area of triangle joining centroids of hexagons
73     DO l = ll_begin,ll_end 
74!$SIMD
75      DO ij=ij_begin_ext,ij_end_ext
76             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
77             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
78       END DO
79    END DO
80
81!$SIMD
82      DO ij=ij_begin,ij_end
83         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
84    ENDDO
85     
86    DO k=1,3
87      DO l =ll_begin,ll_end
88!$SIMD
89      DO ij=ij_begin,ij_end
90              gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k)   +   &
91                                 gradtri(ij+z_rup,l,k) +  gradtri(ij+z_ldown,l,k)   +   & 
92                                 gradtri(ij+z_lup,l,k)+    gradtri(ij+z_rdown,l,k) ) / ar(ij) 
93        END DO
94      END DO
95    ENDDO
96
97    !============================================================================================= LIMITING
98    DO l =ll_begin,ll_end
99!$SIMD
100      DO ij=ij_begin,ij_end
101             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
102             maggrd = sqrt(maggrd) 
103             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
104             minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij)
105             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), &
106                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
107             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), &
108                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
109             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
110             alphamx = max(alphamx,0.0)
111             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
112             alphami = max(alphami,0.0) 
113             alpha   = min(alphamx,alphami,1.0)
114             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
115       END DO
116    END DO
117
118  CALL trace_end("compute_gradq3d")
119  END SUBROUTINE compute_gradq3d
120
121  ! Backward trajectories, for use with Miura approach
122  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
123    USE trace
124    USE omp_para
125    IMPLICIT NONE
126    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
127    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
128    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
129    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory
130    REAL(rstd),INTENT(IN)    :: tau
131
132    REAL(rstd) :: v_e(3), up_e, qe, ed(3)   
133    INTEGER :: ij,l
134
135    CALL trace_start("compute_backward_traj")
136
137    ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
138   
139    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
140    DO l = ll_begin,ll_end
141!$SIMD
142      DO ij=ij_begin,ij_end
143             up_e =1/de(ij+u_right)*(                                                   &
144                  wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
145                  wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
146                  wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+                      &
147                  wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                    &
148                  wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    & 
149                  wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+    &
150                  wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+    &
151                  wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+    &
152                  wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+        &
153                  wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+t_right+u_lup,l)         &
154                  )
155             v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:)
156             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
157
158             up_e=1/de(ij+u_lup)*(                                                      &
159                  wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+                        &
160                  wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                      &
161                  wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                      &
162                  wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+                      &
163                  wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                          & 
164                  wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+          &
165                  wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+              &
166                  wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+              &
167                  wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+            &
168                  wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+t_lup+u_ldown,l)           &
169                  )
170             v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:)
171             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
172           
173
174             up_e=1/de(ij+u_ldown)*(                                                    &
175                  wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    &
176                  wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+                    &
177                  wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
178                  wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
179                  wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+                      & 
180                  wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+        &
181                  wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+      &
182                  wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+    &
183                  wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+    &
184                  wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+t_ldown+u_right,l)     &
185                  )
186             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:)
187             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
188       ENDDO
189    END DO
190
191    CALL trace_end("compute_backward_traj")
192
193  END SUBROUTINE compute_backward_traj
194
195  ! Horizontal transport (S. Dubey, T. Dubos)
196  ! Slope-limited van Leer approach with hexagons
197  SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
198    USE trace
199    USE omp_para
200    IMPLICIT NONE
201    LOGICAL, INTENT(IN)       :: update_mass
202    REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
203    REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
204    REAL(rstd), INTENT(IN)    :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
205    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
206    REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
207
208    REAL(rstd) :: dq,dmass,qe,ed(3), newmass
209    REAL(rstd) :: qflux(3*iim*jjm,llm)
210    INTEGER :: ij,k,l
211
212    CALL trace_start("compute_advect_horiz")
213
214    ! evaluate tracer field at point cc using piecewise linear reconstruction
215    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
216    ! ne*hfluxt>0 iff outgoing
217    DO l = ll_begin,ll_end
218
219!$SIMD
220      DO ij=ij_begin_ext,ij_end_ext
221
222             IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN
223                ed = cc(ij+u_right,l,:) - xyz_i(ij,:)
224                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
225             ELSE
226                ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:)
227                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 
228             ENDIF
229             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
230             
231             IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN
232                ed = cc(ij+u_lup,l,:) - xyz_i(ij,:)
233                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
234             ELSE
235                ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:)
236                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 
237             ENDIF
238             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 
239
240             IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN
241                ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:)
242                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
243             ELSE
244                ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:)
245                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 
246             ENDIF
247             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
248       END DO
249    END DO
250
251    ! update q and, if update_mass, update mass
252    DO l = ll_begin,ll_end
253!$SIMD
254      DO ij=ij_begin,ij_end
255             ! sign convention as Ringler et al. (2010) eq. 21
256             dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   &
257                     + hfluxt(ij+u_lup,l)  *ne(ij,lup)     &
258                     + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   &
259                     + hfluxt(ij+u_rup,l)  *ne(ij,rup)     &
260                     + hfluxt(ij+u_left,l) *ne(ij,left)    &
261                     + hfluxt(ij+u_rdown,l)*ne(ij,rdown)
262
263             dq    =   qflux(ij+u_right,l) *ne(ij,right)   &
264                     + qflux(ij+u_lup,l)   *ne(ij,lup)     &
265                     + qflux(ij+u_ldown,l) *ne(ij,ldown)   &
266                     + qflux(ij+u_rup,l)   *ne(ij,rup)     &
267                     + qflux(ij+u_left,l)  *ne(ij,left)    &
268                     + qflux(ij+u_rdown,l) *ne(ij,rdown)
269
270             
271             newmass = mass(ij,l) - dmass/Ai(ij)
272             qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass
273             IF(update_mass) mass(ij,l) = newmass
274
275       END DO
276    END DO
277
278    CALL trace_end("compute_advect_horiz")
279  CONTAINS
280
281    !====================================================================================
282    PURE REAL(rstd) FUNCTION sum2(a1,a2)
283      IMPLICIT NONE
284      REAL(rstd),INTENT(IN):: a1(3), a2(3)
285      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
286      ! sum2 = 0. ! Godunov scheme
287    END FUNCTION sum2
288
289  END SUBROUTINE compute_advect_horiz
290
291  !==========================================================================
292  PURE SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq,det)
293    IMPLICIT NONE
294    INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
295    REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm)
296    REAL(rstd), INTENT(OUT)    :: dq(3), det
297   
298    REAL(rstd)  ::detx,dety,detz
299    INTEGER    :: info
300    INTEGER    :: IPIV(3)
301
302    REAL(rstd) :: A(3,3)
303    REAL(rstd) :: B(3)
304
305    ! TODO : replace A by A1,A2,A3
306    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) 
307    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) 
308    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3)
309
310    dq(1) = q(n1,l)-q(n0,l)
311    dq(2) = q(n2,l)-q(n0,l)
312    dq(3) = 0.0
313    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
314    CALL determinant(A(:,1),A(:,2),A(:,3),det)
315    CALL determinant(dq,A(:,2),A(:,3),detx)
316    CALL determinant(A(:,1),dq,A(:,3),dety)
317    CALL determinant(A(:,1),A(:,2),dq,detz)
318    dq(1) = detx
319    dq(2) = dety
320    dq(3) = detz 
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
335END MODULE advect_mod
Note: See TracBrowser for help on using the repository browser.