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

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

new parameter use_top_fluxes, Bugfix: ddzw dimensioned 1:nzt+1

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