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

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