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

Last change on this file since 1132 was 1132, checked in by raasch, 8 years ago

r1028 documented

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