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

Last change on this file since 1025 was 1017, checked in by raasch, 12 years ago

last commit documented

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