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

Last change on this file since 183 was 183, checked in by letzel, 13 years ago
  • mrun: (optionally) link local input
  • diffusion_s: bugfix for calculation of horizontal fluxes at vertical walls
  • NCL scripts in trunk/SCRIPTS/NCL updated
  • Property svn:keywords set to Id
File size: 10.8 KB
Line 
1 MODULE diffusion_s_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! bugfix: calculation of fluxes at vertical surfaces
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_s.f90 183 2008-08-04 15:39:12Z letzel $
11!
12! 129 2007-10-30 12:12:24Z letzel
13! replace wall_heatflux by wall_s_flux that is now included in the parameter
14! list, bugfix for assignment of fluxes at walls
15!
16! 20 2007-02-26 00:12:32Z raasch
17! Bugfix: ddzw dimensioned 1:nzt"+1"
18! Calculation extended for gridpoint nzt, fluxes can be given at top,
19! +s_flux_t in parameter list, s_flux renamed s_flux_b
20!
21! RCS Log replace by Id keyword, revision history cleaned up
22!
23! Revision 1.8  2006/02/23 10:34:17  raasch
24! nzb_2d replaced by nzb_s_outer in horizontal diffusion and by nzb_s_inner
25! or nzb_diff_s_inner, respectively, in vertical diffusion, prescribed surface
26! fluxes at vertically oriented topography
27!
28! Revision 1.1  2000/04/13 14:54:02  schroeter
29! Initial revision
30!
31!
32! Description:
33! ------------
34! Diffusion term of scalar quantities (temperature and water content)
35!------------------------------------------------------------------------------!
36
37    PRIVATE
38    PUBLIC diffusion_s
39
40    INTERFACE diffusion_s
41       MODULE PROCEDURE diffusion_s
42       MODULE PROCEDURE diffusion_s_ij
43    END INTERFACE diffusion_s
44
45 CONTAINS
46
47
48!------------------------------------------------------------------------------!
49! Call for all grid points
50!------------------------------------------------------------------------------!
51    SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux_b, s_flux_t, &
52                            wall_s_flux, tend )
53
54       USE control_parameters
55       USE grid_variables
56       USE indices
57
58       IMPLICIT NONE
59
60       INTEGER ::  i, j, k
61       REAL    ::  vertical_gridspace
62       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
63       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
64       REAL    ::  wall_s_flux(0:4)
65       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
66       REAL, DIMENSION(:,:,:), POINTER ::  kh, s
67
68       DO  i = nxl, nxr
69          DO  j = nys,nyn
70!
71!--          Compute horizontal diffusion
72             DO  k = nzb_s_outer(j,i)+1, nzt
73
74                tend(k,j,i) = tend(k,j,i)                                     &
75                                          + 0.5 * (                           &
76                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
77                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
78                                                  ) * ddx2                    &
79                                          + 0.5 * (                           &
80                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
81                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
82                                                  ) * ddy2
83             ENDDO
84
85!
86!--          Apply prescribed horizontal wall heatflux where necessary
87             IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
88             THEN
89                DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
90
91                   tend(k,j,i) = tend(k,j,i)                                  &
92                                                + ( fwxp(j,i) * 0.5 *         &
93                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
94                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
95                                                   -fwxm(j,i) * 0.5 *         &
96                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
97                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
98                                                  ) * ddx2                    &
99                                                + ( fwyp(j,i) * 0.5 *         &
100                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
101                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
102                                                   -fwym(j,i) * 0.5 *         &
103                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
104                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
105                                                  ) * ddy2
106                ENDDO
107             ENDIF
108
109!
110!--          Compute vertical diffusion. In case that surface fluxes have been
111!--          prescribed or computed at bottom and/or top, index k starts/ends at
112!--          nzb+2 or nzt-1, respectively.
113             DO  k = nzb_diff_s_inner(j,i), nzt_diff
114
115                tend(k,j,i) = tend(k,j,i)                                     &
116                                       + 0.5 * (                              &
117            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
118          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
119                                               ) * ddzw(k)
120             ENDDO
121
122!
123!--          Vertical diffusion at the first computational gridpoint along
124!--          z-direction
125             IF ( use_surface_fluxes )  THEN
126
127                k = nzb_s_inner(j,i)+1
128
129                tend(k,j,i) = tend(k,j,i)                                     &
130                                       + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )    &
131                                               * ( s(k+1,j,i)-s(k,j,i) )      &
132                                               * ddzu(k+1)                    &
133                                           + s_flux_b(j,i)                    &
134                                         ) * ddzw(k)
135
136             ENDIF
137
138!
139!--          Vertical diffusion at the last computational gridpoint along
140!--          z-direction
141             IF ( use_top_fluxes )  THEN
142
143                k = nzt
144
145                tend(k,j,i) = tend(k,j,i)                                     &
146                                       + ( - s_flux_t(j,i)                    &
147                                           - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
148                                                 * ( s(k,j,i)-s(k-1,j,i) )    &
149                                                 * ddzu(k)                    &
150                                         ) * ddzw(k)
151
152             ENDIF
153
154          ENDDO
155       ENDDO
156
157    END SUBROUTINE diffusion_s
158
159
160!------------------------------------------------------------------------------!
161! Call for grid point i,j
162!------------------------------------------------------------------------------!
163    SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux_b, s_flux_t, &
164                               wall_s_flux, tend )
165
166       USE control_parameters
167       USE grid_variables
168       USE indices
169
170       IMPLICIT NONE
171
172       INTEGER ::  i, j, k
173       REAL    ::  vertical_gridspace
174       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
175       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
176       REAL    ::  wall_s_flux(0:4)
177       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
178       REAL, DIMENSION(:,:,:), POINTER ::  kh, s
179
180!
181!--    Compute horizontal diffusion
182       DO  k = nzb_s_outer(j,i)+1, nzt
183
184          tend(k,j,i) = tend(k,j,i)                                           &
185                                          + 0.5 * (                           &
186                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
187                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
188                                                  ) * ddx2                    &
189                                          + 0.5 * (                           &
190                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
191                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
192                                                  ) * ddy2
193       ENDDO
194
195!
196!--    Apply prescribed horizontal wall heatflux where necessary
197       IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
198       THEN
199          DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
200
201             tend(k,j,i) = tend(k,j,i)                                        &
202                                                + ( fwxp(j,i) * 0.5 *         &
203                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
204                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
205                                                   -fwxm(j,i) * 0.5 *         &
206                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
207                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
208                                                  ) * ddx2                    &
209                                                + ( fwyp(j,i) * 0.5 *         &
210                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
211                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
212                                                   -fwym(j,i) * 0.5 *         &
213                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
214                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
215                                                  ) * ddy2
216          ENDDO
217       ENDIF
218
219!
220!--    Compute vertical diffusion. In case that surface fluxes have been
221!--    prescribed or computed at bottom and/or top, index k starts/ends at
222!--    nzb+2 or nzt-1, respectively.
223       DO  k = nzb_diff_s_inner(j,i), nzt_diff
224
225          tend(k,j,i) = tend(k,j,i)                                           &
226                                       + 0.5 * (                              &
227            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
228          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
229                                               ) * ddzw(k)
230       ENDDO
231
232!
233!--    Vertical diffusion at the first computational gridpoint along z-direction
234       IF ( use_surface_fluxes )  THEN
235
236          k = nzb_s_inner(j,i)+1
237
238          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
239                                            * ( s(k+1,j,i)-s(k,j,i) )    &
240                                            * ddzu(k+1)                  &
241                                        + s_flux_b(j,i)                  &
242                                      ) * ddzw(k)
243
244       ENDIF
245
246!
247!--    Vertical diffusion at the last computational gridpoint along z-direction
248       IF ( use_top_fluxes )  THEN
249
250          k = nzt
251
252          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                  &
253                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
254                                            * ( s(k,j,i)-s(k-1,j,i) )    &
255                                            * ddzu(k)                    &
256                                      ) * ddzw(k)
257
258       ENDIF
259
260    END SUBROUTINE diffusion_s_ij
261
262 END MODULE diffusion_s_mod
Note: See TracBrowser for help on using the repository browser.