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

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

preliminary update of further changes, advec_particles is not running!

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