Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 76 was 65, checked in by raasch, 17 years ago

bugfix for NEC because -C hopt gives runtime errors
  • Property svn:keywords set to Id
File size: 13.2 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 65 2007-03-13 12:11:43Z 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!
67!--    This if clause must be outside the k-loop because otherwise
68!--    runtime errors occur with -C hopt on NEC
69       IF ( use_pt_reference )  THEN
70
71          DO  i = nxl, nxr
72             DO  j = nys, nyn
73!
74!--             First, calculate phi-function for eventually adjusting the &
75!--             mixing length to the prandtl mixing length
76                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
77                   IF ( rif(j,i) >= 0.0 )  THEN
78                      phi_m = 1.0 + 5.0 * rif(j,i)
79                   ELSE
80                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
81                   ENDIF
82                ENDIF
83
84                DO  k = nzb_s_inner(j,i)+1, nzt
85!
86!--                Calculate the mixing length (for dissipation)
87                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
88                   IF ( dpt_dz > 0.0 ) THEN
89                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
90                                        SQRT( g / pt_reference * dpt_dz ) + 1E-5
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!
113!--          Calculate the tendency terms
114             DO  j = nys, nyn
115                DO  k = nzb_s_inner(j,i)+1, nzt
116
117                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
118                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
119
120                    tend(k,j,i) = tend(k,j,i)                                  &
121                                        + (                                    &
122                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
123                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
124                                          ) * ddx2                             &
125                                        + (                                    &
126                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
127                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
128                                          ) * ddy2                             &
129                                        + (                                    &
130               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
131             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
132                                          ) * ddzw(k)                          &
133                             - dissipation(k,j)
134
135                ENDDO
136             ENDDO
137
138!
139!--          Store dissipation if needed for calculating the sgs particle
140!--          velocities
141             IF ( use_sgs_for_particles )  THEN
142                DO  j = nys, nyn
143                   DO  k = nzb_s_inner(j,i)+1, nzt
144                      diss(k,j,i) = dissipation(k,j)
145                   ENDDO
146                ENDDO
147             ENDIF
148
149          ENDDO
150
151       ELSE
152
153          DO  i = nxl, nxr
154             DO  j = nys, nyn
155!
156!--             First, calculate phi-function for eventually adjusting the &
157!--             mixing length to the prandtl mixing length
158                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
159                   IF ( rif(j,i) >= 0.0 )  THEN
160                      phi_m = 1.0 + 5.0 * rif(j,i)
161                   ELSE
162                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
163                   ENDIF
164                ENDIF
165
166                DO  k = nzb_s_inner(j,i)+1, nzt
167!
168!--                Calculate the mixing length (for dissipation)
169                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
170                   IF ( dpt_dz > 0.0 ) THEN
171                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
172                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
173                   ELSE
174                      l_stable = l_grid(k)
175                   ENDIF
176!
177!--                Adjustment of the mixing length
178                   IF ( wall_adjustment )  THEN
179                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
180                                     l_stable )
181                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
182                   ELSE
183                      l(k,j)  = MIN( l_grid(k), l_stable )
184                      ll(k,j) = l_grid(k)
185                   ENDIF
186                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
187                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
188                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
189                   ENDIF
190
191                ENDDO
192             ENDDO
193
194!
195!--          Calculate the tendency terms
196             DO  j = nys, nyn
197                DO  k = nzb_s_inner(j,i)+1, nzt
198
199                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
200                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
201
202                    tend(k,j,i) = tend(k,j,i)                                  &
203                                        + (                                    &
204                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
205                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
206                                          ) * ddx2                             &
207                                        + (                                    &
208                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
209                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
210                                          ) * ddy2                             &
211                                        + (                                    &
212               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
213             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
214                                          ) * ddzw(k)                          &
215                             - dissipation(k,j)
216
217                ENDDO
218             ENDDO
219
220!
221!--          Store dissipation if needed for calculating the sgs particle
222!--          velocities
223             IF ( use_sgs_for_particles )  THEN
224                DO  j = nys, nyn
225                   DO  k = nzb_s_inner(j,i)+1, nzt
226                      diss(k,j,i) = dissipation(k,j)
227                   ENDDO
228                ENDDO
229             ENDIF
230
231          ENDDO
232
233       ENDIF
234
235!
236!--    Boundary condition for dissipation
237       IF ( use_sgs_for_particles )  THEN
238          DO  i = nxl, nxr
239             DO  j = nys, nyn
240                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
241             ENDDO
242          ENDDO
243       ENDIF
244
245    END SUBROUTINE diffusion_e
246
247
248!------------------------------------------------------------------------------!
249! Call for grid point i,j
250!------------------------------------------------------------------------------!
251    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
252                               theta, rif, tend, zu )
253
254       USE control_parameters
255       USE grid_variables
256       USE indices
257       USE particle_attributes
258
259       IMPLICIT NONE
260
261       INTEGER         ::  i, j, k
262       REAL            ::  dpt_dz, l_stable, phi_m
263       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
264                           l_grid(1:nzt), zu(0:nzt+1)
265       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
266       REAL, DIMENSION(:,:), POINTER   ::  rif
267       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
268       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
269
270
271!
272!--    First, calculate phi-function for eventually adjusting the mixing length
273!--    to the prandtl mixing length
274       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
275          IF ( rif(j,i) >= 0.0 )  THEN
276             phi_m = 1.0 + 5.0 * rif(j,i)
277          ELSE
278             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
279          ENDIF
280       ENDIF
281
282!
283!--    Calculate the mixing length (for dissipation)
284       DO  k = nzb_s_inner(j,i)+1, nzt
285          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
286          IF ( dpt_dz > 0.0 ) THEN
287             IF ( use_pt_reference )  THEN
288                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
289                                  SQRT( g / pt_reference * dpt_dz ) + 1E-5
290             ELSE
291                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
292                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
293             ENDIF
294          ELSE
295             l_stable = l_grid(k)
296          ENDIF
297!
298!--       Adjustment of the mixing length
299          IF ( wall_adjustment )  THEN
300             l(k)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
301             ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
302          ELSE
303             l(k)  = MIN( l_grid(k), l_stable )
304             ll(k) = l_grid(k)
305          ENDIF
306          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
307             l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
308             ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
309          ENDIF
310
311!
312!--       Calculate the tendency term
313          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
314                           SQRT( e(k,j,i) ) / l(k)
315
316          tend(k,j,i) = tend(k,j,i)                                           &
317                                       + (                                    &
318                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
319                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
320                                         ) * ddx2                             &
321                                       + (                                    &
322                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
323                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
324                                         ) * ddy2                             &
325                                       + (                                    &
326              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
327            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
328                                         ) * ddzw(k)                          &
329                                       - dissipation(k)
330
331       ENDDO
332
333!
334!--    Store dissipation if needed for calculating the sgs particle velocities
335       IF ( use_sgs_for_particles )  THEN
336          DO  k = nzb_s_inner(j,i)+1, nzt
337             diss(k,j,i) = dissipation(k)
338          ENDDO
339!
340!--       Boundary condition for dissipation
341          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
342       ENDIF
343
344    END SUBROUTINE diffusion_e_ij
345
346 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |