Changeset 4845 for palm/trunk/SOURCE/dynamics_mod.f90
- Timestamp:
- Jan 18, 2021 11:15:37 AM (3 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/dynamics_mod.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/dynamics_mod.f90
r4843 r4845 24 24 ! ----------------- 25 25 ! $Id$ 26 ! radiation boundary condition always uses maximum phase velocity, respective code for calculating 27 ! the phase velocity is removed 28 ! 29 ! 4843 2021-01-15 15:22:11Z raasch 26 30 ! local namelist parameter added to switch off the module although the respective module namelist 27 31 ! appears in the namelist file … … 83 87 84 88 USE arrays_3d, & 85 ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, & 86 diss, & 89 ONLY: diss, & 87 90 diss_p, & 88 91 dzu, & … … 94 97 q, q_1, q_2, q_p, & 95 98 s, s_1, s_2, s_p, & 96 u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s,&97 v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s,&98 w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s,&99 u, u_1, u_2, u_init, u_p, & 100 v, v_1, v_2, v_p, v_init, & 101 w, w_1, w_2, w_p, & 99 102 zu 100 103 … … 138 141 rans_mode, & 139 142 rans_tke_e, & 140 tsc, & 141 use_cmax 143 tsc 142 144 143 145 USE exchange_horiz_mod, & … … 896 898 INTEGER(iwp) :: m !< running index surface elements 897 899 898 REAL(wp) :: c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition899 REAL(wp) :: denom !< horizontal gradient of velocity component normal to the outflow boundary900 901 900 ! 902 901 !-- Bottom boundary … … 1128 1127 ! 1129 1128 !-- Radiation boundary conditions for the velocities at the respective outflow. 1130 !-- The phase velocity is either assumed to the maximum phase velocity that ensures numerical 1131 !-- stability (CFL-condition) or calculated after Orlanski(1976) and averaged along the outflow 1132 !-- boundary. 1129 !-- The phase velocity is set to the maximum phase velocity that ensures numerical 1130 !-- stability (CFL-condition), i.e. a Courant number of one is assumed. 1133 1131 IF ( bc_radiation_s ) THEN 1134 1135 IF ( use_cmax ) THEN 1136 u_p(:,-1,:) = u(:,0,:) 1137 v_p(:,0,:) = v(:,1,:) 1138 w_p(:,-1,:) = w(:,0,:) 1139 ELSEIF ( .NOT. use_cmax ) THEN 1140 1141 c_max = dy / dt_3d 1142 1143 c_u_m_l = 0.0_wp 1144 c_v_m_l = 0.0_wp 1145 c_w_m_l = 0.0_wp 1146 1147 c_u_m = 0.0_wp 1148 c_v_m = 0.0_wp 1149 c_w_m = 0.0_wp 1150 1151 ! 1152 !-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow 1153 !-- boundary. 1154 DO k = nzb+1, nzt+1 1155 DO i = nxl, nxr 1156 1157 denom = u_m_s(k,0,i) - u_m_s(k,1,i) 1158 1159 IF ( denom /= 0.0_wp ) THEN 1160 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) ) 1161 IF ( c_u(k,i) < 0.0_wp ) THEN 1162 c_u(k,i) = 0.0_wp 1163 ELSEIF ( c_u(k,i) > c_max ) THEN 1164 c_u(k,i) = c_max 1165 ENDIF 1166 ELSE 1167 c_u(k,i) = c_max 1168 ENDIF 1169 1170 denom = v_m_s(k,1,i) - v_m_s(k,2,i) 1171 1172 IF ( denom /= 0.0_wp ) THEN 1173 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) ) 1174 IF ( c_v(k,i) < 0.0_wp ) THEN 1175 c_v(k,i) = 0.0_wp 1176 ELSEIF ( c_v(k,i) > c_max ) THEN 1177 c_v(k,i) = c_max 1178 ENDIF 1179 ELSE 1180 c_v(k,i) = c_max 1181 ENDIF 1182 1183 denom = w_m_s(k,0,i) - w_m_s(k,1,i) 1184 1185 IF ( denom /= 0.0_wp ) THEN 1186 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) ) 1187 IF ( c_w(k,i) < 0.0_wp ) THEN 1188 c_w(k,i) = 0.0_wp 1189 ELSEIF ( c_w(k,i) > c_max ) THEN 1190 c_w(k,i) = c_max 1191 ENDIF 1192 ELSE 1193 c_w(k,i) = c_max 1194 ENDIF 1195 1196 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 1197 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 1198 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 1199 1200 ENDDO 1201 ENDDO 1202 1203 #if defined( __parallel ) 1204 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1205 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, & 1206 ierr ) 1207 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1208 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, & 1209 ierr ) 1210 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1211 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, & 1212 ierr ) 1213 #else 1214 c_u_m = c_u_m_l 1215 c_v_m = c_v_m_l 1216 c_w_m = c_w_m_l 1217 #endif 1218 1219 c_u_m = c_u_m / (nx+1) 1220 c_v_m = c_v_m / (nx+1) 1221 c_w_m = c_w_m / (nx+1) 1222 1223 ! 1224 !-- Save old timelevels for the next timestep 1225 IF ( intermediate_timestep_count == 1 ) THEN 1226 u_m_s(:,:,:) = u(:,0:1,:) 1227 v_m_s(:,:,:) = v(:,1:2,:) 1228 w_m_s(:,:,:) = w(:,0:1,:) 1229 ENDIF 1230 1231 ! 1232 !-- Calculate the new velocities 1233 DO k = nzb+1, nzt+1 1234 DO i = nxlg, nxrg 1235 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * & 1236 ( u(k,-1,i) - u(k,0,i) ) * ddy 1237 1238 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * & 1239 ( v(k,0,i) - v(k,1,i) ) * ddy 1240 1241 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * & 1242 ( w(k,-1,i) - w(k,0,i) ) * ddy 1243 ENDDO 1244 ENDDO 1245 1246 ! 1247 !-- Bottom boundary at the outflow 1248 IF ( ibc_uv_b == 0 ) THEN 1249 u_p(nzb,-1,:) = 0.0_wp 1250 v_p(nzb,0,:) = 0.0_wp 1251 ELSE 1252 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 1253 v_p(nzb,0,:) = v_p(nzb+1,0,:) 1254 ENDIF 1255 w_p(nzb,-1,:) = 0.0_wp 1256 1257 ! 1258 !-- Top boundary at the outflow 1259 IF ( ibc_uv_t == 0 ) THEN 1260 u_p(nzt+1,-1,:) = u_init(nzt+1) 1261 v_p(nzt+1,0,:) = v_init(nzt+1) 1262 ELSE 1263 u_p(nzt+1,-1,:) = u_p(nzt,-1,:) 1264 v_p(nzt+1,0,:) = v_p(nzt,0,:) 1265 ENDIF 1266 w_p(nzt:nzt+1,-1,:) = 0.0_wp 1267 1268 ENDIF 1269 1132 u_p(:,-1,:) = u(:,0,:) 1133 v_p(:,0,:) = v(:,1,:) 1134 w_p(:,-1,:) = w(:,0,:) 1270 1135 ENDIF 1271 1136 1272 1137 IF ( bc_radiation_n ) THEN 1273 1274 IF ( use_cmax ) THEN 1275 u_p(:,ny+1,:) = u(:,ny,:) 1276 v_p(:,ny+1,:) = v(:,ny,:) 1277 w_p(:,ny+1,:) = w(:,ny,:) 1278 ELSEIF ( .NOT. use_cmax ) THEN 1279 1280 c_max = dy / dt_3d 1281 1282 c_u_m_l = 0.0_wp 1283 c_v_m_l = 0.0_wp 1284 c_w_m_l = 0.0_wp 1285 1286 c_u_m = 0.0_wp 1287 c_v_m = 0.0_wp 1288 c_w_m = 0.0_wp 1289 1290 ! 1291 !-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow 1292 !-- boundary. 1293 DO k = nzb+1, nzt+1 1294 DO i = nxl, nxr 1295 1296 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 1297 1298 IF ( denom /= 0.0_wp ) THEN 1299 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) ) 1300 IF ( c_u(k,i) < 0.0_wp ) THEN 1301 c_u(k,i) = 0.0_wp 1302 ELSEIF ( c_u(k,i) > c_max ) THEN 1303 c_u(k,i) = c_max 1304 ENDIF 1305 ELSE 1306 c_u(k,i) = c_max 1307 ENDIF 1308 1309 denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i) 1310 1311 IF ( denom /= 0.0_wp ) THEN 1312 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) ) 1313 IF ( c_v(k,i) < 0.0_wp ) THEN 1314 c_v(k,i) = 0.0_wp 1315 ELSEIF ( c_v(k,i) > c_max ) THEN 1316 c_v(k,i) = c_max 1317 ENDIF 1318 ELSE 1319 c_v(k,i) = c_max 1320 ENDIF 1321 1322 denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i) 1323 1324 IF ( denom /= 0.0_wp ) THEN 1325 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) ) 1326 IF ( c_w(k,i) < 0.0_wp ) THEN 1327 c_w(k,i) = 0.0_wp 1328 ELSEIF ( c_w(k,i) > c_max ) THEN 1329 c_w(k,i) = c_max 1330 ENDIF 1331 ELSE 1332 c_w(k,i) = c_max 1333 ENDIF 1334 1335 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 1336 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 1337 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 1338 1339 ENDDO 1340 ENDDO 1341 1342 #if defined( __parallel ) 1343 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1344 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, & 1345 ierr ) 1346 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1347 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, & 1348 ierr ) 1349 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1350 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, & 1351 ierr ) 1352 #else 1353 c_u_m = c_u_m_l 1354 c_v_m = c_v_m_l 1355 c_w_m = c_w_m_l 1356 #endif 1357 1358 c_u_m = c_u_m / (nx+1) 1359 c_v_m = c_v_m / (nx+1) 1360 c_w_m = c_w_m / (nx+1) 1361 1362 ! 1363 !-- Save old timelevels for the next timestep 1364 IF ( intermediate_timestep_count == 1 ) THEN 1365 u_m_n(:,:,:) = u(:,ny-1:ny,:) 1366 v_m_n(:,:,:) = v(:,ny-1:ny,:) 1367 w_m_n(:,:,:) = w(:,ny-1:ny,:) 1368 ENDIF 1369 1370 ! 1371 !-- Calculate the new velocities 1372 DO k = nzb+1, nzt+1 1373 DO i = nxlg, nxrg 1374 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * & 1375 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 1376 1377 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * & 1378 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 1379 1380 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * & 1381 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 1382 ENDDO 1383 ENDDO 1384 1385 ! 1386 !-- Bottom boundary at the outflow 1387 IF ( ibc_uv_b == 0 ) THEN 1388 u_p(nzb,ny+1,:) = 0.0_wp 1389 v_p(nzb,ny+1,:) = 0.0_wp 1390 ELSE 1391 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 1392 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) 1393 ENDIF 1394 w_p(nzb,ny+1,:) = 0.0_wp 1395 1396 ! 1397 !-- Top boundary at the outflow 1398 IF ( ibc_uv_t == 0 ) THEN 1399 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 1400 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 1401 ELSE 1402 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) 1403 v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) 1404 ENDIF 1405 w_p(nzt:nzt+1,ny+1,:) = 0.0_wp 1406 1407 ENDIF 1408 1138 u_p(:,ny+1,:) = u(:,ny,:) 1139 v_p(:,ny+1,:) = v(:,ny,:) 1140 w_p(:,ny+1,:) = w(:,ny,:) 1409 1141 ENDIF 1410 1142 1411 1143 IF ( bc_radiation_l ) THEN 1412 1413 IF ( use_cmax ) THEN 1414 u_p(:,:,0) = u(:,:,1) 1415 v_p(:,:,-1) = v(:,:,0) 1416 w_p(:,:,-1) = w(:,:,0) 1417 ELSEIF ( .NOT. use_cmax ) THEN 1418 1419 c_max = dx / dt_3d 1420 1421 c_u_m_l = 0.0_wp 1422 c_v_m_l = 0.0_wp 1423 c_w_m_l = 0.0_wp 1424 1425 c_u_m = 0.0_wp 1426 c_v_m = 0.0_wp 1427 c_w_m = 0.0_wp 1428 1429 ! 1430 !-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow 1431 !-- boundary. 1432 DO k = nzb+1, nzt+1 1433 DO j = nys, nyn 1434 1435 denom = u_m_l(k,j,1) - u_m_l(k,j,2) 1436 1437 IF ( denom /= 0.0_wp ) THEN 1438 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) ) 1439 IF ( c_u(k,j) < 0.0_wp ) THEN 1440 c_u(k,j) = 0.0_wp 1441 ELSEIF ( c_u(k,j) > c_max ) THEN 1442 c_u(k,j) = c_max 1443 ENDIF 1444 ELSE 1445 c_u(k,j) = c_max 1446 ENDIF 1447 1448 denom = v_m_l(k,j,0) - v_m_l(k,j,1) 1449 1450 IF ( denom /= 0.0_wp ) THEN 1451 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) ) 1452 IF ( c_v(k,j) < 0.0_wp ) THEN 1453 c_v(k,j) = 0.0_wp 1454 ELSEIF ( c_v(k,j) > c_max ) THEN 1455 c_v(k,j) = c_max 1456 ENDIF 1457 ELSE 1458 c_v(k,j) = c_max 1459 ENDIF 1460 1461 denom = w_m_l(k,j,0) - w_m_l(k,j,1) 1462 1463 IF ( denom /= 0.0_wp ) THEN 1464 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) ) 1465 IF ( c_w(k,j) < 0.0_wp ) THEN 1466 c_w(k,j) = 0.0_wp 1467 ELSEIF ( c_w(k,j) > c_max ) THEN 1468 c_w(k,j) = c_max 1469 ENDIF 1470 ELSE 1471 c_w(k,j) = c_max 1472 ENDIF 1473 1474 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 1475 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 1476 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 1477 1478 ENDDO 1479 ENDDO 1480 1481 #if defined( __parallel ) 1482 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1483 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, & 1484 ierr ) 1485 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1486 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, & 1487 ierr ) 1488 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1489 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, & 1490 ierr ) 1491 #else 1492 c_u_m = c_u_m_l 1493 c_v_m = c_v_m_l 1494 c_w_m = c_w_m_l 1495 #endif 1496 1497 c_u_m = c_u_m / (ny+1) 1498 c_v_m = c_v_m / (ny+1) 1499 c_w_m = c_w_m / (ny+1) 1500 1501 ! 1502 !-- Save old timelevels for the next timestep 1503 IF ( intermediate_timestep_count == 1 ) THEN 1504 u_m_l(:,:,:) = u(:,:,1:2) 1505 v_m_l(:,:,:) = v(:,:,0:1) 1506 w_m_l(:,:,:) = w(:,:,0:1) 1507 ENDIF 1508 1509 ! 1510 !-- Calculate the new velocities 1511 DO k = nzb+1, nzt+1 1512 DO j = nysg, nyng 1513 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * & 1514 ( u(k,j,0) - u(k,j,1) ) * ddx 1515 1516 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * & 1517 ( v(k,j,-1) - v(k,j,0) ) * ddx 1518 1519 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * & 1520 ( w(k,j,-1) - w(k,j,0) ) * ddx 1521 ENDDO 1522 ENDDO 1523 1524 ! 1525 !-- Bottom boundary at the outflow 1526 IF ( ibc_uv_b == 0 ) THEN 1527 u_p(nzb,:,0) = 0.0_wp 1528 v_p(nzb,:,-1) = 0.0_wp 1529 ELSE 1530 u_p(nzb,:,0) = u_p(nzb+1,:,0) 1531 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 1532 ENDIF 1533 w_p(nzb,:,-1) = 0.0_wp 1534 1535 ! 1536 !-- Top boundary at the outflow 1537 IF ( ibc_uv_t == 0 ) THEN 1538 u_p(nzt+1,:,0) = u_init(nzt+1) 1539 v_p(nzt+1,:,-1) = v_init(nzt+1) 1540 ELSE 1541 u_p(nzt+1,:,0) = u_p(nzt,:,0) 1542 v_p(nzt+1,:,-1) = v_p(nzt,:,-1) 1543 ENDIF 1544 w_p(nzt:nzt+1,:,-1) = 0.0_wp 1545 1546 ENDIF 1547 1144 u_p(:,:,0) = u(:,:,1) 1145 v_p(:,:,-1) = v(:,:,0) 1146 w_p(:,:,-1) = w(:,:,0) 1548 1147 ENDIF 1549 1148 1550 1149 IF ( bc_radiation_r ) THEN 1551 1552 IF ( use_cmax ) THEN 1553 u_p(:,:,nx+1) = u(:,:,nx) 1554 v_p(:,:,nx+1) = v(:,:,nx) 1555 w_p(:,:,nx+1) = w(:,:,nx) 1556 ELSEIF ( .NOT. use_cmax ) THEN 1557 1558 c_max = dx / dt_3d 1559 1560 c_u_m_l = 0.0_wp 1561 c_v_m_l = 0.0_wp 1562 c_w_m_l = 0.0_wp 1563 1564 c_u_m = 0.0_wp 1565 c_v_m = 0.0_wp 1566 c_w_m = 0.0_wp 1567 1568 ! 1569 !-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow 1570 !-- boundary. 1571 DO k = nzb+1, nzt+1 1572 DO j = nys, nyn 1573 1574 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) 1575 1576 IF ( denom /= 0.0_wp ) THEN 1577 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) ) 1578 IF ( c_u(k,j) < 0.0_wp ) THEN 1579 c_u(k,j) = 0.0_wp 1580 ELSEIF ( c_u(k,j) > c_max ) THEN 1581 c_u(k,j) = c_max 1582 ENDIF 1583 ELSE 1584 c_u(k,j) = c_max 1585 ENDIF 1586 1587 denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1) 1588 1589 IF ( denom /= 0.0_wp ) THEN 1590 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) ) 1591 IF ( c_v(k,j) < 0.0_wp ) THEN 1592 c_v(k,j) = 0.0_wp 1593 ELSEIF ( c_v(k,j) > c_max ) THEN 1594 c_v(k,j) = c_max 1595 ENDIF 1596 ELSE 1597 c_v(k,j) = c_max 1598 ENDIF 1599 1600 denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1) 1601 1602 IF ( denom /= 0.0_wp ) THEN 1603 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) ) 1604 IF ( c_w(k,j) < 0.0_wp ) THEN 1605 c_w(k,j) = 0.0_wp 1606 ELSEIF ( c_w(k,j) > c_max ) THEN 1607 c_w(k,j) = c_max 1608 ENDIF 1609 ELSE 1610 c_w(k,j) = c_max 1611 ENDIF 1612 1613 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 1614 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 1615 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 1616 1617 ENDDO 1618 ENDDO 1619 1620 #if defined( __parallel ) 1621 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1622 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, & 1623 ierr ) 1624 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1625 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, & 1626 ierr ) 1627 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1628 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, & 1629 ierr ) 1630 #else 1631 c_u_m = c_u_m_l 1632 c_v_m = c_v_m_l 1633 c_w_m = c_w_m_l 1634 #endif 1635 1636 c_u_m = c_u_m / (ny+1) 1637 c_v_m = c_v_m / (ny+1) 1638 c_w_m = c_w_m / (ny+1) 1639 1640 ! 1641 !-- Save old timelevels for the next timestep 1642 IF ( intermediate_timestep_count == 1 ) THEN 1643 u_m_r(:,:,:) = u(:,:,nx-1:nx) 1644 v_m_r(:,:,:) = v(:,:,nx-1:nx) 1645 w_m_r(:,:,:) = w(:,:,nx-1:nx) 1646 ENDIF 1647 1648 ! 1649 !-- Calculate the new velocities 1650 DO k = nzb+1, nzt+1 1651 DO j = nysg, nyng 1652 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * & 1653 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 1654 1655 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * & 1656 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 1657 1658 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * & 1659 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 1660 ENDDO 1661 ENDDO 1662 1663 ! 1664 !-- Bottom boundary at the outflow 1665 IF ( ibc_uv_b == 0 ) THEN 1666 u_p(nzb,:,nx+1) = 0.0_wp 1667 v_p(nzb,:,nx+1) = 0.0_wp 1668 ELSE 1669 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 1670 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) 1671 ENDIF 1672 w_p(nzb,:,nx+1) = 0.0_wp 1673 1674 ! 1675 !-- Top boundary at the outflow 1676 IF ( ibc_uv_t == 0 ) THEN 1677 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 1678 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 1679 ELSE 1680 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) 1681 v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) 1682 ENDIF 1683 w_p(nzt:nzt+1,:,nx+1) = 0.0_wp 1684 1685 ENDIF 1686 1150 u_p(:,:,nx+1) = u(:,:,nx) 1151 v_p(:,:,nx+1) = v(:,:,nx) 1152 w_p(:,:,nx+1) = w(:,:,nx) 1687 1153 ENDIF 1688 1154 1689 1155 END SUBROUTINE dynamics_boundary_conditions 1156 1157 1690 1158 !--------------------------------------------------------------------------------------------------! 1691 1159 ! Description: … … 1780 1248 1781 1249 ENDIF 1782 1783 1250 1784 1251 END SUBROUTINE dynamics_3d_data_averaging
Note: See TracChangeset
for help on using the changeset viewer.
