source: palm/trunk/SOURCE/diffusion_s.f90 @ 1015

Last change on this file since 1015 was 1015, checked in by raasch, 9 years ago

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

  • Property svn:keywords set to Id
File size: 16.8 KB
Line 
1 MODULE diffusion_s_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! ------------------
6! accelerator version (*_acc) added
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_s.f90 1015 2012-09-27 09:23:24Z raasch $
11!
12! 1010 2012-09-20 07:59:54Z raasch
13! cpp switch __nopointer added for pointer free version
14!
15! 1001 2012-09-13 14:08:46Z raasch
16! some arrays comunicated by module instead of parameter list
17!
18! 667 2010-12-23 12:06:00Z suehring/gryschka
19! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
20!
21! 183 2008-08-04 15:39:12Z letzel
22! bugfix: calculation of fluxes at vertical surfaces
23!
24! 129 2007-10-30 12:12:24Z letzel
25! replace wall_heatflux by wall_s_flux that is now included in the parameter
26! list, bugfix for assignment of fluxes at walls
27!
28! 20 2007-02-26 00:12:32Z raasch
29! Bugfix: ddzw dimensioned 1:nzt"+1"
30! Calculation extended for gridpoint nzt, fluxes can be given at top,
31! +s_flux_t in parameter list, s_flux renamed s_flux_b
32!
33! RCS Log replace by Id keyword, revision history cleaned up
34!
35! Revision 1.8  2006/02/23 10:34:17  raasch
36! nzb_2d replaced by nzb_s_outer in horizontal diffusion and by nzb_s_inner
37! or nzb_diff_s_inner, respectively, in vertical diffusion, prescribed surface
38! fluxes at vertically oriented topography
39!
40! Revision 1.1  2000/04/13 14:54:02  schroeter
41! Initial revision
42!
43!
44! Description:
45! ------------
46! Diffusion term of scalar quantities (temperature and water content)
47!------------------------------------------------------------------------------!
48
49    PRIVATE
50    PUBLIC diffusion_s, diffusion_s_acc
51
52    INTERFACE diffusion_s
53       MODULE PROCEDURE diffusion_s
54       MODULE PROCEDURE diffusion_s_ij
55    END INTERFACE diffusion_s
56
57    INTERFACE diffusion_s_acc
58       MODULE PROCEDURE diffusion_s_acc
59    END INTERFACE diffusion_s_acc
60
61 CONTAINS
62
63
64!------------------------------------------------------------------------------!
65! Call for all grid points
66!------------------------------------------------------------------------------!
67    SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
68
69       USE arrays_3d
70       USE control_parameters
71       USE grid_variables
72       USE indices
73
74       IMPLICIT NONE
75
76       INTEGER ::  i, j, k
77       REAL    ::  vertical_gridspace
78       REAL    ::  wall_s_flux(0:4)
79       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
80#if defined( __nopointer )
81       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
82#else
83       REAL, DIMENSION(:,:,:), POINTER ::  s
84#endif
85
86       DO  i = nxl, nxr
87          DO  j = nys,nyn
88!
89!--          Compute horizontal diffusion
90             DO  k = nzb_s_outer(j,i)+1, nzt
91
92                tend(k,j,i) = tend(k,j,i)                                     &
93                                          + 0.5 * (                           &
94                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
95                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
96                                                  ) * ddx2                    &
97                                          + 0.5 * (                           &
98                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
99                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
100                                                  ) * ddy2
101             ENDDO
102
103!
104!--          Apply prescribed horizontal wall heatflux where necessary
105             IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
106             THEN
107                DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
108
109                   tend(k,j,i) = tend(k,j,i)                                  &
110                                                + ( fwxp(j,i) * 0.5 *         &
111                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
112                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
113                                                   -fwxm(j,i) * 0.5 *         &
114                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
115                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
116                                                  ) * ddx2                    &
117                                                + ( fwyp(j,i) * 0.5 *         &
118                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
119                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
120                                                   -fwym(j,i) * 0.5 *         &
121                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
122                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
123                                                  ) * ddy2
124                ENDDO
125             ENDIF
126
127!
128!--          Compute vertical diffusion. In case that surface fluxes have been
129!--          prescribed or computed at bottom and/or top, index k starts/ends at
130!--          nzb+2 or nzt-1, respectively.
131             DO  k = nzb_diff_s_inner(j,i), nzt_diff
132
133                tend(k,j,i) = tend(k,j,i)                                     &
134                                       + 0.5 * (                              &
135            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
136          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
137                                               ) * ddzw(k)
138             ENDDO
139
140!
141!--          Vertical diffusion at the first computational gridpoint along
142!--          z-direction
143             IF ( use_surface_fluxes )  THEN
144
145                k = nzb_s_inner(j,i)+1
146
147                tend(k,j,i) = tend(k,j,i)                                     &
148                                       + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )    &
149                                               * ( s(k+1,j,i)-s(k,j,i) )      &
150                                               * ddzu(k+1)                    &
151                                           + s_flux_b(j,i)                    &
152                                         ) * ddzw(k)
153
154             ENDIF
155
156!
157!--          Vertical diffusion at the last computational gridpoint along
158!--          z-direction
159             IF ( use_top_fluxes )  THEN
160
161                k = nzt
162
163                tend(k,j,i) = tend(k,j,i)                                     &
164                                       + ( - s_flux_t(j,i)                    &
165                                           - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
166                                                 * ( s(k,j,i)-s(k-1,j,i) )    &
167                                                 * ddzu(k)                    &
168                                         ) * ddzw(k)
169
170             ENDIF
171
172          ENDDO
173       ENDDO
174
175    END SUBROUTINE diffusion_s
176
177
178!------------------------------------------------------------------------------!
179! Call for all grid points - accelerator version
180!------------------------------------------------------------------------------!
181    SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
182
183       USE arrays_3d
184       USE control_parameters
185       USE grid_variables
186       USE indices
187
188       IMPLICIT NONE
189
190       INTEGER ::  i, j, k
191       REAL    ::  vertical_gridspace
192       REAL    ::  wall_s_flux(0:4)
193       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
194#if defined( __nopointer )
195       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
196#else
197       REAL, DIMENSION(:,:,:), POINTER ::  s
198#endif
199
200       !$acc kernels present( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, kh )        &
201       !$acc         present( nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, s ) &
202       !$acc         present( s_flux_b, s_flux_t, tend, wall_s_flux )         &
203       !$acc         present( wall_w_x, wall_w_y )
204       !$acc loop
205       DO  i = nxl, nxr
206          DO  j = nys,nyn
207!
208!--          Compute horizontal diffusion
209             !$acc loop vector( 32 )
210             DO  k = 1, nzt
211                IF ( k > nzb_s_outer(j,i) )  THEN
212
213                   tend(k,j,i) = tend(k,j,i)                                  &
214                                          + 0.5 * (                           &
215                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
216                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
217                                                  ) * ddx2                    &
218                                          + 0.5 * (                           &
219                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
220                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
221                                                  ) * ddy2
222                ENDIF
223             ENDDO
224
225!
226!--          Apply prescribed horizontal wall heatflux where necessary
227             !$acc loop vector(32)
228             DO  k = 1, nzt
229                IF ( k > nzb_s_inner(j,i)  .AND.  k <= nzb_s_outer(j,i)  .AND. &
230                     ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 ) )    &
231                THEN
232                   tend(k,j,i) = tend(k,j,i)                                  &
233                                                + ( fwxp(j,i) * 0.5 *         &
234                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
235                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
236                                                   -fwxm(j,i) * 0.5 *         &
237                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
238                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
239                                                  ) * ddx2                    &
240                                                + ( fwyp(j,i) * 0.5 *         &
241                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
242                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
243                                                   -fwym(j,i) * 0.5 *         &
244                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
245                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
246                                                  ) * ddy2
247                ENDIF
248             ENDDO
249
250!
251!--          Compute vertical diffusion. In case that surface fluxes have been
252!--          prescribed or computed at bottom and/or top, index k starts/ends at
253!--          nzb+2 or nzt-1, respectively.
254             !$acc loop vector( 32 )
255             DO  k = 1, nzt_diff
256                IF ( k >= nzb_diff_s_inner(j,i) )  THEN
257                   tend(k,j,i) = tend(k,j,i)                                  &
258                                       + 0.5 * (                              &
259            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
260          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
261                                               ) * ddzw(k)
262                ENDIF
263             ENDDO
264
265!
266!--          Vertical diffusion at the first computational gridpoint along
267!--          z-direction
268             !$acc loop vector( 32 )
269             DO  k = 1, nzt
270                IF ( use_surface_fluxes  .AND.  k == nzb_s_inner(j,i)+1 )  THEN
271                   tend(k,j,i) = tend(k,j,i)                                  &
272                                          + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
273                                                  * ( s(k+1,j,i)-s(k,j,i) )   &
274                                                  * ddzu(k+1)                 &
275                                              + s_flux_b(j,i)                 &
276                                            ) * ddzw(k)
277                ENDIF
278
279!
280!--             Vertical diffusion at the last computational gridpoint along
281!--             z-direction
282                IF ( use_top_fluxes  .AND.  k == nzt )  THEN
283                   tend(k,j,i) = tend(k,j,i)                                   &
284                                          + ( - s_flux_t(j,i)                  &
285                                              - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
286                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
287                                                    * ddzu(k)                  &
288                                            ) * ddzw(k)
289                ENDIF
290             ENDDO
291
292          ENDDO
293       ENDDO
294       !$acc end kernels
295
296    END SUBROUTINE diffusion_s_acc
297
298
299!------------------------------------------------------------------------------!
300! Call for grid point i,j
301!------------------------------------------------------------------------------!
302    SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
303
304       USE arrays_3d
305       USE control_parameters
306       USE grid_variables
307       USE indices
308
309       IMPLICIT NONE
310
311       INTEGER ::  i, j, k
312       REAL    ::  vertical_gridspace
313       REAL    ::  wall_s_flux(0:4)
314       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
315#if defined( __nopointer )
316       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
317#else
318       REAL, DIMENSION(:,:,:), POINTER ::  s
319#endif
320
321!
322!--    Compute horizontal diffusion
323       DO  k = nzb_s_outer(j,i)+1, nzt
324
325          tend(k,j,i) = tend(k,j,i)                                           &
326                                          + 0.5 * (                           &
327                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
328                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
329                                                  ) * ddx2                    &
330                                          + 0.5 * (                           &
331                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
332                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
333                                                  ) * ddy2
334       ENDDO
335
336!
337!--    Apply prescribed horizontal wall heatflux where necessary
338       IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
339       THEN
340          DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
341
342             tend(k,j,i) = tend(k,j,i)                                        &
343                                                + ( fwxp(j,i) * 0.5 *         &
344                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
345                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
346                                                   -fwxm(j,i) * 0.5 *         &
347                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
348                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
349                                                  ) * ddx2                    &
350                                                + ( fwyp(j,i) * 0.5 *         &
351                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
352                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
353                                                   -fwym(j,i) * 0.5 *         &
354                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
355                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
356                                                  ) * ddy2
357          ENDDO
358       ENDIF
359
360!
361!--    Compute vertical diffusion. In case that surface fluxes have been
362!--    prescribed or computed at bottom and/or top, index k starts/ends at
363!--    nzb+2 or nzt-1, respectively.
364       DO  k = nzb_diff_s_inner(j,i), nzt_diff
365
366          tend(k,j,i) = tend(k,j,i)                                           &
367                                       + 0.5 * (                              &
368            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
369          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
370                                               ) * ddzw(k)
371       ENDDO
372
373!
374!--    Vertical diffusion at the first computational gridpoint along z-direction
375       IF ( use_surface_fluxes )  THEN
376
377          k = nzb_s_inner(j,i)+1
378
379          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
380                                            * ( s(k+1,j,i)-s(k,j,i) )    &
381                                            * ddzu(k+1)                  &
382                                        + s_flux_b(j,i)                  &
383                                      ) * ddzw(k)
384
385       ENDIF
386
387!
388!--    Vertical diffusion at the last computational gridpoint along z-direction
389       IF ( use_top_fluxes )  THEN
390
391          k = nzt
392
393          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                  &
394                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
395                                            * ( s(k,j,i)-s(k-1,j,i) )    &
396                                            * ddzu(k)                    &
397                                      ) * ddzw(k)
398
399       ENDIF
400
401    END SUBROUTINE diffusion_s_ij
402
403 END MODULE diffusion_s_mod
Note: See TracBrowser for help on using the repository browser.