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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago

Code annotations made doxygen readable

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