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

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

preliminary uncomplete changes for ocean version

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