Changeset 1257 for palm/trunk/SOURCE/timestep.f90
- Timestamp:
- Nov 8, 2013 3:18:40 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/timestep.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/timestep.f90
r1093 r1257 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! openacc porting 23 ! bugfix for calculation of advective timestep in case of vertically stretched 24 ! grids 23 25 ! 24 26 ! Former revisions: … … 93 95 IMPLICIT NONE 94 96 95 INTEGER :: i, j, k, u_max_cfl_ijk(3), v_max_cfl_ijk(3) 96 97 REAL :: div, dt_diff, dt_diff_l, dt_plant_canopy, dt_plant_canopy_l, & 98 dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w, & 99 dt_u, dt_v, dt_w, u_max_cfl, value, v_max_cfl 97 INTEGER :: i, j, k 98 99 REAL :: div, dt_diff, dt_diff_l, dt_plant_canopy, dt_plant_canopy_l, & 100 dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w, & 101 dt_u, dt_u_l, dt_v, dt_v_l, dt_w, dt_w_l, u_gtrans_l, u_max_l, & 102 u_min_l, value, v_gtrans_l, v_max_l, v_min_l, w_max_l, w_min_l 100 103 101 104 REAL, DIMENSION(2) :: uv_gtrans, uv_gtrans_l 105 REAL, DIMENSION(3) :: reduce, reduce_l 102 106 REAL, DIMENSION(nzb+1:nzt) :: dxyz2_min 103 107 … … 127 131 ! 128 132 !-- Averaging over the entire model domain. 129 uv_gtrans_l = 0.0 133 u_gtrans_l = 0.0 134 v_gtrans_l = 0.0 135 !$acc parallel present( u, v ) 130 136 DO i = nxl, nxr 131 137 DO j = nys, nyn 132 138 DO k = nzb+1, nzt 133 u v_gtrans_l(1) = uv_gtrans_l(1)+ u(k,j,i)134 uv_gtrans_l(2) = uv_gtrans_l(2)+ v(k,j,i)139 u_gtrans_l = u_gtrans_l + u(k,j,i) 140 v_gtrans_l = v_gtrans_l + v(k,j,i) 135 141 ENDDO 136 142 ENDDO 137 143 ENDDO 138 uv_gtrans_l = uv_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) ) 144 !$acc end parallel 145 uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) ) 146 uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) ) 139 147 #if defined( __parallel ) 140 148 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) … … 151 159 152 160 ! 153 !-- Determine the maxima of the velocity components. 161 !-- Determine the maxima of the velocity components, including their 162 !-- grid index positions. 163 #if defined( __openacc ) 164 IF ( dt_fixed ) THEN ! otherwise do it further below for better cache usage 165 u_max_l = -999999.9 166 u_min_l = 999999.9 167 v_max_l = -999999.9 168 v_min_l = 999999.9 169 w_max_l = -999999.9 170 w_min_l = 999999.9 171 !$acc parallel present( u, v, w ) 172 DO i = nxl, nxr 173 DO j = nys, nyn 174 DO k = nzb+1, nzt 175 u_max_l = MAX( u_max_l, u(k,j,i) ) 176 u_min_l = MIN( u_min_l, u(k,j,i) ) 177 v_max_l = MAX( v_max_l, v(k,j,i) ) 178 v_min_l = MIN( v_min_l, v(k,j,i) ) 179 w_max_l = MAX( w_max_l, w(k,j,i) ) 180 w_min_l = MIN( w_min_l, w(k,j,i) ) 181 ENDDO 182 ENDDO 183 ENDDO 184 !$acc end parallel 185 #if defined( __parallel ) 186 reduce_l(1) = u_max_l 187 reduce_l(2) = v_max_l 188 reduce_l(3) = w_max_l 189 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 190 CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr ) 191 u_max = reduce(1) 192 v_max = reduce(2) 193 w_max = reduce(3) 194 reduce_l(1) = u_min_l 195 reduce_l(2) = v_min_l 196 reduce_l(3) = w_min_l 197 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 198 CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr ) 199 IF ( ABS( reduce(1) ) > u_max ) u_max = reduce(1) 200 IF ( ABS( reduce(2) ) > v_max ) v_max = reduce(2) 201 IF ( ABS( reduce(3) ) > w_max ) w_max = reduce(3) 202 #else 203 IF ( ABS( u_min_l ) > u_max_l ) THEN 204 u_max = u_min_l 205 ELSE 206 u_max = u_max_l 207 ENDIF 208 IF ( ABS( v_min_l ) > v_max_l ) THEN 209 v_max = v_min_l 210 ELSE 211 v_max = v_max_l 212 ENDIF 213 IF ( ABS( w_min_l ) > w_max_l ) THEN 214 w_max = w_min_l 215 ELSE 216 w_max = w_max_l 217 ENDIF 218 #endif 219 ENDIF 220 #else 154 221 CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0, & 155 222 u_max, u_max_ijk ) … … 158 225 CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0, & 159 226 w_max, w_max_ijk ) 160 161 ! 162 !-- In case of Galilei transformation, the horizontal velocity maxima have 163 !-- to be calculated from the transformed horizontal velocities 164 IF ( galilei_transformation ) THEN 165 CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'absoff', & 166 u_gtrans, u_max_cfl, u_max_cfl_ijk ) 167 CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'absoff', & 168 v_gtrans, v_max_cfl, v_max_cfl_ijk ) 169 ELSE 170 u_max_cfl = u_max 171 v_max_cfl = v_max 172 u_max_cfl_ijk = u_max_ijk 173 v_max_cfl_ijk = v_max_ijk 174 ENDIF 175 227 #endif 176 228 177 229 IF ( .NOT. dt_fixed ) THEN 230 #if defined( __openacc ) 178 231 ! 179 232 !-- Variable time step: 180 ! 181 !-- For each component, compute the maximum time step according to the 182 !-- CFL-criterion. 183 dt_u = dx / ( ABS( u_max_cfl ) + 1.0E-10 ) 184 dt_v = dy / ( ABS( v_max_cfl ) + 1.0E-10 ) 185 dt_w = dzu(MAX( 1, w_max_ijk(1) )) / ( ABS( w_max ) + 1.0E-10 ) 233 !-- Calculate the maximum time step according to the CFL-criterion, 234 !-- individually for each velocity component 235 dt_u_l = 999999.9 236 dt_v_l = 999999.9 237 dt_w_l = 999999.9 238 u_max_l = -999999.9 239 u_min_l = 999999.9 240 v_max_l = -999999.9 241 v_min_l = 999999.9 242 w_max_l = -999999.9 243 w_min_l = 999999.9 244 !$acc parallel loop collapse(3) present( u, v, w ) 245 DO i = nxl, nxr 246 DO j = nys, nyn 247 DO k = nzb+1, nzt 248 dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) ) 249 dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) ) 250 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10 ) ) ) 251 u_max_l = MAX( u_max_l, u(k,j,i) ) 252 u_min_l = MIN( u_min_l, u(k,j,i) ) 253 v_max_l = MAX( v_max_l, v(k,j,i) ) 254 v_min_l = MIN( v_min_l, v(k,j,i) ) 255 w_max_l = MAX( w_max_l, w(k,j,i) ) 256 w_min_l = MIN( w_min_l, w(k,j,i) ) 257 ENDDO 258 ENDDO 259 ENDDO 260 !$acc end parallel 261 262 #if defined( __parallel ) 263 reduce_l(1) = dt_u_l 264 reduce_l(2) = dt_v_l 265 reduce_l(3) = dt_w_l 266 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 267 CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr ) 268 dt_u = reduce(1) 269 dt_v = reduce(2) 270 dt_w = reduce(3) 271 272 reduce_l(1) = u_max_l 273 reduce_l(2) = v_max_l 274 reduce_l(3) = w_max_l 275 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 276 CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr ) 277 u_max = reduce(1) 278 v_max = reduce(2) 279 w_max = reduce(3) 280 reduce_l(1) = u_min_l 281 reduce_l(2) = v_min_l 282 reduce_l(3) = w_min_l 283 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 284 CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr ) 285 IF ( ABS( reduce(1) ) > u_max ) u_max = reduce(1) 286 IF ( ABS( reduce(2) ) > v_max ) v_max = reduce(2) 287 IF ( ABS( reduce(3) ) > w_max ) w_max = reduce(3) 288 #else 289 dt_u = dt_u_l 290 dt_v = dt_v_l 291 dt_w = dt_w_l 292 293 IF ( ABS( u_min_l ) > u_max_l ) THEN 294 u_max = u_min_l 295 ELSE 296 u_max = u_max_l 297 ENDIF 298 IF ( ABS( v_min_l ) > v_max_l ) THEN 299 v_max = v_min_l 300 ELSE 301 v_max = v_max_l 302 ENDIF 303 IF ( ABS( w_min_l ) > w_max_l ) THEN 304 w_max = w_min_l 305 ELSE 306 w_max = w_max_l 307 ENDIF 308 #endif 309 310 #else 311 ! 312 !-- Variable time step: 313 !-- Calculate the maximum time step according to the CFL-criterion, 314 !-- individually for each velocity component 315 dt_u_l = 999999.9 316 dt_v_l = 999999.9 317 dt_w_l = 999999.9 318 DO i = nxl, nxr 319 DO j = nys, nyn 320 DO k = nzb+1, nzt 321 dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) ) 322 dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) ) 323 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10 ) ) ) 324 ENDDO 325 ENDDO 326 ENDDO 327 328 #if defined( __parallel ) 329 reduce_l(1) = dt_u_l 330 reduce_l(2) = dt_v_l 331 reduce_l(3) = dt_w_l 332 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 333 CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr ) 334 dt_u = reduce(1) 335 dt_v = reduce(2) 336 dt_w = reduce(3) 337 #else 338 dt_u = dt_u_l 339 dt_v = dt_v_l 340 dt_w = dt_w_l 341 #endif 342 343 #endif 186 344 187 345 ! … … 199 357 !$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l) 200 358 !$OMP DO 359 !$acc parallel loop collapse(3) present( kh, km ) 201 360 DO i = nxl, nxr 202 361 DO j = nys, nyn 203 362 DO k = nzb+1, nzt 204 value = dxyz2_min(k) / ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 ) 205 206 dt_diff_l = MIN( value, dt_diff_l ) 363 dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / & 364 ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 ) ) 207 365 ENDDO 208 366 ENDDO 209 367 ENDDO 368 !$acc end parallel 210 369 !$OMP END PARALLEL 211 370 #if defined( __parallel ) … … 334 493 '&dt_diff = ', dt_diff, ' s', & 335 494 '&dt_plant_canopy = ', dt_plant_canopy, ' s', & 336 '&u_max _cfl = ', u_max_cfl, ' m/s k=', u_max_cfl_ijk(1),&495 '&u_max = ', u_max, ' m/s k=', u_max_ijk(1), & 337 496 ' j=', u_max_ijk(2), ' i=', u_max_ijk(3), & 338 '&v_max _cfl = ', v_max_cfl, ' m/s k=', v_max_cfl_ijk(1),&497 '&v_max = ', v_max, ' m/s k=', v_max_ijk(1), & 339 498 ' j=', v_max_ijk(2), ' i=', v_max_ijk(3), & 340 '&w_max = ', w_max, ' m/s k=', w_max_ijk(1),&499 '&w_max = ', w_max, ' m/s k=', w_max_ijk(1), & 341 500 ' j=', w_max_ijk(2), ' i=', w_max_ijk(3) 342 501 CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
Note: See TracChangeset
for help on using the changeset viewer.
