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

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

Set constant sign for wind way :
ne(ij,right)==ne_right=1
ne(ij,rup)==ne_rup=-1
ne(ij,lup)==ne_lup=1
ne(ij,left)==ne_left=-1
ne(ij,ldown)==ne_ldown=1
ne(ij,rdown)==ne_rdown=-1

Modified transfert function to be compliant for this convention.

YM
File size: 8.3 KB
Line 
1MODULE advect_tracer_mod
2  USE icosa
3  IMPLICIT NONE
4  PRIVATE
5
6  TYPE(t_field),POINTER :: f_normal(:)
7  TYPE(t_field),POINTER :: f_tangent(:)
8  TYPE(t_field),POINTER :: f_gradq3d(:)
9  TYPE(t_field),POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach)
10
11  REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz
12
13  PUBLIC init_advect_tracer, advect_tracer
14
15CONTAINS
16
17  SUBROUTINE init_advect_tracer
18    USE advect_mod
19    REAL(rstd),POINTER :: tangent(:,:)
20    REAL(rstd),POINTER :: normal(:,:)
21    INTEGER :: ind
22
23    CALL allocate_field(f_normal,field_u,type_real,3, name='normal')
24    CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent')
25    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d')
26    CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
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_hfluxt, f_wfluxt,f_u, f_q,f_rhodz)
39    USE advect_mod
40    USE mpipara
41    USE trace
42    USE write_field
43    IMPLICIT NONE
44   
45    TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux
46    TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux
47    TYPE(t_field),POINTER :: f_u(:)        ! velocity (for back-trajectories)
48    TYPE(t_field),POINTER :: f_q(:)        ! tracer
49    TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step
50
51    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), gradq3d(:,:,:), cc(:,:,:)
52    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
53    REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
54    INTEGER :: ind,k
55
56    CALL trace_start("advect_tracer") 
57
58    CALL transfert_request(f_u,req_e1_vect)
59!    CALL transfert_request(f_hfluxt,req_e1)  ! BUG : This (unnecessary) transfer makes the computation go wrong
60    CALL transfert_request(f_wfluxt,req_i1)
61    CALL transfert_request(f_q,req_i1)
62    CALL transfert_request(f_rhodz,req_i1)
63       
64    IF (is_mpi_root) PRINT *, 'Advection scheme'
65
66!    DO ind=1,ndomain
67!       CALL swap_dimensions(ind)
68!       CALL swap_geometry(ind)
69!       normal  = f_normal(ind)
70!       tangent = f_tangent(ind)
71!       cc      = f_cc(ind)
72!       u       = f_u(ind)
73!       q       = f_q(ind)
74!       rhodz   = f_rhodz(ind)
75!       hfluxt  = f_hfluxt(ind)
76!       wfluxt  = f_wfluxt(ind)
77!       gradq3d = f_gradq3d(ind)
78!
79!       ! 1/2 vertical transport
80!       DO k = 1, nqtot
81!          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k))
82!       END DO
83!
84!       ! horizontal transport
85!       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc)
86!       DO k = 1,nqtot
87!          CALL compute_gradq3d(q(:,:,k),gradq3d)
88!          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
89!       END DO
90!
91!       ! 1/2 vertical transport
92!       DO k = 1,nqtot
93!          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
94!       END DO
95!    END DO
96
97    ! 1/2 vertical transport + back-trajectories
98    DO ind=1,ndomain
99       CALL swap_dimensions(ind)
100       CALL swap_geometry(ind)
101       normal  = f_normal(ind)
102       tangent = f_tangent(ind)
103       cc      = f_cc(ind)
104       u       = f_u(ind)
105       q       = f_q(ind)
106       rhodz   = f_rhodz(ind)
107       wfluxt  = f_wfluxt(ind) 
108       DO k = 1, nqtot
109          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k))
110       END DO
111       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
112    END DO
113
114    CALL transfert_request(f_q,req_i1)      ! necessary ?
115    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
116
117    ! horizontal transport - split in two to place transfer of gradq3d
118
119    DO k = 1, nqtot
120       DO ind=1,ndomain
121          CALL swap_dimensions(ind)
122          CALL swap_geometry(ind)
123          q       = f_q(ind)
124          gradq3d = f_gradq3d(ind)
125          CALL compute_gradq3d(q(:,:,k),gradq3d)
126       END DO
127
128       CALL transfert_request(f_gradq3d,req_i1)
129
130       DO ind=1,ndomain
131          CALL swap_dimensions(ind)
132          CALL swap_geometry(ind)
133          cc      = f_cc(ind)
134          q       = f_q(ind)
135          rhodz   = f_rhodz(ind)
136          hfluxt  = f_hfluxt(ind) 
137          gradq3d = f_gradq3d(ind)
138          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
139       END DO
140    END DO
141   
142    CALL transfert_request(f_q,req_i1)      ! necessary ?
143    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
144   
145    ! 1/2 vertical transport
146    DO ind=1,ndomain
147       CALL swap_dimensions(ind)
148       CALL swap_geometry(ind)
149       q       = f_q(ind)
150       rhodz   = f_rhodz(ind)
151       wfluxt  = f_wfluxt(ind) 
152       DO k = 1,nqtot
153          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
154       END DO
155    END DO
156
157    CALL trace_end("advect_tracer")
158
159  END SUBROUTINE advect_tracer
160
161  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q)
162    !
163    !     Auteurs:   P.Le Van, F.Hourdin, F.Forget, T. Dubos
164    !
165    !    ********************************************************************
166    !     Update tracers using vertical mass flux only
167    !     Van Leer scheme with minmod limiter
168    !     wfluxt >0 for upward transport
169    !    ********************************************************************
170    IMPLICIT NONE
171    LOGICAL, INTENT(IN)       :: update_mass
172    REAL(rstd), INTENT(IN)    :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux
173    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
174    REAL(rstd), INTENT(INOUT) :: q(iim*jjm,llm)
175
176    REAL(rstd) :: dq(iim*jjm,llm), & ! increase of q
177         dzqw(iim*jjm,llm),        & ! vertical finite difference of q
178         adzqw(iim*jjm,llm),       & ! abs(dzqw)
179         dzq(iim*jjm,llm),         & ! limited slope of q
180         wq(iim*jjm,llm+1)           ! time-integrated flux of q
181    REAL(rstd) :: dzqmax, newmass, sigw, qq, w
182    INTEGER :: i,ij,l,j
183
184    ! finite difference of q
185    DO l=2,llm
186       DO j=jj_begin-1,jj_end+1
187          DO i=ii_begin-1,ii_end+1
188             ij=(j-1)*iim+i
189             dzqw(ij,l)=q(ij,l)-q(ij,l-1)
190             adzqw(ij,l)=abs(dzqw(ij,l))
191          ENDDO
192       ENDDO
193    ENDDO
194
195    ! minmod-limited slope of q
196    ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l
197    DO l=2,llm-1
198       DO j=jj_begin-1,jj_end+1
199          DO i=ii_begin-1,ii_end+1
200             ij=(j-1)*iim+i
201             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
202                dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) )
203                dzqmax    = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) )
204                dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) )  ! NB : sign(a,b)=a*sign(b)
205             ELSE
206                dzq(ij,l)=0.
207             ENDIF
208          ENDDO
209       ENDDO
210    ENDDO
211
212    ! 0 slope in top and bottom layers
213    DO j=jj_begin-1,jj_end+1
214       DO i=ii_begin-1,ii_end+1
215          ij=(j-1)*iim+i
216          dzq(ij,1)=0.
217          dzq(ij,llm)=0.
218       ENDDO
219    ENDDO
220
221    ! sigw = fraction of mass that leaves level l/l+1
222    ! then amount of q leaving level l/l+1 = wq = w * qq
223    DO l = 1,llm-1
224       DO j=jj_begin-1,jj_end+1
225          DO i=ii_begin-1,ii_end+1
226             ij=(j-1)*iim+i
227             w = fac*wfluxt(ij,l+1)
228             IF(w>0.) THEN  ! upward transport, upwind side is at level l
229                sigw       = w/mass(ij,l)
230                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
231             ELSE           ! downward transport, upwind side is at level l+1
232                sigw       = w/mass(ij,l+1)
233                qq         = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0               
234             ENDIF
235             wq(ij,l+1) = w*qq
236          ENDDO
237       ENDDO
238    END DO
239    ! wq = 0 at top and bottom
240    DO j=jj_begin-1,jj_end+1
241       DO i=ii_begin-1,ii_end+1
242          ij=(j-1)*iim+i
243          wq(ij,llm+1)=0.
244          wq(ij,1)=0.
245       ENDDO
246    END DO
247
248    ! update q, mass is updated only after all q's have been updated
249    DO l=1,llm
250       DO j=jj_begin-1,jj_end+1
251          DO i=ii_begin-1,ii_end+1
252             ij=(j-1)*iim+i
253             newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1))
254             q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass
255             IF(update_mass) mass(ij,l)=newmass
256          ENDDO
257       ENDDO
258    END DO
259
260  END SUBROUTINE vlz
261
262END MODULE advect_tracer_mod
Note: See TracBrowser for help on using the repository browser.