source: palm/trunk/SOURCE/diffusivities.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: 5.8 KB
Line 
1 SUBROUTINE diffusivities( var, var_reference )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Adjustment of mixing length calculation for the ocean version.
7! This is also a bugfix, because the height above the topography is now
8! used instead of the height above level k=0.
9! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
10! use_pt_reference renamed use_reference
11!
12! Former revisions:
13! -----------------
14! $Id: diffusivities.f90 97 2007-06-21 08:23:15Z raasch $
15!
16! 57 2007-03-09 12:05:41Z raasch
17! Reference temperature pt_reference can be used in buoyancy term
18!
19! RCS Log replace by Id keyword, revision history cleaned up
20!
21! Revision 1.24  2006/04/26 12:16:26  raasch
22! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
23!
24! Revision 1.1  1997/09/19 07:41:10  raasch
25! Initial revision
26!
27!
28! Description:
29! ------------
30! Computation of the turbulent diffusion coefficients for momentum and heat
31! according to Prandtl-Kolmogorov
32!------------------------------------------------------------------------------!
33
34    USE arrays_3d
35    USE control_parameters
36    USE grid_variables
37    USE indices
38    USE pegrid
39    USE statistics
40
41    IMPLICIT NONE
42
43    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
44
45    REAL    ::  dvar_dz, l_stable, var_reference
46
47    REAL, SAVE ::  phi_m = 1.0
48
49    REAL    ::  var(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
50
51    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
52
53
54!
55!-- Default thread number in case of one thread
56    tn = 0
57
58!
59!-- Initialization for calculation of the mixing length profile
60    sums_l_l = 0.0
61
62!
63!-- Compute the turbulent diffusion coefficient for momentum
64    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
65!$  tn = omp_get_thread_num()
66
67    !$OMP DO
68    DO  i = nxl-1, nxr+1
69       DO  j = nys-1, nyn+1
70
71!
72!--       Compute the Phi-function for a possible adaption of the mixing length
73!--       to the Prandtl mixing length
74          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
75             IF ( rif(j,i) >= 0.0 )  THEN
76                phi_m = 1.0 + 5.0 * rif(j,i)
77             ELSE
78                phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
79             ENDIF
80          ENDIF
81         
82!
83!--       Introduce an optional minimum tke
84          IF ( e_min > 0.0 )  THEN
85             DO  k = nzb_s_inner(j,i)+1, nzt
86                e(k,j,i) = MAX( e(k,j,i), e_min )
87             ENDDO
88          ENDIF
89
90!
91!--       Calculate square root of e in a seperate loop, because it is used
92!--       twice in the next loop (better vectorization)
93          DO  k = nzb_s_inner(j,i)+1, nzt
94             sqrt_e(k) = SQRT( e(k,j,i) )
95          ENDDO
96
97!
98!--       Determine the mixing length
99          DO  k = nzb_s_inner(j,i)+1, nzt
100             dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
101                       ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
102             IF ( dvar_dz > 0.0 ) THEN
103                IF ( use_reference )  THEN
104                   l_stable = 0.76 * sqrt_e(k) / &
105                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
106                ELSE
107                   l_stable = 0.76 * sqrt_e(k) / &
108                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
109                ENDIF
110             ELSE
111                l_stable = l_grid(k)
112             ENDIF
113!
114!--          Adjustment of the mixing length
115             IF ( wall_adjustment )  THEN
116                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
117                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
118             ELSE
119                l(k)  = MIN( l_grid(k), l_stable )
120                ll(k) = l_grid(k)
121             ENDIF
122             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
123                l(k)  = MIN( l(k),  kappa * &
124                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
125                ll(k) = MIN( ll(k), kappa * &
126                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
127             ENDIF
128
129!
130!--          Compute diffusion coefficients for momentum and heat
131             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
132             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
133
134          ENDDO
135
136!
137!--       Summation for averaged profile (cf. flow_statistics)
138          DO  sr = 0, statistic_regions
139             IF ( rmask(j,i,sr) /= 0.0 )  THEN
140                DO  k = nzb_s_outer(j,i)+1, nzt
141                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
142                ENDDO
143             ENDIF
144          ENDDO
145
146       ENDDO
147    ENDDO
148
149    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
150                                                ! data output
151
152    !$OMP END PARALLEL
153
154!
155!-- Set vertical boundary values (Neumann conditions both at bottom and top).
156!-- Horizontal boundary conditions at vertical walls are not set because
157!-- so far vertical walls require usage of a Prandtl-layer where the boundary
158!-- values of the diffusivities are not needed
159    !$OMP PARALLEL DO
160    DO  i = nxl-1, nxr+1
161       DO  j = nys-1, nyn+1
162          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
163          km(nzt+1,j,i)            = km(nzt,j,i)
164          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
165          kh(nzt+1,j,i)            = kh(nzt,j,i)
166       ENDDO
167    ENDDO
168
169!
170!-- Set Neumann boundary conditions at the outflow boundaries in case of
171!-- non-cyclic lateral boundaries
172    IF ( outflow_l )  THEN
173       km(:,:,nxl-1) = km(:,:,nxl)
174       kh(:,:,nxl-1) = kh(:,:,nxl)
175    ENDIF
176    IF ( outflow_r )  THEN
177       km(:,:,nxr+1) = km(:,:,nxr)
178       kh(:,:,nxr+1) = kh(:,:,nxr)
179    ENDIF
180    IF ( outflow_s )  THEN
181       km(:,nys-1,:) = km(:,nys,:)
182       kh(:,nys-1,:) = kh(:,nys,:)
183    ENDIF
184    IF ( outflow_n )  THEN
185       km(:,nyn+1,:) = km(:,nyn,:)
186       kh(:,nyn+1,:) = kh(:,nyn,:)
187    ENDIF
188
189
190 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.