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

Last change on this file since 1025 was 1025, checked in by millour, 4 years ago

Change "tracer_mod" to "tracer_icosa_mod" to avoid conflicting module
name with some physics packages.
EM
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_icosa_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    CALL wait_message(req_cc)
151
152    ! horizontal transport - split in two to place transfer of gradq3d
153    DO k = 1, nqtot
154      IF (advection_scheme(k)==advect_vanleer) THEN
155       
156        DO ind=1,ndomain
157          IF (.NOT. assigned_domain(ind)) CYCLE
158          CALL swap_dimensions(ind)
159          CALL swap_geometry(ind)
160          q       = f_q(ind)
161          gradq3d = f_gradq3d(ind)
162          sqrt_leng=f_sqrt_leng(ind)
163          CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v)
164        END DO
165
166        CALL send_message(f_gradq3d,req_gradq3d)
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.