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

Last change on this file since 1310 was 1310, checked in by raasch, 10 years ago

update of GPL copyright

  • Property svn:keywords set to Id
File size: 17.7 KB
RevLine 
[1]1 MODULE diffusion_s_mod
2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[484]20! Current revisions:
[1001]21! ------------------
[1]22!
[1258]23!
[1]24! Former revisions:
25! -----------------
[3]26! $Id: diffusion_s.f90 1310 2014-03-14 08:01:56Z raasch $
[39]27!
[1258]28! 1257 2013-11-08 15:18:40Z raasch
29! openacc loop and loop vector clauses removed
30!
[1132]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!
[1093]35! 1092 2013-02-02 11:24:22Z raasch
36! unused variables removed
37!
[1037]38! 1036 2012-10-22 13:43:42Z raasch
39! code put under GPL (PALM 3.9)
40!
[1017]41! 1015 2012-09-27 09:23:24Z raasch
42! accelerator version (*_acc) added
43!
[1011]44! 1010 2012-09-20 07:59:54Z raasch
45! cpp switch __nopointer added for pointer free version
46!
[1002]47! 1001 2012-09-13 14:08:46Z raasch
48! some arrays comunicated by module instead of parameter list
49!
[668]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!
[198]53! 183 2008-08-04 15:39:12Z letzel
54! bugfix: calculation of fluxes at vertical surfaces
55!
[139]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!
[39]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!
[3]65! RCS Log replace by Id keyword, revision history cleaned up
66!
[1]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
[1015]82    PUBLIC diffusion_s, diffusion_s_acc
[1]83
84    INTERFACE diffusion_s
85       MODULE PROCEDURE diffusion_s
86       MODULE PROCEDURE diffusion_s_ij
87    END INTERFACE diffusion_s
88
[1015]89    INTERFACE diffusion_s_acc
90       MODULE PROCEDURE diffusion_s_acc
91    END INTERFACE diffusion_s_acc
92
[1]93 CONTAINS
94
95
96!------------------------------------------------------------------------------!
97! Call for all grid points
98!------------------------------------------------------------------------------!
[1001]99    SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
[1]100
[1001]101       USE arrays_3d
[1]102       USE control_parameters
103       USE grid_variables
104       USE indices
105
106       IMPLICIT NONE
107
108       INTEGER ::  i, j, k
[129]109       REAL    ::  wall_s_flux(0:4)
[1001]110       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
[1010]111#if defined( __nopointer )
112       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
113#else
[1001]114       REAL, DIMENSION(:,:,:), POINTER ::  s
[1010]115#endif
[1]116
117       DO  i = nxl, nxr
118          DO  j = nys,nyn
119!
120!--          Compute horizontal diffusion
[19]121             DO  k = nzb_s_outer(j,i)+1, nzt
[1]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)                                  &
[183]141                                                + ( fwxp(j,i) * 0.5 *         &
[1]142                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
[129]143                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
[183]144                                                   -fwxm(j,i) * 0.5 *         &
[1]145                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
[129]146                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
[1]147                                                  ) * ddx2                    &
[183]148                                                + ( fwyp(j,i) * 0.5 *         &
[1]149                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
[129]150                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
[183]151                                                   -fwym(j,i) * 0.5 *         &
[1]152                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
[129]153                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
[1]154                                                  ) * ddy2
155                ENDDO
156             ENDIF
157
158!
159!--          Compute vertical diffusion. In case that surface fluxes have been
[19]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
[1]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!
[19]172!--          Vertical diffusion at the first computational gridpoint along
[1]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)                    &
[19]182                                           + s_flux_b(j,i)                    &
[1]183                                         ) * ddzw(k)
184
185             ENDIF
186
[19]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)                    &
[20]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)                    &
[19]199                                         ) * ddzw(k)
200
201             ENDIF
202
[1]203          ENDDO
204       ENDDO
205
206    END SUBROUTINE diffusion_s
207
208
209!------------------------------------------------------------------------------!
[1015]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 )
[1128]234       DO  i = i_left, i_right
235          DO  j = j_south, j_north
[1015]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!------------------------------------------------------------------------------!
[1]325! Call for grid point i,j
326!------------------------------------------------------------------------------!
[1001]327    SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
[1]328
[1001]329       USE arrays_3d
[1]330       USE control_parameters
331       USE grid_variables
332       USE indices
333
334       IMPLICIT NONE
335
336       INTEGER ::  i, j, k
[129]337       REAL    ::  wall_s_flux(0:4)
[1001]338       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
[1010]339#if defined( __nopointer )
340       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
341#else
[1001]342       REAL, DIMENSION(:,:,:), POINTER ::  s
[1010]343#endif
[1]344
345!
346!--    Compute horizontal diffusion
[19]347       DO  k = nzb_s_outer(j,i)+1, nzt
[1]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)                                        &
[183]367                                                + ( fwxp(j,i) * 0.5 *         &
[1]368                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
[129]369                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
[183]370                                                   -fwxm(j,i) * 0.5 *         &
[1]371                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
[129]372                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
[1]373                                                  ) * ddx2                    &
[183]374                                                + ( fwyp(j,i) * 0.5 *         &
[1]375                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
[129]376                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
[183]377                                                   -fwym(j,i) * 0.5 *         &
[1]378                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
[129]379                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
[1]380                                                  ) * ddy2
381          ENDDO
382       ENDIF
383
384!
385!--    Compute vertical diffusion. In case that surface fluxes have been
[19]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
[1]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!
[19]398!--    Vertical diffusion at the first computational gridpoint along z-direction
[1]399       IF ( use_surface_fluxes )  THEN
400
401          k = nzb_s_inner(j,i)+1
402
[19]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)
[1]408
409       ENDIF
410
[19]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
[1]425    END SUBROUTINE diffusion_s_ij
426
427 END MODULE diffusion_s_mod
Note: See TracBrowser for help on using the repository browser.