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

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

leapfrog timestep scheme and upstream-spline advection scheme completely removed from the code,
reading of dt_fixed from restart file removed

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