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

Last change on this file since 856 was 668, checked in by suehring, 14 years ago

last commit documented

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