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

Last change on this file since 1691 was 1691, checked in by maronga, 8 years ago

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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