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

Last change on this file since 327 was 327, checked in by ymipsl, 9 years ago

Merge recent developments from saturn branch onto trunk.

  • lmdz generic physics interface
  • performance improvment on mix mpi/openmp
  • asynchrone and overlaping communication
  • best domain distribution between process and threads
  • ....

This version is compatible with the actual saturn version and the both branches are considered merged on dynamico component.

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