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

Last change on this file since 750 was 668, checked in by suehring, 14 years ago

last commit documented

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