Context Navigation

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

Last change on this file since 252 was 252, checked in by dubos, 10 years ago
Fix recently introduced bug : slope limiter
File size: 9.8 KB

Line
1	MODULE advect_tracer_mod
2	USE icosa
3	IMPLICIT NONE
4	PRIVATE
5
6	TYPE(t_field),SAVE,POINTER :: f_normal(:)
7	TYPE(t_field),SAVE,POINTER :: f_tangent(:)
8	TYPE(t_field),SAVE,POINTER :: f_gradq3d(:)
9	TYPE(t_field),SAVE,POINTER :: f_cc(:) ! starting point of backward-trajectory (Miura approach)
10	TYPE(t_field),SAVE,POINTER :: f_sqrt_leng(:)
11
12	TYPE(t_message),SAVE :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d
13
14	REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz
15
16	! temporary shared variable for vlz
17	TYPE(t_field),SAVE,POINTER :: f_dzqw(:) ! vertical finite difference of q
18	TYPE(t_field),SAVE,POINTER :: f_adzqw(:) ! abs(dzqw)
19	TYPE(t_field),SAVE,POINTER :: f_dzq(:) ! limited slope of q
20	TYPE(t_field),SAVE,POINTER :: f_wq(:) ! time-integrated flux of q
21
22	PUBLIC init_advect_tracer, advect_tracer
23
24	CONTAINS
25
26	SUBROUTINE init_advect_tracer
27	USE advect_mod
28	REAL(rstd),POINTER :: tangent(:,:)
29	REAL(rstd),POINTER :: normal(:,:)
30	REAL(rstd),POINTER :: sqrt_leng(:)
31	INTEGER :: ind
32
33	CALL allocate_field(f_normal,field_u,type_real,3, name='normal')
34	CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent')
35	CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d')
36	CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
37	CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng')
38	CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw')
39	CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw')
40	CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq')
41	CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq')
42
43	DO ind=1,ndomain
44	IF (.NOT. assigned_domain(ind)) CYCLE
45	CALL swap_dimensions(ind)
46	CALL swap_geometry(ind)
47	normal=f_normal(ind)
48	tangent=f_tangent(ind)
49	sqrt_leng=f_sqrt_leng(ind)
50	CALL init_advect(normal,tangent,sqrt_leng)
51	END DO
52
53	END SUBROUTINE init_advect_tracer
54
55	SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz)
56	USE advect_mod
57	USE mpipara
58	USE trace
59	USE write_field
60	IMPLICIT NONE
61
62	TYPE(t_field),POINTER :: f_hfluxt(:) ! time-integrated horizontal mass flux
63	TYPE(t_field),POINTER :: f_wfluxt(:) ! time-integrated vertical mass flux
64	TYPE(t_field),POINTER :: f_u(:) ! velocity (for back-trajectories)
65	TYPE(t_field),POINTER :: f_q(:) ! tracer
66	TYPE(t_field),POINTER :: f_rhodz(:) ! mass field at beginning of macro time step
67
68	REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
69	REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
70	REAL(rstd),POINTER :: rhodz(:,:), u(:,:)
71	! temporary shared variable for vlz
72	REAL(rstd),POINTER :: dzqw(:,:) ! vertical finite difference of q
73	REAL(rstd),POINTER :: adzqw(:,:) ! abs(dzqw)
74	REAL(rstd),POINTER :: dzq(:,:) ! limited slope of q
75	REAL(rstd),POINTER :: wq(:,:) ! time-integrated flux of q
76
77	INTEGER :: ind,k
78	LOGICAL,SAVE :: first=.TRUE.
79	!$OMP THREADPRIVATE(first)
80
81	IF (first) THEN
82	first=.FALSE.
83	CALL init_message(f_u,req_e1_vect,req_u)
84	CALL init_message(f_cc,req_e1_scal,req_cc)
85	CALL init_message(f_wfluxt,req_i1,req_wfluxt)
86	CALL init_message(f_q,req_i1,req_q)
87	CALL init_message(f_rhodz,req_i1,req_rhodz)
88	CALL init_message(f_gradq3d,req_i1,req_gradq3d)
89	ENDIF
90
91	!!$OMP BARRIER
92
93	CALL trace_start("advect_tracer")
94
95	CALL send_message(f_u,req_u)
96	CALL wait_message(req_u)
97	CALL send_message(f_wfluxt,req_wfluxt)
98	CALL wait_message(req_wfluxt)
99	CALL send_message(f_q,req_q)
100	CALL wait_message(req_q)
101	CALL send_message(f_rhodz,req_rhodz)
102	CALL wait_message(req_rhodz)
103
104	! CALL wait_message(req_u)
105	! CALL wait_message(req_wfluxt)
106	! CALL wait_message(req_q)
107	! CALL wait_message(req_rhodz)
108
109	! 1/2 vertical transport + back-trajectories
110	DO ind=1,ndomain
111	IF (.NOT. assigned_domain(ind)) CYCLE
112	CALL swap_dimensions(ind)
113	CALL swap_geometry(ind)
114	normal = f_normal(ind)
115	tangent = f_tangent(ind)
116	cc = f_cc(ind)
117	u = f_u(ind)
118	q = f_q(ind)
119	rhodz = f_rhodz(ind)
120	wfluxt = f_wfluxt(ind)
121	dzqw = f_dzqw(ind)
122	adzqw = f_adzqw(ind)
123	dzq = f_dzq(ind)
124	wq = f_wq(ind)
125
126	DO k = 1, nqtot
127	CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k),1,dzqw, adzqw, dzq, wq)
128	END DO
129
130	CALL compute_backward_traj(tangent,normal,u,0.5dtitau_adv, cc)
131
132	END DO
133
134	CALL send_message(f_cc,req_cc)
135	CALL wait_message(req_cc)
136
137
138	! horizontal transport - split in two to place transfer of gradq3d
139	!!$OMP BARRIER
140	DO k = 1, nqtot
141	DO ind=1,ndomain
142	IF (.NOT. assigned_domain(ind)) CYCLE
143	CALL swap_dimensions(ind)
144	CALL swap_geometry(ind)
145	q = f_q(ind)
146	gradq3d = f_gradq3d(ind)
147	sqrt_leng=f_sqrt_leng(ind)
148	CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v)
149	END DO
150
151	CALL send_message(f_gradq3d,req_gradq3d)
152	! CALL wait_message(req_cc)
153	CALL wait_message(req_gradq3d)
154
155
156	DO ind=1,ndomain
157	IF (.NOT. assigned_domain(ind)) CYCLE
158	CALL swap_dimensions(ind)
159	CALL swap_geometry(ind)
160	cc = f_cc(ind)
161	q = f_q(ind)
162	rhodz = f_rhodz(ind)
163	hfluxt = f_hfluxt(ind)
164	gradq3d = f_gradq3d(ind)
165	CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
166	END DO
167	END DO
168
169	! 1/2 vertical transport
170	!!$OMP BARRIER
171
172	DO ind=1,ndomain
173	IF (.NOT. assigned_domain(ind)) CYCLE
174	CALL swap_dimensions(ind)
175	CALL swap_geometry(ind)
176	q = f_q(ind)
177	rhodz = f_rhodz(ind)
178	wfluxt = f_wfluxt(ind)
179	dzqw = f_dzqw(ind)
180	adzqw = f_adzqw(ind)
181	dzq = f_dzq(ind)
182	wq = f_wq(ind)
183
184	DO k = 1,nqtot
185	CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k),0, dzqw, adzqw, dzq, wq)
186	END DO
187
188	END DO
189
190	CALL trace_end("advect_tracer")
191
192	!!$OMP BARRIER
193
194	END SUBROUTINE advect_tracer
195
196	SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q, halo, dzqw, adzqw, dzq, wq)
197	!
198	! Auteurs: P.Le Van, F.Hourdin, F.Forget, T. Dubos
199	!
200	! ********************************************************************
201	! Update tracers using vertical mass flux only
202	! Van Leer scheme with minmod limiter
203	! wfluxt >0 for upward transport
204	! ********************************************************************
205	USE trace
206	USE omp_para
207	IMPLICIT NONE
208	LOGICAL, INTENT(IN) :: update_mass
209	REAL(rstd), INTENT(IN) :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux
210	REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
211	REAL(rstd), INTENT(INOUT) :: q(iim*jjm,llm)
212	INTEGER, INTENT(IN) :: halo
213
214	! temporary shared variable
215	REAL(rstd),INTENT(INOUT) :: dzqw(iim*jjm,llm), & ! vertical finite difference of q
216	adzqw(iim*jjm,llm), & ! abs(dzqw)
217	dzq(iim*jjm,llm), & ! limited slope of q
218	wq(iim*jjm,llm+1) ! time-integrated flux of q
219
220
221	REAL(rstd) :: dzqmax, newmass, sigw, qq, w
222	INTEGER :: i,ij,l,j,ijb,ije
223
224	CALL trace_start("vlz")
225
226	ijb=((jj_begin-halo)-1)*iim+ii_begin-halo
227	ije = ((jj_end+halo)-1)*iim+ii_end+halo
228
229	! finite difference of q
230
231	DO l=ll_beginp1,ll_end
232	!$SIMD
233	DO ij=ijb,ije
234	dzqw(ij,l)=q(ij,l)-q(ij,l-1)
235	adzqw(ij,l)=abs(dzqw(ij,l))
236	ENDDO
237	ENDDO
238
239	!--> flush dzqw, adzqw
240	!!$OMP BARRIER
241
242	! minmod-limited slope of q
243	! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l
244
245	DO l=ll_beginp1,ll_endm1
246	!$SIMD
247	DO ij=ijb,ije
248	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
249	dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) )
250	dzqmax = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) )
251	dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) ) ! NB : sign(a,b)=a*sign(b)
252	ELSE
253	dzq(ij,l)=0.
254	ENDIF
255	ENDDO
256	ENDDO
257
258
259	! 0 slope in top and bottom layers
260	IF (omp_first) THEN
261	DO ij=ijb,ije
262	dzq(ij,1)=0.
263	ENDDO
264	ENDIF
265
266	IF (omp_last) THEN
267	DO ij=ijb,ije
268	dzq(ij,llm)=0.
269	ENDDO
270	ENDIF
271
272	!---> flush dzq
273	!!$OMP BARRIER
274
275	! sigw = fraction of mass that leaves level l/l+1
276	! then amount of q leaving level l/l+1 = wq = w * qq
277	DO l=ll_beginp1,ll_end
278	!$SIMD
279	DO ij=ijb,ije
280	w = fac*wfluxt(ij,l)
281	IF(w>0.) THEN ! upward transport, upwind side is at level l
282	sigw = w/mass(ij,l-1)
283	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
284	ELSE ! downward transport, upwind side is at level l+1
285	sigw = w/mass(ij,l)
286	qq = q(ij,l)-0.5(1.+sigw)dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0
287	ENDIF
288	wq(ij,l) = w*qq
289	ENDDO
290	END DO
291	! wq = 0 at top and bottom
292	IF (omp_first) THEN
293	DO ij=ijb,ije
294	wq(ij,1)=0.
295	END DO
296	ENDIF
297
298	IF (omp_last) THEN
299	DO ij=ijb,ije
300	wq(ij,llm+1)=0.
301	END DO
302	ENDIF
303
304	! --> flush wq
305	!!$OMP BARRIER
306
307
308	! update q, mass is updated only after all q's have been updated
309	DO l=ll_begin,ll_end
310	!$SIMD
311	DO ij=ijb,ije
312	newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1))
313	q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass
314	IF(update_mass) mass(ij,l)=newmass
315	ENDDO
316	END DO
317
318	CALL trace_end("vlz")
319
320	END SUBROUTINE vlz
321
322	END MODULE advect_tracer_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: