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

Last change on this file since 21 was 19, checked in by ymipsl, 12 years ago

Simplify the management of the module.

YM
File size: 13.6 KB
Line 
1      MODULE advect_mod
2   USE icosa
3        IMPLICIT NONE
4
5
6  TYPE(t_field),POINTER :: f_normal(:)
7  TYPE(t_field),POINTER :: f_tangent(:)
8  TYPE(t_field),POINTER ::  f_gradq3d(:)
9  REAL(rstd),POINTER :: gradq3d(:,:,:) 
10  REAL(rstd),POINTER :: normal(:,:)
11  REAL(rstd),POINTER :: tangent(:,:)
12                       
13CONTAINS
14         
15       
16  SUBROUTINE allocate_advect
17  USE field_mod
18  USE domain_mod
19  USE dimensions
20  USE geometry
21  USE metric
22  IMPLICIT NONE
23
24    CALL allocate_field(f_normal,field_u,type_real,3)
25    CALL allocate_field(f_tangent,field_u,type_real,3)
26    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3)
27
28  END SUBROUTINE allocate_advect
29!==========================================================================
30  SUBROUTINE swap_advect(ind)
31  USE domain_mod
32  USE dimensions
33  USE geometry
34  USE metric
35  IMPLICIT NONE
36  INTEGER,INTENT(IN) :: ind
37 
38      normal=f_normal(ind)
39      tangent=f_tangent(ind)
40      gradq3d = f_gradq3d(ind) 
41
42  END SUBROUTINE swap_advect
43!==========================================================================
44   
45  SUBROUTINE init_advect
46  USE domain_mod
47  USE dimensions
48  USE geometry
49  USE metric
50  USE vector
51  IMPLICIT NONE
52   INTEGER :: ind,i,j,n
53   
54   CALL allocate_advect
55   
56    DO j=jj_begin-1,jj_end+1
57      DO i=ii_begin-1,ii_end+1
58        n=(j-1)*iim+i
59           
60        CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:))
61        normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right)
62         
63        CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:))
64         normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup)
65       
66        CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))   
67        normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown)
68                     
69        tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:) 
70        tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50)
71             
72        tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:) 
73        tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50)
74     
75        tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:)
76        tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50)
77      END DO
78     ENDDO     
79  END SUBROUTINE init_advect
80!=======================================================================================
81
82
83
84!=======================================================================================
85  SUBROUTINE advect1(qi)
86  USE domain_mod
87  USE dimensions
88  USE geometry
89  USE metric
90  IMPLICIT NONE
91    REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm)
92    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
93    REAL(rstd) :: alphamx,alphami,alpha ,maggrd,leng
94    REAL(rstd) :: leng1,leng2
95    REAL(rstd) :: arr(2*iim*jjm)
96    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
97    REAL(rstd) :: ar
98    INTEGER :: i,j,n,k,ind,l
99!========================================================================================== GRADIENT
100       Do l = 1,llm 
101        DO j=jj_begin-1,jj_end+1
102         DO i=ii_begin-1,ii_end+1
103          n=(j-1)*iim+i   
104          CALL gradq(n,n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up))
105          CALL gradq(n,n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down))
106         END DO
107        END DO
108       END DO
109
110!      Do l =1,llm
111        DO j=jj_begin,jj_end
112          DO i=ii_begin,ii_end
113             n=(j-1)*iim+i
114             gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:)   +   &
115                           gradtri(n+z_rup,:,:) +  gradtri(n+z_ldown,:,:)   +   & 
116                          gradtri(n+z_lup,:,:)+    gradtri(n+z_rdown,:,:) 
117             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)
118             gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50) 
119           END DO
120         END DO   
121!        END DO
122 
123!============================================================================================= LIMITING
124! GO TO 120
125  Do l =1,llm
126   DO j=jj_begin,jj_end
127     DO i=ii_begin,ii_end
128       n=(j-1)*iim+i
129       maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:))
130       maggrd = sqrt(maggrd) 
131       leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), &
132       sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  &
133       sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2))
134       maxq_c = qi(n,l) + maggrd*sqrt(leng)
135       minq_c = qi(n,l) - maggrd*sqrt(leng)
136       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), &
137             qi(n+t_rdown,l),qi(n+t_ldown,l))
138       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), &
139             qi(n+t_rdown,l),qi(n+t_ldown,l))
140       alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) )
141       alphamx = max(alphamx,0.0)
142       alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l))
143       alphami = max(alphami,0.0) 
144       alpha   = min(alphamx,alphami,1.0)
145       gradq3d(n,l,:) = alpha*gradq3d(n,l,:)
146      END DO
147    END DO
148   END DO
149         END SUBROUTINE ADVECT1         
150!===================================================================================================     
151         SUBROUTINE advect2(qi,him,ue,he,bigt)
152  USE domain_mod
153  USE dimensions
154  USE geometry
155  USE metric
156  IMPLICIT NONE
157    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm)
158    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm)
159    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,llm)
160    REAL(rstd),INTENT(IN) :: he(3*iim*jjm,llm)
161    REAL(rstd),INTENT(IN)::bigt
162    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 
163    REAL(rstd) :: cc(3*iim*jjm,llm,3) 
164    REAL(rstd) :: v_e(3*iim*jjm,llm,3)
165    REAL(rstd) :: up_e
166   
167    REAL(rstd) :: qe(3*iim*jjm,llm)
168    REAL(rstd) :: ed(3)   
169    REAL(rstd) :: flxx(3*iim*jjm,llm)
170    INTEGER :: i,j,n,k,l
171    REAL(rstd):: him_old(llm) 
172
173
174!go to  123
175   DO l = 1,llm
176    DO j=jj_begin,jj_end
177       DO i=ii_begin,ii_end
178        n=(j-1)*iim+i
179
180                    up_e =1/de(n+u_right)*(                                                   &
181                         wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
182                         wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
183                         wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      &
184                         wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    &
185                         wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
186                         wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    &
187                         wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    &
188                         wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    &
189                         wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        &
190                         wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         &
191                                        )
192       v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:)
193 
194                up_e=1/de(n+u_lup)*(                                                        &
195                         wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        &
196                         wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      &
197                         wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      &
198                         wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      &
199                         wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          & 
200                         wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          &
201                         wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              &
202                         wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              &
203                         wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            &
204                         wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           &
205                                  )
206       v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:)
207   
208                  up_e=1/de(n+u_ldown)*(                                                    &
209                         wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &
210                         wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    &
211                         wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
212                         wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
213                         wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
214                         wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        &
215                         wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      &
216                         wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    &
217                         wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    &
218                         wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     &
219                                      )
220         v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:)
221
222          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
223          cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt 
224          cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt 
225       ENDDO
226     ENDDO 
227    END DO
228!123     continue           
229!==========================================================================================================
230    DO l = 1,llm
231        DO j=jj_begin-1,jj_end+1
232         DO i=ii_begin-1,ii_end+1
233            n=(j-1)*iim+i
234        if (ne(n,right)*ue(n+u_right,l)>0) then
235          ed = cc(n+u_right,l,:) - xyz_i(n,:)
236         qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
237        else
238           ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:)
239         qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
240        endif
241        if (ne(n,lup)*ue(n+u_lup,l)>0) then
242           ed = cc(n+u_lup,l,:) - xyz_i(n,:)
243         qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)
244        else
245           ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
246         qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
247        endif
248        if (ne(n,ldown)*ue(n+u_ldown,l)>0) then
249          ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
250         qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed) 
251        else
252         ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
253        qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
254        endif
255        end do
256    end do
257   END DO
258
259
260         DO j=jj_begin-1,jj_end+1
261           DO i=ii_begin-1,ii_end+1
262            n=(j-1)*iim+i 
263            flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right)
264            flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup) 
265            flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown)
266           ENDDO
267         ENDDO
268                               
269     DO j=jj_begin,jj_end
270       DO i=ii_begin,ii_end
271        n=(j-1)*iim+i 
272
273         dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) &
274                       + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) &
275                       + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) ) 
276
277         dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       & 
278                         +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    &
279                         - flxx(n+u_left,:) - flxx(n+u_rdown,:) )   
280         him_old(:) = him(n,:) 
281         him(n,:) = him(n,:) + dhi(n,:) 
282         qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:) 
283                       
284       END DO
285     END DO
286
287        CONTAINS
288!====================================================================================
289          REAL FUNCTION sum2(a1,a2)
290        IMPLICIT NONE
291        REAL,INTENT(IN):: a1(3), a2(3)
292                sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
293           END FUNCTION sum2
294 
295                END SUBROUTINE advect2 
296!==========================================================================
297  SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det)
298  USE geometry
299  USE metric
300  USE dimensions
301  IMPLICIT NONE
302    INTEGER, INTENT(IN) :: n0,n1,n2,n3
303    REAL,INTENT(IN)     :: q(iim*jjm)
304    REAL,INTENT(OUT)    :: dq(3)
305    REAL(rstd)  ::det,detx,dety,detz
306    INTEGER    :: info
307    INTEGER    :: IPIV(3)
308
309    REAL(rstd) :: A(3,3)
310    REAL(rstd) :: B(3)
311
312    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) 
313    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) 
314    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);             A(3,3)= xyz_v(n3,3)
315   
316 
317    dq(1) = q(n1)-q(n0)
318    dq(2) = q(n2)-q(n0)
319    dq(3) = 0.0
320!    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
321         CALL determinant(A(:,1),A(:,2),A(:,3),det)
322         CALL determinant(dq,A(:,2),A(:,3),detx)
323         CALL determinant(A(:,1),dq,A(:,3),dety)
324         CALL determinant(A(:,1),A(:,2),dq,detz)
325         dq(1) = detx
326         dq(2) = dety
327         dq(3) = detz 
328  END SUBROUTINE gradq
329!==========================================================================
330          SUBROUTINE determinant(a1,a2,a3,det)
331        IMPLICIT NONE
332        REAL, DIMENSION(3) :: a1, a2,a3
333        REAL ::  det,x1,x2,x3
334        x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
335        x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
336        x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
337        det =  x1 - x2 + x3
338        END SUBROUTINE   
339        END MODULE 
Note: See TracBrowser for help on using the repository browser.