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

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

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

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