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

Last change on this file since 821 was 791, 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 791 2011-11-29 03:33:42Z maronga $
11!
12! 790 2011-11-29 03:11:20Z raasch
13! diss is also calculated in case that the Wang kernel is used
14!
15! 667 2010-12-23 12:06:00Z suehring/gryschka
16! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
17!
18! 97 2007-06-21 08:23:15Z raasch
19! Adjustment of mixing length calculation for the ocean version. zw added to
20! argument list.
21! This is also a bugfix, because the height above the topography is now
22! used instead of the height above level k=0.
23! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
24! use_pt_reference renamed use_reference
25!
26! 65 2007-03-13 12:11:43Z raasch
27! Reference temperature pt_reference can be used in buoyancy term
28!
29! 20 2007-02-26 00:12:32Z raasch
30! Bugfix: ddzw dimensioned 1:nzt"+1"
31! Calculation extended for gridpoint nzt
32!
33! RCS Log replace by Id keyword, revision history cleaned up
34!
35! Revision 1.18  2006/08/04 14:29:43  raasch
36! dissipation is stored in extra array diss if needed later on for calculating
37! the sgs particle velocities
38!
39! Revision 1.1  1997/09/19 07:40:24  raasch
40! Initial revision
41!
42!
43! Description:
44! ------------
45! Diffusion- and dissipation terms for the TKE
46!------------------------------------------------------------------------------!
47
48    PRIVATE
49    PUBLIC diffusion_e
50   
51
52    INTERFACE diffusion_e
53       MODULE PROCEDURE diffusion_e
54       MODULE PROCEDURE diffusion_e_ij
55    END INTERFACE diffusion_e
56 
57 CONTAINS
58
59
60!------------------------------------------------------------------------------!
61! Call for all grid points
62!------------------------------------------------------------------------------!
63    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, var, &
64                            var_reference, rif, tend, zu, zw )
65
66       USE control_parameters
67       USE grid_variables
68       USE indices
69       USE particle_attributes
70
71       IMPLICIT NONE
72
73       INTEGER ::  i, j, k
74       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
75       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
76                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
77       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
78       REAL, DIMENSION(:,:), POINTER   ::  rif
79       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
80       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
81 
82
83!
84!--    This if clause must be outside the k-loop because otherwise
85!--    runtime errors occur with -C hopt on NEC
86       IF ( use_reference )  THEN
87
88          DO  i = nxl, nxr
89             DO  j = nys, nyn
90!
91!--             First, calculate phi-function for eventually adjusting the &
92!--             mixing length to the prandtl mixing length
93                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
94                   IF ( rif(j,i) >= 0.0 )  THEN
95                      phi_m = 1.0 + 5.0 * rif(j,i)
96                   ELSE
97                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
98                   ENDIF
99                ENDIF
100
101                DO  k = nzb_s_inner(j,i)+1, nzt
102!
103!--                Calculate the mixing length (for dissipation)
104                   dvar_dz = atmos_ocean_sign * &
105                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
106                   IF ( dvar_dz > 0.0 ) THEN
107                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
108                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
109                   ELSE
110                      l_stable = l_grid(k)
111                   ENDIF
112!
113!--                Adjustment of the mixing length
114                   IF ( wall_adjustment )  THEN
115                      l(k,j)  = MIN( wall_adjustment_factor *          &
116                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
117                                     l_grid(k), l_stable )
118                      ll(k,j) = MIN( wall_adjustment_factor *          &
119                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
120                                     l_grid(k) )
121                   ELSE
122                      l(k,j)  = MIN( l_grid(k), l_stable )
123                      ll(k,j) = l_grid(k)
124                   ENDIF
125                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
126                      l(k,j)  = MIN( l(k,j),  kappa *                          &
127                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
128                                              / phi_m )
129                      ll(k,j) = MIN( ll(k,j), kappa *                          &
130                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
131                                              / phi_m )
132                   ENDIF
133
134                ENDDO
135             ENDDO
136
137!
138!--          Calculate the tendency terms
139             DO  j = nys, nyn
140                DO  k = nzb_s_inner(j,i)+1, nzt
141
142                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
143                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
144
145                    tend(k,j,i) = tend(k,j,i)                                  &
146                                        + (                                    &
147                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
148                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
149                                          ) * ddx2                             &
150                                        + (                                    &
151                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
152                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
153                                          ) * ddy2                             &
154                                        + (                                    &
155               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
156             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
157                                          ) * ddzw(k)                          &
158                             - dissipation(k,j)
159
160                ENDDO
161             ENDDO
162
163!
164!--          Store dissipation if needed for calculating the sgs particle
165!--          velocities
166             IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
167                DO  j = nys, nyn
168                   DO  k = nzb_s_inner(j,i)+1, nzt
169                      diss(k,j,i) = dissipation(k,j)
170                   ENDDO
171                ENDDO
172             ENDIF
173
174          ENDDO
175
176       ELSE
177
178          DO  i = nxl, nxr
179             DO  j = nys, nyn
180!
181!--             First, calculate phi-function for eventually adjusting the &
182!--             mixing length to the prandtl mixing length
183                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
184                   IF ( rif(j,i) >= 0.0 )  THEN
185                      phi_m = 1.0 + 5.0 * rif(j,i)
186                   ELSE
187                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
188                   ENDIF
189                ENDIF
190
191                DO  k = nzb_s_inner(j,i)+1, nzt
192!
193!--                Calculate the mixing length (for dissipation)
194                   dvar_dz = atmos_ocean_sign * &
195                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
196                   IF ( dvar_dz > 0.0 ) THEN
197                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
198                                        SQRT( g / var(k,j,i) * dvar_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,j)  = MIN( wall_adjustment_factor *          &
206                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
207                                     l_grid(k), l_stable )
208                      ll(k,j) = MIN( wall_adjustment_factor *          &
209                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
210                                     l_grid(k) )
211                   ELSE
212                      l(k,j)  = MIN( l_grid(k), l_stable )
213                      ll(k,j) = l_grid(k)
214                   ENDIF
215                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
216                      l(k,j)  = MIN( l(k,j),  kappa *                          &
217                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
218                                              / phi_m )
219                      ll(k,j) = MIN( ll(k,j), kappa *                          &
220                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
221                                              / phi_m )
222                   ENDIF
223
224                ENDDO
225             ENDDO
226
227!
228!--          Calculate the tendency terms
229             DO  j = nys, nyn
230                DO  k = nzb_s_inner(j,i)+1, nzt
231
232                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
233                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
234
235                    tend(k,j,i) = tend(k,j,i)                                  &
236                                        + (                                    &
237                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
238                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
239                                          ) * ddx2                             &
240                                        + (                                    &
241                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
242                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
243                                          ) * ddy2                             &
244                                        + (                                    &
245               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
246             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
247                                          ) * ddzw(k)                          &
248                             - dissipation(k,j)
249
250                ENDDO
251             ENDDO
252
253!
254!--          Store dissipation if needed for calculating the sgs particle
255!--          velocities
256             IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
257                DO  j = nys, nyn
258                   DO  k = nzb_s_inner(j,i)+1, nzt
259                      diss(k,j,i) = dissipation(k,j)
260                   ENDDO
261                ENDDO
262             ENDIF
263
264          ENDDO
265
266       ENDIF
267
268!
269!--    Boundary condition for dissipation
270       IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
271          DO  i = nxl, nxr
272             DO  j = nys, nyn
273                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
274             ENDDO
275          ENDDO
276       ENDIF
277
278    END SUBROUTINE diffusion_e
279
280
281!------------------------------------------------------------------------------!
282! Call for grid point i,j
283!------------------------------------------------------------------------------!
284    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
285                               var, var_reference, rif, tend, zu, zw )
286
287       USE control_parameters
288       USE grid_variables
289       USE indices
290       USE particle_attributes
291
292       IMPLICIT NONE
293
294       INTEGER         ::  i, j, k
295       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
296       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
297                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
298       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
299       REAL, DIMENSION(:,:), POINTER   ::  rif
300       REAL, DIMENSION(:,:,:), POINTER ::  e, km, 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_collision_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.