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

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