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

Last change on this file since 146 was 141, checked in by dubos, 11 years ago

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