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

Last change on this file since 186 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
File size: 9.9 KB
Line 
1MODULE 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_one_over_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
24CONTAINS
25
26  SUBROUTINE init_advect_tracer
27    USE advect_mod
28    REAL(rstd),POINTER :: tangent(:,:)
29    REAL(rstd),POINTER :: normal(:,:)
30    REAL(rstd),POINTER :: one_over_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_one_over_sqrt_leng,field_t,type_real, name='one_over_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       one_over_sqrt_leng=f_one_over_sqrt_leng(ind)
50       CALL init_advect(normal,tangent,one_over_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(:,:), one_over_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.5*dt*itau_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          one_over_sqrt_leng=f_one_over_sqrt_leng(ind)
148          CALL compute_gradq3d(q(:,:,k),one_over_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
322END MODULE advect_tracer_mod
Note: See TracBrowser for help on using the repository browser.