source: palm/trunk/SOURCE/diffusion_w.f90 @ 1004

Last change on this file since 1004 was 1002, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1 MODULE diffusion_w_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_w.f90 1002 2012-09-13 15:12:24Z raasch $
11!
12! 1001 2012-09-13 14:08:46Z raasch
13! arrays comunicated by module instead of parameter list
14!
15! 978 2012-08-09 08:28:32Z fricke
16! outflow damping layer removed
17! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
18! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
19!
20! 667 2010-12-23 12:06:00Z suehring/gryschka
21! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
22!
23! 366 2009-08-25 08:06:27Z raasch
24! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
25!
26! 75 2007-03-22 09:54:05Z raasch
27! Wall functions now include diabatic conditions, call of routine wall_fluxes,
28! z0 removed from argument list
29!
30! 20 2007-02-26 00:12:32Z raasch
31! Bugfix: ddzw dimensioned 1:nzt"+1"
32!
33! RCS Log replace by Id keyword, revision history cleaned up
34!
35! Revision 1.12  2006/02/23 10:38:03  raasch
36! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
37! +z0 in argument list
38! WARNING: loops containing the MAX function are still not properly vectorized!
39!
40! Revision 1.1  1997/09/12 06:24:11  raasch
41! Initial revision
42!
43!
44! Description:
45! ------------
46! Diffusion term of the w-component
47!------------------------------------------------------------------------------!
48
49    USE wall_fluxes_mod
50
51    PRIVATE
52    PUBLIC diffusion_w
53
54    INTERFACE diffusion_w
55       MODULE PROCEDURE diffusion_w
56       MODULE PROCEDURE diffusion_w_ij
57    END INTERFACE diffusion_w
58
59 CONTAINS
60
61
62!------------------------------------------------------------------------------!
63! Call for all grid points
64!------------------------------------------------------------------------------!
65    SUBROUTINE diffusion_w
66
67       USE arrays_3d
68       USE control_parameters
69       USE grid_variables
70       USE indices
71
72       IMPLICIT NONE
73
74       INTEGER ::  i, j, k
75       REAL    ::  kmxm, kmxp, kmym, kmyp
76
77       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
78
79
80!
81!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
82!--    walls, if neccessary
83       IF ( topography /= 'flat' )  THEN
84          CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
85                            nzb_w_outer, wall_w_x )
86          CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
87                            nzb_w_outer, wall_w_y )
88       ENDIF
89
90       DO  i = nxl, nxr
91          DO  j = nys, nyn
92             DO  k = nzb_w_outer(j,i)+1, nzt-1
93!
94!--             Interpolate eddy diffusivities on staggered gridpoints
95                kmxp = 0.25 * &
96                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
97                kmxm = 0.25 * &
98                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
99                kmyp = 0.25 * &
100                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
101                kmym = 0.25 * &
102                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
103
104                tend(k,j,i) = tend(k,j,i)                                      &
105                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
106                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
107                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
108                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
109                      &   ) * ddx                                              &
110                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
111                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
112                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
113                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
114                      &   ) * ddy                                              &
115                      & + 2.0 * (                                              &
116                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
117                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
118                      &         ) * ddzu(k+1)
119             ENDDO
120
121!
122!--          Wall functions at all vertical walls, where necessary
123             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
124
125                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
126!
127!--                Interpolate eddy diffusivities on staggered gridpoints
128                   kmxp = 0.25 * &
129                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
130                   kmxm = 0.25 * &
131                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
132                   kmyp = 0.25 * &
133                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
134                   kmym = 0.25 * &
135                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
136
137                   tend(k,j,i) = tend(k,j,i)                                   &
138                                 + (   fwxp(j,i) * (                           &
139                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
140                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
141                                                   )                           &
142                                     - fwxm(j,i) * (                           &
143                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
144                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
145                                                   )                           &
146                                     + wall_w_x(j,i) * wsus(k,j,i)             &
147                                   ) * ddx                                     &
148                                 + (   fwyp(j,i) * (                           &
149                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
150                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
151                                                   )                           &
152                                     - fwym(j,i) * (                           &
153                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
154                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
155                                                   )                           &
156                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
157                                   ) * ddy                                     &
158                                 + 2.0 * (                                     &
159                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
160                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
161                                         ) * ddzu(k+1)
162                ENDDO
163             ENDIF
164
165          ENDDO
166       ENDDO
167
168    END SUBROUTINE diffusion_w
169
170
171!------------------------------------------------------------------------------!
172! Call for grid point i,j
173!------------------------------------------------------------------------------!
174    SUBROUTINE diffusion_w_ij( i, j )
175
176       USE arrays_3d
177       USE control_parameters
178       USE grid_variables
179       USE indices
180
181       IMPLICIT NONE
182
183       INTEGER ::  i, j, k
184       REAL    ::  kmxm, kmxp, kmym, kmyp
185
186       REAL, DIMENSION(nzb:nzt+1) ::  wsus, wsvs
187
188
189       DO  k = nzb_w_outer(j,i)+1, nzt-1
190!
191!--       Interpolate eddy diffusivities on staggered gridpoints
192          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
193          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
194          kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
195          kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
196
197          tend(k,j,i) = tend(k,j,i)                                            &
198                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
199                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
200                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
201                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
202                      &   ) * ddx                                              &
203                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
204                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
205                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
206                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
207                      &   ) * ddy                                              &
208                      & + 2.0 * (                                              &
209                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
210                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
211                      &         ) * ddzu(k+1)
212       ENDDO
213
214!
215!--    Wall functions at all vertical walls, where necessary
216       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
217
218!
219!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
220          IF ( wall_w_x(j,i) /= 0.0 )  THEN
221             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
222                               wsus, 0.0, 0.0, 0.0, 1.0 )
223          ELSE
224             wsus = 0.0
225          ENDIF
226
227          IF ( wall_w_y(j,i) /= 0.0 )  THEN
228             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),  &
229                               wsvs, 0.0, 0.0, 1.0, 0.0 )
230          ELSE
231             wsvs = 0.0
232          ENDIF
233
234          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
235!
236!--          Interpolate eddy diffusivities on staggered gridpoints
237             kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
238             kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
239             kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
240             kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
241
242             tend(k,j,i) = tend(k,j,i)                                         &
243                                 + (   fwxp(j,i) * (                           &
244                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
245                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
246                                                   )                           &
247                                     - fwxm(j,i) * (                           &
248                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
249                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
250                                                   )                           &
251                                     + wall_w_x(j,i) * wsus(k)                 &
252                                   ) * ddx                                     &
253                                 + (   fwyp(j,i) * (                           &
254                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
255                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
256                                                   )                           &
257                                     - fwym(j,i) * (                           &
258                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
259                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
260                                                   )                           &
261                                     + wall_w_y(j,i) * wsvs(k)                 &
262                                   ) * ddy                                     &
263                                 + 2.0 * (                                     &
264                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
265                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
266                                         ) * ddzu(k+1)
267          ENDDO
268       ENDIF
269
270    END SUBROUTINE diffusion_w_ij
271
272 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.