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

Last change on this file since 19 was 19, checked in by raasch, 17 years ago

preliminary version of modified boundary conditions at top

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