source: palm/trunk/SOURCE/diffusivities.f90 @ 1017

Last change on this file since 1017 was 1017, checked in by raasch, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 5.9 KB
Line 
1 SUBROUTINE diffusivities( var, var_reference )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusivities.f90 1017 2012-09-27 11:28:50Z raasch $
11!
12! 1015 2012-09-27 09:23:24Z raasch
13! OpenACC statements added + code changes required for GPU optimization,
14! adjustment of mixing length to the Prandtl mixing length at first grid point
15! above ground removed
16!
17! 667 2010-12-23 12:06:00Z suehring/gryschka
18! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
19!
20! 137 2007-11-28 08:50:10Z letzel
21! Bugfix for summation of sums_l_l for flow_statistics
22! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
23!
24! 97 2007-06-21 08:23:15Z raasch
25! Adjustment of mixing length calculation for the ocean version.
26! This is also a bugfix, because the height above the topography is now
27! used instead of the height above level k=0.
28! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
29! use_pt_reference renamed use_reference
30!
31! 57 2007-03-09 12:05:41Z raasch
32! Reference temperature pt_reference can be used in buoyancy term
33!
34! RCS Log replace by Id keyword, revision history cleaned up
35!
36! Revision 1.24  2006/04/26 12:16:26  raasch
37! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
38!
39! Revision 1.1  1997/09/19 07:41:10  raasch
40! Initial revision
41!
42!
43! Description:
44! ------------
45! Computation of the turbulent diffusion coefficients for momentum and heat
46! according to Prandtl-Kolmogorov
47!------------------------------------------------------------------------------!
48
49    USE arrays_3d
50    USE control_parameters
51    USE grid_variables
52    USE indices
53    USE pegrid
54    USE statistics
55
56    IMPLICIT NONE
57
58    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
59
60    REAL    ::  dvar_dz, l, ll, l_stable, sqrt_e, var_reference
61
62    REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
63
64
65!
66!-- Default thread number in case of one thread
67    tn = 0
68
69!
70!-- Initialization for calculation of the mixing length profile
71    sums_l_l = 0.0
72
73!
74!-- Compute the turbulent diffusion coefficient for momentum
75    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn)
76!$  tn = omp_get_thread_num()
77
78!
79!-- Data declerations for accelerators
80    !$acc data present( dd2zu, e, km, kh, l_grid, l_wall, nzb_s_inner, rif, var )
81    !$acc kernels
82
83!
84!-- Introduce an optional minimum tke
85    IF ( e_min > 0.0 )  THEN
86       !$OMP DO
87       !$acc loop
88       DO  i = nxlg, nxrg
89          DO  j = nysg, nyng
90             !$acc loop vector( 32 )
91             DO  k = 1, nzt
92                IF ( k > nzb_s_inner(j,i) )  THEN
93                   e(k,j,i) = MAX( e(k,j,i), e_min )
94                ENDIF
95             ENDDO
96          ENDDO
97       ENDDO
98    ENDIF
99
100    !$OMP DO
101    !$acc loop
102    DO  i = nxlg, nxrg
103       DO  j = nysg, nyng
104          !$acc loop vector( 32 )
105          DO  k = 1, nzt
106
107             IF ( k > nzb_s_inner(j,i) )  THEN
108
109                sqrt_e = SQRT( e(k,j,i) )
110!
111!--             Determine the mixing length
112                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
113                          ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
114                IF ( dvar_dz > 0.0 ) THEN
115                   IF ( use_reference )  THEN
116                      l_stable = 0.76 * sqrt_e / &
117                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
118                   ELSE
119                      l_stable = 0.76 * sqrt_e / &
120                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
121                   ENDIF
122                ELSE
123                   l_stable = l_grid(k)
124                ENDIF
125!
126!--             Adjustment of the mixing length
127                IF ( wall_adjustment )  THEN
128                   = MIN( l_wall(k,j,i), l_grid(k), l_stable )
129                   ll = MIN( l_wall(k,j,i), l_grid(k) )
130                ELSE
131                   = MIN( l_grid(k), l_stable )
132                   ll = l_grid(k)
133                ENDIF
134
135      !
136      !--       Compute diffusion coefficients for momentum and heat
137                km(k,j,i) = 0.1 * l * sqrt_e
138                kh(k,j,i) = ( 1.0 + 2.0 * l / ll ) * km(k,j,i)
139
140             ENDIF
141
142#if ! defined( __openacc )
143!
144!++          Statistics still have to be realized for accelerators
145!--          Summation for averaged profile (cf. flow_statistics)
146             DO  sr = 0, statistic_regions
147                sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)
148             ENDDO
149#endif
150
151          ENDDO
152       ENDDO
153    ENDDO
154
155#if ! defined( __openacc )
156!
157!++ Statistics still have to be realized for accelerators
158    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
159                                                  ! data output
160#endif
161    !$OMP END PARALLEL
162
163!
164!-- Set vertical boundary values (Neumann conditions both at bottom and top).
165!-- Horizontal boundary conditions at vertical walls are not set because
166!-- so far vertical walls require usage of a Prandtl-layer where the boundary
167!-- values of the diffusivities are not needed
168    !$OMP PARALLEL DO
169    !$acc loop
170    DO  i = nxlg, nxrg
171       DO  j = nysg, nyng
172          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
173          km(nzt+1,j,i)            = km(nzt,j,i)
174          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
175          kh(nzt+1,j,i)            = kh(nzt,j,i)
176       ENDDO
177    ENDDO
178
179!
180!-- Set Neumann boundary conditions at the outflow boundaries in case of
181!-- non-cyclic lateral boundaries
182    IF ( outflow_l )  THEN
183       km(:,:,nxl-1) = km(:,:,nxl)
184       kh(:,:,nxl-1) = kh(:,:,nxl)
185    ENDIF
186    IF ( outflow_r )  THEN
187       km(:,:,nxr+1) = km(:,:,nxr)
188       kh(:,:,nxr+1) = kh(:,:,nxr)
189    ENDIF
190    IF ( outflow_s )  THEN
191       km(:,nys-1,:) = km(:,nys,:)
192       kh(:,nys-1,:) = kh(:,nys,:)
193    ENDIF
194    IF ( outflow_n )  THEN
195       km(:,nyn+1,:) = km(:,nyn,:)
196       kh(:,nyn+1,:) = kh(:,nyn,:)
197    ENDIF
198
199    !$acc end kernels
200    !$acc end data
201
202 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.