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

Last change on this file since 1002 was 1002, checked in by raasch, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 14.8 KB
Line 
1 MODULE diffusion_e_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_e.f90 1002 2012-09-13 15:12:24Z raasch $
11!
12! 1001 2012-09-13 14:08:46Z raasch
13! most arrays comunicated by module instead of parameter list
14!
15! 825 2012-02-19 03:03:44Z raasch
16! wang_collision_kernel renamed wang_kernel
17!
18! 790 2011-11-29 03:11:20Z raasch
19! diss is also calculated in case that the Wang kernel is used
20!
21! 667 2010-12-23 12:06:00Z suehring/gryschka
22! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
23!
24! 97 2007-06-21 08:23:15Z raasch
25! Adjustment of mixing length calculation for the ocean version. zw added to
26! argument list.
27! This is also a bugfix, because the height above the topography is now
28! used instead of the height above level k=0.
29! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
30! use_pt_reference renamed use_reference
31!
32! 65 2007-03-13 12:11:43Z raasch
33! Reference temperature pt_reference can be used in buoyancy term
34!
35! 20 2007-02-26 00:12:32Z raasch
36! Bugfix: ddzw dimensioned 1:nzt"+1"
37! Calculation extended for gridpoint nzt
38!
39! RCS Log replace by Id keyword, revision history cleaned up
40!
41! Revision 1.18  2006/08/04 14:29:43  raasch
42! dissipation is stored in extra array diss if needed later on for calculating
43! the sgs particle velocities
44!
45! Revision 1.1  1997/09/19 07:40:24  raasch
46! Initial revision
47!
48!
49! Description:
50! ------------
51! Diffusion- and dissipation terms for the TKE
52!------------------------------------------------------------------------------!
53
54    PRIVATE
55    PUBLIC diffusion_e
56   
57
58    INTERFACE diffusion_e
59       MODULE PROCEDURE diffusion_e
60       MODULE PROCEDURE diffusion_e_ij
61    END INTERFACE diffusion_e
62 
63 CONTAINS
64
65
66!------------------------------------------------------------------------------!
67! Call for all grid points
68!------------------------------------------------------------------------------!
69    SUBROUTINE diffusion_e( var, var_reference )
70
71       USE arrays_3d
72       USE control_parameters
73       USE grid_variables
74       USE indices
75       USE particle_attributes
76
77       IMPLICIT NONE
78
79       INTEGER ::  i, j, k
80       REAL    ::  dvar_dz, l_stable, phi_m, var_reference
81
82       REAL, DIMENSION(:,:,:), POINTER ::  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, var, var_reference )
288
289       USE arrays_3d
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
300       REAL, DIMENSION(:,:,:), POINTER ::  var
301       REAL, DIMENSION(nzb+1:nzt) ::  dissipation, l, ll
302
303
304!
305!--    First, calculate phi-function for eventually adjusting the mixing length
306!--    to the prandtl mixing length
307       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
308          IF ( rif(j,i) >= 0.0 )  THEN
309             phi_m = 1.0 + 5.0 * rif(j,i)
310          ELSE
311             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
312          ENDIF
313       ENDIF
314
315!
316!--    Calculate the mixing length (for dissipation)
317       DO  k = nzb_s_inner(j,i)+1, nzt
318          dvar_dz = atmos_ocean_sign * &
319                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
320          IF ( dvar_dz > 0.0 ) THEN
321             IF ( use_reference )  THEN
322                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
323                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
324             ELSE
325                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
326                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
327             ENDIF
328          ELSE
329             l_stable = l_grid(k)
330          ENDIF
331!
332!--       Adjustment of the mixing length
333          IF ( wall_adjustment )  THEN
334             l(k)  = MIN( wall_adjustment_factor *                     &
335                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
336                          l_stable )
337             ll(k) = MIN( wall_adjustment_factor *                     &
338                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
339          ELSE
340             l(k)  = MIN( l_grid(k), l_stable )
341             ll(k) = l_grid(k)
342          ENDIF
343          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
344             l(k)  = MIN( l(k),  kappa * &
345                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
346             ll(k) = MIN( ll(k), kappa * &
347                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
348          ENDIF
349
350!
351!--       Calculate the tendency term
352          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
353                           SQRT( e(k,j,i) ) / l(k)
354
355          tend(k,j,i) = tend(k,j,i)                                           &
356                                       + (                                    &
357                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
358                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
359                                         ) * ddx2                             &
360                                       + (                                    &
361                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
362                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
363                                         ) * ddy2                             &
364                                       + (                                    &
365              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
366            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
367                                         ) * ddzw(k)                          &
368                                       - dissipation(k)
369
370       ENDDO
371
372!
373!--    Store dissipation if needed for calculating the sgs particle velocities
374       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
375          DO  k = nzb_s_inner(j,i)+1, nzt
376             diss(k,j,i) = dissipation(k)
377          ENDDO
378!
379!--       Boundary condition for dissipation
380          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
381       ENDIF
382
383    END SUBROUTINE diffusion_e_ij
384
385 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.