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

Last change on this file since 888 was 826, checked in by raasch, 13 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 15.4 KB
Line 
1 MODULE diffusion_e_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_e.f90 826 2012-02-19 03:41:34Z suehring $
11!
12! 825 2012-02-19 03:03:44Z raasch
13! wang_collision_kernel renamed wang_kernel
14!
15! 790 2011-11-29 03:11:20Z raasch
16! diss is also calculated in case that the Wang kernel is used
17!
18! 667 2010-12-23 12:06:00Z suehring/gryschka
19! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
20!
21! 97 2007-06-21 08:23:15Z raasch
22! Adjustment of mixing length calculation for the ocean version. zw added to
23! argument list.
24! This is also a bugfix, because the height above the topography is now
25! used instead of the height above level k=0.
26! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
27! use_pt_reference renamed use_reference
28!
29! 65 2007-03-13 12:11:43Z raasch
30! Reference temperature pt_reference can be used in buoyancy term
31!
32! 20 2007-02-26 00:12:32Z raasch
33! Bugfix: ddzw dimensioned 1:nzt"+1"
34! Calculation extended for gridpoint nzt
35!
36! RCS Log replace by Id keyword, revision history cleaned up
37!
38! Revision 1.18  2006/08/04 14:29:43  raasch
39! dissipation is stored in extra array diss if needed later on for calculating
40! the sgs particle velocities
41!
42! Revision 1.1  1997/09/19 07:40:24  raasch
43! Initial revision
44!
45!
46! Description:
47! ------------
48! Diffusion- and dissipation terms for the TKE
49!------------------------------------------------------------------------------!
50
51    PRIVATE
52    PUBLIC diffusion_e
53   
54
55    INTERFACE diffusion_e
56       MODULE PROCEDURE diffusion_e
57       MODULE PROCEDURE diffusion_e_ij
58    END INTERFACE diffusion_e
59 
60 CONTAINS
61
62
63!------------------------------------------------------------------------------!
64! Call for all grid points
65!------------------------------------------------------------------------------!
66    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, var, &
67                            var_reference, rif, tend, zu, zw )
68
69       USE control_parameters
70       USE grid_variables
71       USE indices
72       USE particle_attributes
73
74       IMPLICIT NONE
75
76       INTEGER ::  i, j, k
77       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
78       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
79                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
80       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
81       REAL, DIMENSION(:,:), POINTER   ::  rif
82       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
83       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
84 
85
86!
87!--    This if clause must be outside the k-loop because otherwise
88!--    runtime errors occur with -C hopt on NEC
89       IF ( use_reference )  THEN
90
91          DO  i = nxl, nxr
92             DO  j = nys, nyn
93!
94!--             First, calculate phi-function for eventually adjusting the &
95!--             mixing length to the prandtl mixing length
96                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
97                   IF ( rif(j,i) >= 0.0 )  THEN
98                      phi_m = 1.0 + 5.0 * rif(j,i)
99                   ELSE
100                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
101                   ENDIF
102                ENDIF
103
104                DO  k = nzb_s_inner(j,i)+1, nzt
105!
106!--                Calculate the mixing length (for dissipation)
107                   dvar_dz = atmos_ocean_sign * &
108                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
109                   IF ( dvar_dz > 0.0 ) THEN
110                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
111                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
112                   ELSE
113                      l_stable = l_grid(k)
114                   ENDIF
115!
116!--                Adjustment of the mixing length
117                   IF ( wall_adjustment )  THEN
118                      l(k,j)  = MIN( wall_adjustment_factor *          &
119                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
120                                     l_grid(k), l_stable )
121                      ll(k,j) = MIN( wall_adjustment_factor *          &
122                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
123                                     l_grid(k) )
124                   ELSE
125                      l(k,j)  = MIN( l_grid(k), l_stable )
126                      ll(k,j) = l_grid(k)
127                   ENDIF
128                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
129                      l(k,j)  = MIN( l(k,j),  kappa *                          &
130                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
131                                              / phi_m )
132                      ll(k,j) = MIN( ll(k,j), kappa *                          &
133                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
134                                              / phi_m )
135                   ENDIF
136
137                ENDDO
138             ENDDO
139
140!
141!--          Calculate the tendency terms
142             DO  j = nys, nyn
143                DO  k = nzb_s_inner(j,i)+1, nzt
144
145                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
146                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
147
148                    tend(k,j,i) = tend(k,j,i)                                  &
149                                        + (                                    &
150                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
151                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
152                                          ) * ddx2                             &
153                                        + (                                    &
154                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
155                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
156                                          ) * ddy2                             &
157                                        + (                                    &
158               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
159             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
160                                          ) * ddzw(k)                          &
161                             - dissipation(k,j)
162
163                ENDDO
164             ENDDO
165
166!
167!--          Store dissipation if needed for calculating the sgs particle
168!--          velocities
169             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
170                DO  j = nys, nyn
171                   DO  k = nzb_s_inner(j,i)+1, nzt
172                      diss(k,j,i) = dissipation(k,j)
173                   ENDDO
174                ENDDO
175             ENDIF
176
177          ENDDO
178
179       ELSE
180
181          DO  i = nxl, nxr
182             DO  j = nys, nyn
183!
184!--             First, calculate phi-function for eventually adjusting the &
185!--             mixing length to the prandtl mixing length
186                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
187                   IF ( rif(j,i) >= 0.0 )  THEN
188                      phi_m = 1.0 + 5.0 * rif(j,i)
189                   ELSE
190                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
191                   ENDIF
192                ENDIF
193
194                DO  k = nzb_s_inner(j,i)+1, nzt
195!
196!--                Calculate the mixing length (for dissipation)
197                   dvar_dz = atmos_ocean_sign * &
198                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
199                   IF ( dvar_dz > 0.0 ) THEN
200                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
201                                        SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
202                   ELSE
203                      l_stable = l_grid(k)
204                   ENDIF
205!
206!--                Adjustment of the mixing length
207                   IF ( wall_adjustment )  THEN
208                      l(k,j)  = MIN( wall_adjustment_factor *          &
209                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
210                                     l_grid(k), l_stable )
211                      ll(k,j) = MIN( wall_adjustment_factor *          &
212                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
213                                     l_grid(k) )
214                   ELSE
215                      l(k,j)  = MIN( l_grid(k), l_stable )
216                      ll(k,j) = l_grid(k)
217                   ENDIF
218                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
219                      l(k,j)  = MIN( l(k,j),  kappa *                          &
220                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
221                                              / phi_m )
222                      ll(k,j) = MIN( ll(k,j), kappa *                          &
223                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
224                                              / phi_m )
225                   ENDIF
226
227                ENDDO
228             ENDDO
229
230!
231!--          Calculate the tendency terms
232             DO  j = nys, nyn
233                DO  k = nzb_s_inner(j,i)+1, nzt
234
235                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
236                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
237
238                    tend(k,j,i) = tend(k,j,i)                                  &
239                                        + (                                    &
240                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
241                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
242                                          ) * ddx2                             &
243                                        + (                                    &
244                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
245                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
246                                          ) * ddy2                             &
247                                        + (                                    &
248               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
249             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
250                                          ) * ddzw(k)                          &
251                             - dissipation(k,j)
252
253                ENDDO
254             ENDDO
255
256!
257!--          Store dissipation if needed for calculating the sgs particle
258!--          velocities
259             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
260                DO  j = nys, nyn
261                   DO  k = nzb_s_inner(j,i)+1, nzt
262                      diss(k,j,i) = dissipation(k,j)
263                   ENDDO
264                ENDDO
265             ENDIF
266
267          ENDDO
268
269       ENDIF
270
271!
272!--    Boundary condition for dissipation
273       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
274          DO  i = nxl, nxr
275             DO  j = nys, nyn
276                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
277             ENDDO
278          ENDDO
279       ENDIF
280
281    END SUBROUTINE diffusion_e
282
283
284!------------------------------------------------------------------------------!
285! Call for grid point i,j
286!------------------------------------------------------------------------------!
287    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
288                               var, var_reference, rif, tend, zu, zw )
289
290       USE control_parameters
291       USE grid_variables
292       USE indices
293       USE particle_attributes
294
295       IMPLICIT NONE
296
297       INTEGER         ::  i, j, k
298       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
299       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
300                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
301       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
302       REAL, DIMENSION(:,:), POINTER   ::  rif
303       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
304       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
305
306
307!
308!--    First, calculate phi-function for eventually adjusting the mixing length
309!--    to the prandtl mixing length
310       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
311          IF ( rif(j,i) >= 0.0 )  THEN
312             phi_m = 1.0 + 5.0 * rif(j,i)
313          ELSE
314             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
315          ENDIF
316       ENDIF
317
318!
319!--    Calculate the mixing length (for dissipation)
320       DO  k = nzb_s_inner(j,i)+1, nzt
321          dvar_dz = atmos_ocean_sign * &
322                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
323          IF ( dvar_dz > 0.0 ) THEN
324             IF ( use_reference )  THEN
325                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
326                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
327             ELSE
328                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
329                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
330             ENDIF
331          ELSE
332             l_stable = l_grid(k)
333          ENDIF
334!
335!--       Adjustment of the mixing length
336          IF ( wall_adjustment )  THEN
337             l(k)  = MIN( wall_adjustment_factor *                     &
338                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
339                          l_stable )
340             ll(k) = MIN( wall_adjustment_factor *                     &
341                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
342          ELSE
343             l(k)  = MIN( l_grid(k), l_stable )
344             ll(k) = l_grid(k)
345          ENDIF
346          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
347             l(k)  = MIN( l(k),  kappa * &
348                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
349             ll(k) = MIN( ll(k), kappa * &
350                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
351          ENDIF
352
353!
354!--       Calculate the tendency term
355          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
356                           SQRT( e(k,j,i) ) / l(k)
357
358          tend(k,j,i) = tend(k,j,i)                                           &
359                                       + (                                    &
360                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
361                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
362                                         ) * ddx2                             &
363                                       + (                                    &
364                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
365                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
366                                         ) * ddy2                             &
367                                       + (                                    &
368              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
369            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
370                                         ) * ddzw(k)                          &
371                                       - dissipation(k)
372
373       ENDDO
374
375!
376!--    Store dissipation if needed for calculating the sgs particle velocities
377       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
378          DO  k = nzb_s_inner(j,i)+1, nzt
379             diss(k,j,i) = dissipation(k)
380          ENDDO
381!
382!--       Boundary condition for dissipation
383          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
384       ENDIF
385
386    END SUBROUTINE diffusion_e_ij
387
388 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.