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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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