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

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

code has been put under the GNU General Public License (v3)

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