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

Last change on this file since 156 was 151, checked in by ymipsl, 11 years ago

Implementation of mixte parallelism MPI/OpenMP into src directory

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