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

Last change on this file since 449 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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