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

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

Various optimisations

  • dissipation is not called every timestep (similar way than LMDZ)
  • transfert size of halos has been reduced : need to synchronise redondant data between tiles at itau_sync timestep

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