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

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

preliminary version of modified boundary conditions at top

  • Property svn:keywords set to Id
File size: 9.1 KB
Line 
1 MODULE diffusion_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Calculation extended for gridpoint nzt
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_e.f90 19 2007-02-23 04:53:48Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.18  2006/08/04 14:29:43  raasch
14! dissipation is stored in extra array diss if needed later on for calculating
15! the sgs particle velocities
16!
17! Revision 1.1  1997/09/19 07:40:24  raasch
18! Initial revision
19!
20!
21! Description:
22! ------------
23! Diffusion- and dissipation terms for the TKE
24!------------------------------------------------------------------------------!
25
26    PRIVATE
27    PUBLIC diffusion_e
28   
29
30    INTERFACE diffusion_e
31       MODULE PROCEDURE diffusion_e
32       MODULE PROCEDURE diffusion_e_ij
33    END INTERFACE diffusion_e
34 
35 CONTAINS
36
37
38!------------------------------------------------------------------------------!
39! Call for all grid points
40!------------------------------------------------------------------------------!
41    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
42                            rif, tend, zu )
43
44       USE control_parameters
45       USE grid_variables
46       USE indices
47       USE particle_attributes
48
49       IMPLICIT NONE
50
51       INTEGER ::  i, j, k
52       REAL            ::  dpt_dz, l_stable, phi_m
53       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt), &
54                           l_grid(1:nzt), zu(0:nzt+1)
55       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
56       REAL, DIMENSION(:,:), POINTER   ::  rif
57       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
58       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
59 
60
61       DO  i = nxl, nxr
62          DO  j = nys, nyn
63!
64!--          First, calculate phi-function for eventually adjusting the &
65!--          mixing length to the prandtl mixing length
66             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
67                IF ( rif(j,i) >= 0.0 )  THEN
68                   phi_m = 1.0 + 5.0 * rif(j,i)
69                ELSE
70                   phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
71                ENDIF
72             ENDIF
73
74             DO  k = nzb_s_inner(j,i)+1, nzt
75!
76!--             Calculate the mixing length (for dissipation)
77                dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
78                IF ( dpt_dz > 0.0 ) THEN
79                   l_stable = 0.76 * SQRT( e(k,j,i) ) / &
80                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
81                ELSE
82                   l_stable = l_grid(k)
83                ENDIF
84!
85!--             Adjustment of the mixing length
86                IF ( wall_adjustment )  THEN
87                   l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
88                                  l_stable )
89                   ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
90                ELSE
91                   l(k,j)  = MIN( l_grid(k), l_stable )
92                   ll(k,j) = l_grid(k)
93                ENDIF
94                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
95                   l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
96                   ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
97                ENDIF
98
99             ENDDO
100          ENDDO
101!
102!--       Calculate the tendency terms
103          DO  j = nys, nyn
104             DO  k = nzb_s_inner(j,i)+1, nzt
105
106                 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
107                                    e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
108
109                 tend(k,j,i) = tend(k,j,i)                                     &
110                                        + (                                    &
111                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
112                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
113                                          ) * ddx2                             &
114                                        + (                                    &
115                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
116                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
117                                          ) * ddy2                             &
118                                        + (                                    &
119               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
120             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
121                                          ) * ddzw(k)                          &
122                             - dissipation(k,j)
123
124             ENDDO
125          ENDDO
126
127!
128!--       Store dissipation if needed for calculating the sgs particle
129!--       velocities
130          IF ( use_sgs_for_particles )  THEN
131             DO  j = nys, nyn
132                DO  k = nzb_s_inner(j,i)+1, nzt
133                   diss(k,j,i) = dissipation(k,j)
134                ENDDO
135             ENDDO
136          ENDIF
137
138       ENDDO
139
140!
141!--    Boundary condition for dissipation
142       IF ( use_sgs_for_particles )  THEN
143          DO  i = nxl, nxr
144             DO  j = nys, nyn
145                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
146             ENDDO
147          ENDDO
148       ENDIF
149
150    END SUBROUTINE diffusion_e
151
152
153!------------------------------------------------------------------------------!
154! Call for grid point i,j
155!------------------------------------------------------------------------------!
156    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
157                               theta, rif, tend, zu )
158
159       USE control_parameters
160       USE grid_variables
161       USE indices
162       USE particle_attributes
163
164       IMPLICIT NONE
165
166       INTEGER         ::  i, j, k
167       REAL            ::  dpt_dz, l_stable, phi_m
168       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt), &
169                           l_grid(1:nzt), zu(0:nzt+1)
170       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
171       REAL, DIMENSION(:,:), POINTER   ::  rif
172       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
173       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
174
175
176!
177!--    First, calculate phi-function for eventually adjusting the mixing length
178!--    to the prandtl mixing length
179       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
180          IF ( rif(j,i) >= 0.0 )  THEN
181             phi_m = 1.0 + 5.0 * rif(j,i)
182          ELSE
183             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
184          ENDIF
185       ENDIF
186
187!
188!--    Calculate the mixing length (for dissipation)
189       DO  k = nzb_s_inner(j,i)+1, nzt
190          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
191          IF ( dpt_dz > 0.0 ) THEN
192             l_stable = 0.76 * SQRT( e(k,j,i) ) / &
193                               SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
194          ELSE
195             l_stable = l_grid(k)
196          ENDIF
197!
198!--       Adjustment of the mixing length
199          IF ( wall_adjustment )  THEN
200             l(k)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
201             ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
202          ELSE
203             l(k)  = MIN( l_grid(k), l_stable )
204             ll(k) = l_grid(k)
205          ENDIF
206          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
207             l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
208             ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
209          ENDIF
210
211!
212!--       Calculate the tendency term
213          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
214                           SQRT( e(k,j,i) ) / l(k)
215
216          tend(k,j,i) = tend(k,j,i)                                           &
217                                       + (                                    &
218                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
219                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
220                                         ) * ddx2                             &
221                                       + (                                    &
222                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
223                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
224                                         ) * ddy2                             &
225                                       + (                                    &
226              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
227            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
228                                         ) * ddzw(k)                          &
229                                       - dissipation(k)
230
231       ENDDO
232
233!
234!--    Store dissipation if needed for calculating the sgs particle velocities
235       IF ( use_sgs_for_particles )  THEN
236          DO  k = nzb_s_inner(j,i)+1, nzt
237             diss(k,j,i) = dissipation(k)
238          ENDDO
239!
240!--       Boundary condition for dissipation
241          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
242       ENDIF
243
244    END SUBROUTINE diffusion_e_ij
245
246 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.