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

Last change on this file since 1660 was 1375, checked in by raasch, 11 years ago

last commit documented

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