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

Last change on this file since 550 was 484, checked in by raasch, 14 years ago

typo in file headers removed

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