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

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

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:


inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:


Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

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