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

Last change on this file since 1108 was 1093, checked in by raasch, 12 years ago

last commit documented

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