source: codes/icosagcm/trunk/src/advect_tracer.f90 @ 67

Last change on this file since 67 was 30, checked in by sdubey, 12 years ago

bug fix for advection, with transfer
File size: 12.0 KB
Line 
1MODULE advect_tracer_mod
2  USE icosa
3  PRIVATE
4  INTEGER,PARAMETER::iapp_tracvl= 3
5  REAL(rstd),SAVE :: dt
6
7  TYPE(t_field),POINTER :: f_normal(:)
8  TYPE(t_field),POINTER :: f_tangent(:)
9  TYPE(t_field),POINTER :: f_gradq3d(:)
10
11  PUBLIC init_advect_tracer, advect_tracer
12
13CONTAINS
14
15  SUBROUTINE init_advect_tracer(dt_in)
16    USE advect_mod
17    IMPLICIT NONE
18    REAL(rstd),INTENT(IN) :: dt_in
19    REAL(rstd),POINTER :: tangent(:,:)
20    REAL(rstd),POINTER :: normal(:,:)
21    INTEGER :: ind
22
23    dt=dt_in
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    DO ind=1,ndomain
29       CALL swap_dimensions(ind)
30       CALL swap_geometry(ind)
31       normal=f_normal(ind)
32       tangent=f_tangent(ind)
33       CALL init_advect(normal,tangent)
34    END DO
35
36  END SUBROUTINE init_advect_tracer
37
38  SUBROUTINE advect_tracer(f_ps,f_u,f_q)
39    USE icosa
40    USE advect_mod
41    USE disvert_mod
42    IMPLICIT NONE
43    TYPE(t_field),POINTER :: f_ps(:)
44    TYPE(t_field),POINTER :: f_u(:)
45    TYPE(t_field),POINTER :: f_q(:)
46    REAL(rstd),POINTER :: q(:,:,:)
47    REAL(rstd),POINTER :: u(:,:) 
48    REAL(rstd),POINTER :: ps(:) 
49    REAL(rstd),POINTER :: massflx(:,:)
50    REAL(rstd),POINTER :: rhodz(:,:)
51    TYPE(t_field),POINTER,SAVE :: f_massflxc(:)
52    TYPE(t_field),POINTER,SAVE :: f_massflx(:)
53    TYPE(t_field),POINTER,SAVE :: f_uc(:)
54    TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:)
55    TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
56    REAL(rstd),POINTER,SAVE :: massflxc(:,:)
57    REAL(rstd),POINTER,SAVE :: uc(:,:)
58    REAL(rstd),POINTER,SAVE :: rhodzm1(:,:)
59    REAL(rstd):: bigt
60    INTEGER :: ind,it,i,j,l,ij
61    INTEGER,SAVE :: iadvtr=0
62    LOGICAL,SAVE:: first=.TRUE. 
63!------------------------------------------------------sarvesh
64       CALL transfert_request(f_ps,req_i1)
65       CALL transfert_request(f_u,req_e1)
66       CALL transfert_request(f_u,req_e1)
67       CALL transfert_request(f_q,req_i1)
68       
69    IF ( first ) THEN
70       CALL allocate_field(f_rhodz,field_t,type_real,llm) 
71       CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 
72       CALL allocate_field(f_massflxc,field_u,type_real,llm) 
73       CALL allocate_field(f_massflx,field_u,type_real,llm) 
74       CALL allocate_field(f_uc,field_u,type_real,llm) 
75       first = .FALSE.
76    END IF
77
78    DO ind=1,ndomain
79       CALL swap_dimensions(ind)
80       CALL swap_geometry(ind)
81       rhodz=f_rhodz(ind)
82       massflx=f_massflx(ind)
83       ps=f_ps(ind)
84       u=f_u(ind)
85
86       DO l = 1, llm
87          DO j=jj_begin-1,jj_end+1
88             DO i=ii_begin-1,ii_end+1
89                ij=(j-1)*iim+i
90                rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
91             ENDDO
92          ENDDO
93       ENDDO
94
95       DO l = 1, llm
96          DO j=jj_begin-1,jj_end+1
97             DO i=ii_begin-1,ii_end+1
98                ij=(j-1)*iim+i
99                massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
100                massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
101                massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
102             ENDDO
103          ENDDO
104       ENDDO
105    ENDDO
106
107    IF ( iadvtr == 0 ) THEN
108       DO ind=1,ndomain         
109          CALL swap_dimensions(ind)
110          CALL swap_geometry(ind)
111          rhodz=f_rhodz(ind) 
112          rhodzm1 = f_rhodzm1(ind) 
113          massflxc = f_massflxc(ind)  ! accumulated mass fluxes
114          uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura)
115          rhodzm1 = rhodz
116          massflxc = 0.0 
117          uc = 0.0 
118       END DO
119       CALL transfert_request(f_rhodzm1,req_i1)   !---->
120       CALL transfert_request(f_massflxc,req_e1) !---->
121       CALL transfert_request(f_massflxc,req_e1) !------>
122       CALL transfert_request(f_uc,req_e1) !---->
123       CALL transfert_request(f_uc,req_e1) 
124    END IF
125
126    iadvtr = iadvtr + 1 
127
128    DO ind=1,ndomain
129       CALL swap_dimensions(ind)
130       CALL swap_geometry(ind)
131       massflx  = f_massflx(ind)
132       massflxc = f_massflxc(ind) 
133       uc = f_uc(ind) 
134       u  = f_u(ind) 
135       massflxc = massflxc + massflx 
136       uc = uc + u 
137    END DO
138
139    IF ( iadvtr == iapp_tracvl ) THEN
140       PRINT *, 'Advection scheme'
141       bigt = dt*iapp_tracvl
142       DO ind=1,ndomain
143          CALL swap_dimensions(ind)
144          CALL swap_geometry(ind)
145          uc = f_uc(ind) 
146          uc = uc/real(iapp_tracvl) 
147          massflxc = f_massflxc(ind)
148          massflxc = massflxc*dt
149          ! now massflx is time-integrated
150       END DO
151
152       CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt) 
153       iadvtr = 0
154    END IF
155  END SUBROUTINE advect_tracer
156
157  SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)
158    USE field_mod
159    USE domain_mod
160    USE dimensions
161    USE grid_param
162    USE geometry
163    USE metric
164    USE advect_mod
165    IMPLICIT NONE
166
167    TYPE(t_field),POINTER :: f_q(:)
168    TYPE(t_field),POINTER :: f_u(:)
169    TYPE(t_field),POINTER :: f_rhodz(:) 
170    TYPE(t_field),POINTER :: f_massflx(:)
171
172    TYPE(t_field),POINTER,SAVE :: f_wg(:)
173    TYPE(t_field),POINTER,SAVE :: f_zm(:)
174    TYPE(t_field),POINTER,SAVE :: f_zq(:)
175
176    REAL(rstd)::bigt,dt
177    REAL(rstd),POINTER :: q(:,:,:)
178    REAL(rstd),POINTER :: u(:,:)
179    REAL(rstd),POINTER :: rhodz(:,:) 
180    REAL(rstd),POINTER :: massflx(:,:) 
181    REAL(rstd),POINTER :: wg(:,:)
182    REAL(rstd),POINTER :: zq(:,:,:) 
183    REAL(rstd),POINTER :: zm(:,:) 
184    REAL(rstd),POINTER :: tangent(:,:)
185    REAL(rstd),POINTER :: normal(:,:)
186    REAL(rstd),POINTER :: gradq3d(:,:,:)
187
188    REAL(rstd)::pente_max
189    LOGICAL,SAVE::first = .true. 
190    INTEGER :: i,ij,l,j,ind,k
191    REAL(rstd) :: zzpbar, zzw 
192    REAL::qvmax,qvmin 
193
194    IF ( first ) THEN
195       CALL allocate_field(f_wg,field_t,type_real,llm)
196       CALL allocate_field(f_zm,field_t,type_real,llm)
197       CALL allocate_field(f_zq,field_t,type_real,llm,nqtot) 
198       first = .FALSE.
199    END IF
200
201    DO ind=1,ndomain
202       CALL swap_dimensions(ind)
203       CALL swap_geometry(ind)
204       q=f_q(ind) 
205       rhodz=f_rhodz(ind) 
206       zq=f_zq(ind)
207       zm=f_zm(ind) 
208       zm = rhodz ; zq = q 
209       wg = f_wg(ind)
210       wg = 0.0
211       massflx=f_massflx(ind) 
212       CALL advtrac(massflx,wg) ! compute vertical mass fluxes
213    END DO
214
215    DO ind=1,ndomain
216       CALL swap_dimensions(ind)
217       CALL swap_geometry(ind)
218       zq=f_zq(ind)
219       zm=f_zm(ind) 
220       wg=f_wg(ind)
221       wg=wg*0.5
222       DO k = 1, nqtot
223          CALL vlz(zq(:,:,k),2.,zm,wg)
224       END DO
225    END DO
226
227    DO ind=1,ndomain
228       CALL swap_dimensions(ind)
229       CALL swap_geometry(ind)
230       zq=f_zq(ind)
231       zq = f_zq(ind)
232       zm = f_zm(ind)
233       massflx =f_massflx(ind)
234       u = f_u(ind)
235       tangent = f_tangent(ind)
236       normal = f_normal(ind)
237       gradq3d = f_gradq3d(ind)
238
239       DO k = 1,nqtot
240          CALL compute_gradq3d(zq(:,:,k),gradq3d)
241          CALL compute_advect_horiz(tangent,normal,zq(:,:,k),gradq3d,zm,u,massflx,bigt) 
242       END DO
243    END DO
244
245    DO ind=1,ndomain
246       CALL swap_dimensions(ind)
247       CALL swap_geometry(ind)
248       q = f_q(ind)
249       zq = f_zq(ind) 
250       zm = f_zm(ind) 
251       wg = f_wg(ind)
252       DO k = 1,nqtot
253          CALL vlz(zq(:,:,k),2.,zm,wg)
254       END DO
255       q = zq 
256    END DO
257
258  END SUBROUTINE vlsplt
259
260  !============================================================================== 
261  SUBROUTINE advtrac(massflx,wgg)
262    USE domain_mod
263    USE dimensions
264    USE grid_param
265    USE geometry
266    USE metric
267    USE disvert_mod
268    IMPLICIT NONE
269    REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm)
270    REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm)
271
272    INTEGER :: i,j,ij,l
273    REAL(rstd) :: convm(iim*jjm,llm) 
274
275    DO l = 1, llm
276       DO j=jj_begin,jj_end
277          DO i=ii_begin,ii_end
278             ij=(j-1)*iim+i
279             ! divergence of horizontal flux
280             convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  &
281                  ne(ij,rup)*massflx(ij+u_rup,l)       +  &
282                  ne(ij,lup)*massflx(ij+u_lup,l)       +  &
283                  ne(ij,left)*massflx(ij+u_left,l)     +  &
284                  ne(ij,ldown)*massflx(ij+u_ldown,l)   +  &
285                  ne(ij,rdown)*massflx(ij+u_rdown,l))
286          ENDDO
287       ENDDO
288    ENDDO
289
290    ! accumulate divergence from top of model
291    DO  l = llm-1, 1, -1
292       DO j=jj_begin,jj_end
293          DO i=ii_begin,ii_end
294             ij=(j-1)*iim+i
295             convm(ij,l) = convm(ij,l) + convm(ij,l+1)
296          ENDDO
297       ENDDO
298    ENDDO
299
300!!! Compute vertical velocity
301    DO l = 1,llm-1
302       DO j=jj_begin,jj_end
303          DO i=ii_begin,ii_end
304             ij=(j-1)*iim+i
305             wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )) 
306          ENDDO
307       ENDDO
308    ENDDO
309
310    DO j=jj_begin,jj_end
311       DO i=ii_begin,ii_end
312          ij=(j-1)*iim+i
313          wgg(ij,1)  = 0.
314       ENDDO
315    ENDDO
316  END SUBROUTINE advtrac
317
318  SUBROUTINE vlz(q,pente_max,masse,wgg)
319    !c
320    !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
321    !c
322    !c    ********************************************************************
323    !c     Shema  d'advection " pseudo amont " .
324    !c    ********************************************************************
325    USE icosa
326    IMPLICIT NONE
327    !c
328    !c   Arguments:
329    !c   ----------
330    REAL masse(iim*jjm,llm),pente_max
331    REAL q(iim*jjm,llm)
332    REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1) 
333    REAL dq(iim*jjm,llm)
334    INTEGER i,ij,l,j,ii
335    !c
336    REAL wq(iim*jjm,llm+1),newmasse
337
338    REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax
339    REAL sigw
340
341    REAL      SSUM
342
343
344    w(:,1:llm) = -wgg(:,:)  ! w>0 for downward transport
345    w(:,llm+1) = 0.0 
346
347    !c    On oriente tout dans le sens de la pression c'est a dire dans le
348    !c    sens de W
349
350    DO l=2,llm
351       DO j=jj_begin,jj_end
352          DO i=ii_begin,ii_end
353             ij=(j-1)*iim+i
354             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
355             adzqw(ij,l)=abs(dzqw(ij,l))
356          ENDDO
357       ENDDO
358    ENDDO
359
360    DO l=2,llm-1
361       DO j=jj_begin,jj_end
362          DO i=ii_begin,ii_end
363             ij=(j-1)*iim+i
364             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
365                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
366             ELSE
367                dzq(ij,l)=0.
368             ENDIF
369             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
370             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
371          ENDDO
372       ENDDO
373    ENDDO
374
375    DO l=2,llm-1
376       DO j=jj_begin,jj_end
377          DO i=ii_begin,ii_end
378             ij=(j-1)*iim+i
379             dzq(ij,1)=0.
380             dzq(ij,llm)=0.
381          ENDDO
382       ENDDO
383    ENDDO
384    !c ---------------------------------------------------------------
385    !c   .... calcul des termes d'advection verticale  .......
386    !c ---------------------------------------------------------------
387
388    !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
389
390    DO l = 1,llm-1
391       DO j=jj_begin,jj_end
392          DO i=ii_begin,ii_end
393             ij=(j-1)*iim+i
394             IF(w(ij,l+1).gt.0.) THEN
395                ! upwind only if downward transport
396                sigw=w(ij,l+1)/masse(ij,l+1)
397                wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
398             ELSE
399                ! upwind only if upward transport
400                sigw=w(ij,l+1)/masse(ij,l)
401                wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
402             ENDIF
403          ENDDO
404       ENDDO
405    END DO
406
407    DO j=jj_begin,jj_end
408       DO i=ii_begin,ii_end
409          ij=(j-1)*iim+i
410          wq(ij,llm+1)=0.
411          wq(ij,1)=0.
412       ENDDO
413    END DO
414
415    DO l=1,llm
416       DO j=jj_begin,jj_end
417          DO i=ii_begin,ii_end
418             ij=(j-1)*iim+i
419             ! masse -= dw/dz but w>0 <=> downward
420             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l)) 
421!             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>>
422             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse   
423!             masse(ij,l)=newmasse
424          ENDDO
425       ENDDO
426    END DO
427    RETURN
428  END SUBROUTINE vlz
429
430END MODULE advect_tracer_mod
Note: See TracBrowser for help on using the repository browser.