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

Last change on this file since 44 was 39, checked in by raasch, 18 years ago

comments prepared for 3.1c

  • Property svn:keywords set to Id
File size: 9.2 KB
Line 
1 MODULE diffusion_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_e.f90 39 2007-03-01 12:46:59Z 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                   l_stable = 0.76 * SQRT( e(k,j,i) ) / &
85                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
86                ELSE
87                   l_stable = l_grid(k)
88                ENDIF
89!
90!--             Adjustment of the mixing length
91                IF ( wall_adjustment )  THEN
92                   l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
93                                  l_stable )
94                   ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
95                ELSE
96                   l(k,j)  = MIN( l_grid(k), l_stable )
97                   ll(k,j) = l_grid(k)
98                ENDIF
99                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
100                   l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
101                   ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
102                ENDIF
103
104             ENDDO
105          ENDDO
106!
107!--       Calculate the tendency terms
108          DO  j = nys, nyn
109             DO  k = nzb_s_inner(j,i)+1, nzt
110
111                 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
112                                    e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
113
114                 tend(k,j,i) = tend(k,j,i)                                     &
115                                        + (                                    &
116                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
117                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
118                                          ) * ddx2                             &
119                                        + (                                    &
120                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
121                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
122                                          ) * ddy2                             &
123                                        + (                                    &
124               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
125             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
126                                          ) * ddzw(k)                          &
127                             - dissipation(k,j)
128
129             ENDDO
130          ENDDO
131
132!
133!--       Store dissipation if needed for calculating the sgs particle
134!--       velocities
135          IF ( use_sgs_for_particles )  THEN
136             DO  j = nys, nyn
137                DO  k = nzb_s_inner(j,i)+1, nzt
138                   diss(k,j,i) = dissipation(k,j)
139                ENDDO
140             ENDDO
141          ENDIF
142
143       ENDDO
144
145!
146!--    Boundary condition for dissipation
147       IF ( use_sgs_for_particles )  THEN
148          DO  i = nxl, nxr
149             DO  j = nys, nyn
150                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
151             ENDDO
152          ENDDO
153       ENDIF
154
155    END SUBROUTINE diffusion_e
156
157
158!------------------------------------------------------------------------------!
159! Call for grid point i,j
160!------------------------------------------------------------------------------!
161    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
162                               theta, rif, tend, zu )
163
164       USE control_parameters
165       USE grid_variables
166       USE indices
167       USE particle_attributes
168
169       IMPLICIT NONE
170
171       INTEGER         ::  i, j, k
172       REAL            ::  dpt_dz, l_stable, phi_m
173       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
174                           l_grid(1:nzt), zu(0:nzt+1)
175       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
176       REAL, DIMENSION(:,:), POINTER   ::  rif
177       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
178       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
179
180
181!
182!--    First, calculate phi-function for eventually adjusting the mixing length
183!--    to the prandtl mixing length
184       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
185          IF ( rif(j,i) >= 0.0 )  THEN
186             phi_m = 1.0 + 5.0 * rif(j,i)
187          ELSE
188             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
189          ENDIF
190       ENDIF
191
192!
193!--    Calculate the mixing length (for dissipation)
194       DO  k = nzb_s_inner(j,i)+1, nzt
195          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
196          IF ( dpt_dz > 0.0 ) THEN
197             l_stable = 0.76 * SQRT( e(k,j,i) ) / &
198                               SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
199          ELSE
200             l_stable = l_grid(k)
201          ENDIF
202!
203!--       Adjustment of the mixing length
204          IF ( wall_adjustment )  THEN
205             l(k)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
206             ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
207          ELSE
208             l(k)  = MIN( l_grid(k), l_stable )
209             ll(k) = l_grid(k)
210          ENDIF
211          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
212             l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
213             ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
214          ENDIF
215
216!
217!--       Calculate the tendency term
218          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
219                           SQRT( e(k,j,i) ) / l(k)
220
221          tend(k,j,i) = tend(k,j,i)                                           &
222                                       + (                                    &
223                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
224                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
225                                         ) * ddx2                             &
226                                       + (                                    &
227                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
228                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
229                                         ) * ddy2                             &
230                                       + (                                    &
231              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
232            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
233                                         ) * ddzw(k)                          &
234                                       - dissipation(k)
235
236       ENDDO
237
238!
239!--    Store dissipation if needed for calculating the sgs particle velocities
240       IF ( use_sgs_for_particles )  THEN
241          DO  k = nzb_s_inner(j,i)+1, nzt
242             diss(k,j,i) = dissipation(k)
243          ENDDO
244!
245!--       Boundary condition for dissipation
246          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
247       ENDIF
248
249    END SUBROUTINE diffusion_e_ij
250
251 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.