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

Last change on this file since 150 was 148, checked in by ymipsl, 11 years ago

Various optimisations

  • dissipation is not called every timestep (similar way than LMDZ)
  • transfert size of halos has been reduced : need to synchronise redondant data between tiles at itau_sync timestep

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