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

Last change on this file since 78 was 69, checked in by dubos, 12 years ago

Bugfix in advect.f90 (halo)
File size: 13.1 KB
Line 
1MODULE advect_mod
2  USE icosa
3  IMPLICIT NONE
4
5CONTAINS
6
7  !==========================================================================
8
9  SUBROUTINE init_advect(normal,tangent)
10    USE domain_mod
11    USE dimensions
12    USE geometry
13    USE metric
14    USE vector
15    IMPLICIT NONE
16    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
17    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
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       END DO
43    ENDDO
44
45  END SUBROUTINE init_advect
46
47  !=======================================================================================
48
49  SUBROUTINE compute_gradq3d(qi,gradq3d)
50    USE domain_mod
51    USE dimensions
52    USE geometry
53    USE metric
54    IMPLICIT NONE
55    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
56    REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3) 
57    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
58    REAL(rstd) :: alphamx,alphami,alpha ,maggrd,leng
59    REAL(rstd) :: leng1,leng2
60    REAL(rstd) :: arr(2*iim*jjm)
61    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
62    REAL(rstd) :: ar
63    INTEGER :: i,j,n,k,ind,l
64    !========================================================================================== GRADIENT
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  !===================================================================================================   
117  SUBROUTINE compute_advect_horiz(normal,tangent,qi,gradq3d,him,ue,he,bigt)
118    USE domain_mod
119    USE dimensions
120    USE geometry
121    USE metric
122    IMPLICIT NONE
123    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
124    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
125    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm)
126    REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
127    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm)
128    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
129    REAL(rstd),INTENT(IN)    :: he(3*iim*jjm,llm) ! mass flux
130    REAL(rstd),INTENT(IN)    :: bigt
131
132    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 
133    REAL(rstd) :: cc(3*iim*jjm,llm,3) 
134    REAL(rstd) :: v_e(3*iim*jjm,llm,3)
135    REAL(rstd) :: up_e
136
137    REAL(rstd) :: qe(3*iim*jjm,llm)
138    REAL(rstd) :: ed(3)   
139    REAL(rstd) :: flxx(3*iim*jjm,llm)
140    INTEGER :: i,j,n,k,l
141    REAL(rstd):: him_old(llm) 
142
143
144    !go to  123
145    DO l = 1,llm
146       DO j=jj_begin-1,jj_end+1
147          DO i=ii_begin-1,ii_end+1
148             n=(j-1)*iim+i
149
150             up_e =1/de(n+u_right)*(                                                   &
151                  wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
152                  wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
153                  wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      &
154                  wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    &
155                  wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
156                  wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    &
157                  wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    &
158                  wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    &
159                  wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        &
160                  wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         &
161                  )
162             v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:)
163
164             up_e=1/de(n+u_lup)*(                                                        &
165                  wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        &
166                  wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      &
167                  wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      &
168                  wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      &
169                  wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          & 
170                  wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          &
171                  wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              &
172                  wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              &
173                  wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            &
174                  wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           &
175                  )
176             v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:)
177
178             up_e=1/de(n+u_ldown)*(                                                    &
179                  wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &
180                  wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    &
181                  wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
182                  wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
183                  wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
184                  wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        &
185                  wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      &
186                  wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    &
187                  wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    &
188                  wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     &
189                  )
190             v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:)
191
192             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i
193             cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt 
194             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt 
195          ENDDO
196       ENDDO
197    END DO
198    !123         continue           
199    !==========================================================================================================
200    DO l = 1,llm
201       DO j=jj_begin-1,jj_end+1
202          DO i=ii_begin-1,ii_end+1
203             n=(j-1)*iim+i
204             if (ne(n,right)*ue(n+u_right,l)>0) then
205                ed = cc(n+u_right,l,:) - xyz_i(n,:)
206                qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
207             else
208                ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:)
209                qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
210             endif
211             if (ne(n,lup)*ue(n+u_lup,l)>0) then
212                ed = cc(n+u_lup,l,:) - xyz_i(n,:)
213                qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)
214             else
215                ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
216                qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
217             endif
218             if (ne(n,ldown)*ue(n+u_ldown,l)>0) then
219                ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
220                qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed) 
221             else
222                ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
223                qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
224             endif
225          end do
226       end do
227    END DO
228
229
230    DO j=jj_begin-1,jj_end+1
231       DO i=ii_begin-1,ii_end+1
232          n=(j-1)*iim+i 
233          flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right)
234          flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup) 
235          flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown)
236       ENDDO
237    ENDDO
238
239    DO j=jj_begin,jj_end
240       DO i=ii_begin,ii_end
241          n=(j-1)*iim+i 
242
243          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) &
244               + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) &
245               + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) ) 
246
247          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       & 
248               +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    &
249               - flxx(n+u_left,:) - flxx(n+u_rdown,:) )   
250          him_old(:) = him(n,:) 
251          him(n,:) = him(n,:) + dhi(n,:) 
252          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:) 
253
254       END DO
255    END DO
256
257  CONTAINS
258    !====================================================================================
259    REAL FUNCTION sum2(a1,a2)
260      IMPLICIT NONE
261      REAL,INTENT(IN):: a1(3), a2(3)
262      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
263    END FUNCTION sum2
264
265  END SUBROUTINE compute_advect_horiz
266  !==========================================================================
267  SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det)
268    USE geometry
269    USE metric
270    USE dimensions
271    IMPLICIT NONE
272    INTEGER, INTENT(IN) :: n0,n1,n2,n3
273    REAL,INTENT(IN)     :: q(iim*jjm)
274    REAL,INTENT(OUT)    :: dq(3)
275    REAL(rstd)  ::det,detx,dety,detz
276    INTEGER    :: info
277    INTEGER    :: IPIV(3)
278
279    REAL(rstd) :: A(3,3)
280    REAL(rstd) :: B(3)
281
282    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) 
283    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) 
284    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);             A(3,3)= xyz_v(n3,3)
285
286
287    dq(1) = q(n1)-q(n0)
288    dq(2) = q(n2)-q(n0)
289    dq(3) = 0.0
290    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
291    CALL determinant(A(:,1),A(:,2),A(:,3),det)
292    CALL determinant(dq,A(:,2),A(:,3),detx)
293    CALL determinant(A(:,1),dq,A(:,3),dety)
294    CALL determinant(A(:,1),A(:,2),dq,detz)
295    dq(1) = detx
296    dq(2) = dety
297    dq(3) = detz 
298  END SUBROUTINE gradq
299  !==========================================================================
300  SUBROUTINE determinant(a1,a2,a3,det)
301    IMPLICIT NONE
302    REAL, DIMENSION(3) :: a1, a2,a3
303    REAL ::  det,x1,x2,x3
304    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
305    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
306    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
307    det =  x1 - x2 + x3
308  END SUBROUTINE determinant
309
310END MODULE advect_mod
Note: See TracBrowser for help on using the repository browser.