source: codes/icosagcm/trunk/src/transport/advect_tracer.F90 @ 954

Last change on this file since 954 was 954, checked in by adurocher, 5 years ago

trunk : Added metric terms to kernels parameters to avoid Host/GPU transferts

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