source: codes/icosagcm/devel/src/transport/advect_tracer.f90 @ 915

Last change on this file since 915 was 856, checked in by dubos, 5 years ago

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