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

Last change on this file since 583 was 583, checked in by dubos, 7 years ago

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