Changeset 1797 for palm/trunk/SOURCE/pmc_interface.f90
- Timestamp:
- Mar 21, 2016 4:50:28 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/pmc_interface.f90 (modified) (87 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pmc_interface.f90
r1792 r1797 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! introduction of different datatransfer modes, 23 ! further formatting cleanup, parameter checks added (including mismatches 24 ! between root and client model settings), 25 ! +routine pmci_check_setting_mismatches 26 ! comm_world_nesting introduced 23 27 ! 24 28 ! Former revisions: … … 145 149 ! 146 150 !-- Coupler setup 151 INTEGER(iwp), SAVE :: comm_world_nesting !: 147 152 INTEGER(iwp), SAVE :: cpl_id = 1 !: 148 153 CHARACTER(LEN=32), SAVE :: cpl_name !: … … 152 157 ! 153 158 !-- Control parameters, will be made input parameters later 159 CHARACTER(LEN=7), SAVE :: nesting_datatransfer_mode = 'mixed' !: steering 160 !: parameter for data- 161 !: transfer mode 154 162 CHARACTER(LEN=7), SAVE :: nesting_mode = 'two-way' !: steering parameter 155 163 !: for 1- or 2-way nesting … … 308 316 309 317 310 INTERFACE pmci_c lient_datatrans311 MODULE PROCEDURE pmci_c lient_datatrans318 INTERFACE pmci_check_setting_mismatches 319 MODULE PROCEDURE pmci_check_setting_mismatches 312 320 END INTERFACE 313 321 … … 320 328 END INTERFACE 321 329 330 INTERFACE pmci_datatrans 331 MODULE PROCEDURE pmci_datatrans 332 END INTERFACE pmci_datatrans 333 322 334 INTERFACE pmci_ensure_nest_mass_conservation 323 335 MODULE PROCEDURE pmci_ensure_nest_mass_conservation … … 346 358 PUBLIC anterp_relax_length_l, anterp_relax_length_r, & 347 359 anterp_relax_length_s, anterp_relax_length_n, & 348 anterp_relax_length_t, client_to_server, c pl_id, nested_run, &349 nesting_mode, server_to_client350 PUBLIC pmci_client_datatrans360 anterp_relax_length_t, client_to_server, comm_world_nesting, & 361 cpl_id, nested_run, nesting_datatransfer_mode, nesting_mode, & 362 server_to_client 351 363 PUBLIC pmci_client_initialize 352 364 PUBLIC pmci_client_synchronize 365 PUBLIC pmci_datatrans 353 366 PUBLIC pmci_ensure_nest_mass_conservation 354 367 PUBLIC pmci_init 355 368 PUBLIC pmci_modelconfiguration 356 PUBLIC pmci_server_datatrans357 369 PUBLIC pmci_server_initialize 358 370 PUBLIC pmci_server_synchronize … … 365 377 SUBROUTINE pmci_init( world_comm ) 366 378 379 USE control_parameters, & 380 ONLY: message_string 381 367 382 IMPLICIT NONE 368 383 … … 376 391 377 392 378 CALL pmc_init_model( world_comm, nesting_mode, pmc_status ) 393 CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode, & 394 pmc_status ) 379 395 380 396 IF ( pmc_status == pmc_no_namelist_found ) THEN … … 390 406 391 407 ! 408 !-- Check steering parameter values 409 IF ( TRIM( nesting_mode ) /= 'one-way' .AND. & 410 TRIM( nesting_mode ) /= 'two-way' ) & 411 THEN 412 message_string = 'illegal nesting mode: ' // TRIM( nesting_mode ) 413 CALL message( 'pmci_init', 'PA0417', 3, 2, 0, 6, 0 ) 414 ENDIF 415 416 IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade' .AND. & 417 TRIM( nesting_datatransfer_mode ) /= 'mixed' .AND. & 418 TRIM( nesting_datatransfer_mode ) /= 'overlap' ) & 419 THEN 420 message_string = 'illegal nesting datatransfer mode: ' & 421 // TRIM( nesting_datatransfer_mode ) 422 CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 ) 423 ENDIF 424 425 ! 392 426 !-- Set the general steering switch which tells PALM that its a nested run 393 427 nested_run = .TRUE. … … 396 430 !-- Get some variables required by the pmc-interface (and in some cases in the 397 431 !-- PALM code out of the pmci) out of the pmc-core 398 CALL pmc_get_model_info( cpl_id = cpl_id, cpl_parent_id = cpl_parent_id, & 432 CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting, & 433 cpl_id = cpl_id, cpl_parent_id = cpl_parent_id, & 399 434 cpl_name = cpl_name, npe_total = cpl_npe_total, & 400 435 lower_left_x = lower_left_coord_x, & … … 445 480 !-- Initialize PMC Server 446 481 CALL pmci_setup_server 482 ! 483 !-- Check for mismatches between settings of master and client variables 484 !-- (e.g., all clients have to follow the end_time settings of the root master) 485 CALL pmci_check_setting_mismatches 447 486 448 487 CALL location_message( 'finished', .TRUE. ) … … 701 740 ELSE 702 741 ! 703 !-- TO_DO: Klaus: comment why thi edummy allocation is required742 !-- TO_DO: Klaus: comment why this dummy allocation is required 704 743 ALLOCATE( index_list(6,1) ) 705 744 CALL pmc_s_set_2d_index_list( client_id, index_list ) … … 789 828 SIZE(define_coarse_grid_real), 0, 21, ierr ) 790 829 CALL pmc_recv_from_server( define_coarse_grid_int, 3, 0, 22, ierr ) 791 792 !793 !-- Receive also the dz-,zu- and zw-arrays here.794 !-- TO_DO: what is the meaning of above comment795 830 ! 796 831 !-- Debug-printouts - keep them … … 852 887 853 888 ! 854 !-- TO_DO: give comments what is happening here889 !-- Find the index bounds for the nest domain in the coarse-grid index space 855 890 CALL pmci_map_fine_to_coarse_grid 891 ! 892 !-- TO_DO: Klaus give a comment what is happening here 856 893 CALL pmc_c_get_2d_index_list 857 894 … … 879 916 880 917 ! 881 !-- Two-way coupling 882 !-- TO_DO: comment what is happening here 918 !-- Two-way coupling. 919 !-- Precompute the index arrays and relaxation functions for the 920 !-- anterpolation 883 921 IF ( nesting_mode == 'two-way' ) THEN 884 922 CALL pmci_init_anterp_tophat … … 895 933 896 934 SUBROUTINE pmci_map_fine_to_coarse_grid 897 935 ! 936 !-- Determine index bounds of interpolation/anterpolation area in the coarse 937 !-- grid index space 898 938 IMPLICIT NONE 899 939 … … 907 947 908 948 ! 909 !-- Determine indices of interpolation/anterpolation area in the coarse grid 910 !-- If the fine- and coarse grid nodes do not match. 949 !-- If the fine- and coarse grid nodes do not match: 911 950 loffset = MOD( coord_x(nxl), cg%dx ) 912 951 xexl = cg%dx + loffset … … 991 1030 !-- Precomputation of the interpolation coefficients and client-array indices 992 1031 !-- to be used by the interpolation routines interp_tril_lr, interp_tril_ns 993 !-- and interp_tril_t. Constant dz is still assumed.1032 !-- and interp_tril_t. 994 1033 995 1034 IMPLICIT NONE … … 1199 1238 !-- Left boundary 1200 1239 IF ( nest_bound_l ) THEN 1240 1201 1241 ALLOCATE( logc_u_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) ) 1202 1242 ALLOCATE( logc_v_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) ) … … 1235 1275 1236 1276 ENDDO 1277 1237 1278 ENDIF 1238 1279 … … 1240 1281 !-- Right boundary 1241 1282 IF ( nest_bound_r ) THEN 1283 1242 1284 ALLOCATE( logc_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) ) 1243 1285 ALLOCATE( logc_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) ) … … 1252 1294 DO j = nys, nyn 1253 1295 ! 1254 !-- Right boundary for u .1296 !-- Right boundary for u 1255 1297 i = nxr + 1 1256 1298 kb = nzb_u_inner(j,i) … … 1262 1304 logc_ratio_u_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1) 1263 1305 lcr(0:ncorr-1) = 1.0_wp 1264 1265 ! 1266 !-- Right boundary for v. 1306 ! 1307 !-- Right boundary for v 1267 1308 i = nxr + 1 1268 1309 kb = nzb_v_inner(j,i) … … 1276 1317 1277 1318 ENDDO 1319 1278 1320 ENDIF 1279 1321 … … 1281 1323 !-- South boundary 1282 1324 IF ( nest_bound_s ) THEN 1325 1283 1326 ALLOCATE( logc_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) ) 1284 1327 ALLOCATE( logc_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) ) … … 1291 1334 direction = 1 1292 1335 inc = 1 1336 1293 1337 DO i = nxl, nxr 1294 1338 ! 1295 !-- South boundary for u .1339 !-- South boundary for u 1296 1340 j = -1 1297 1341 kb = nzb_u_inner(j,i) … … 1303 1347 logc_ratio_u_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1) 1304 1348 lcr(0:ncorr-1) = 1.0_wp 1305 1306 1349 ! 1307 1350 !-- South boundary for v … … 1315 1358 logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1) 1316 1359 lcr(0:ncorr-1) = 1.0_wp 1317 ENDDO 1360 1361 ENDDO 1362 1318 1363 ENDIF 1319 1364 … … 1321 1366 !-- North boundary 1322 1367 IF ( nest_bound_n ) THEN 1368 1323 1369 ALLOCATE( logc_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) ) 1324 1370 ALLOCATE( logc_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) ) … … 1331 1377 direction = 1 1332 1378 inc = 1 1379 1333 1380 DO i = nxl, nxr 1334 1381 ! 1335 !-- North boundary for u .1382 !-- North boundary for u 1336 1383 j = nyn + 1 1337 1384 kb = nzb_u_inner(j,i) … … 1343 1390 logc_ratio_u_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1) 1344 1391 lcr(0:ncorr-1) = 1.0_wp 1345 1346 ! 1347 !-- North boundary for v. 1392 ! 1393 !-- North boundary for v 1348 1394 j = nyn + 1 1349 1395 kb = nzb_v_inner(j,i) … … 1357 1403 1358 1404 ENDDO 1405 1359 1406 ENDIF 1360 1407 1361 1408 ! 1362 !-- Then vertical walls and corners if necessary .1409 !-- Then vertical walls and corners if necessary 1363 1410 IF ( topography /= 'flat' ) THEN 1411 1364 1412 kb = 0 ! kb is not used when direction > 1 1365 1413 ! 1366 1414 !-- Left boundary 1367 1415 IF ( nest_bound_l ) THEN 1416 1368 1417 ALLOCATE( logc_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) ) 1369 1418 ALLOCATE( logc_ratio_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2, & … … 1374 1423 DO j = nys, nyn 1375 1424 DO k = nzb, nzt_topo_nestbc_l 1376 1377 1425 ! 1378 1426 !-- Wall for u on the south side, but not on the north side … … 1392 1440 lcr(0:ncorr-1) = 1.0_wp 1393 1441 ENDIF 1442 1394 1443 ! 1395 1444 !-- Wall for u on the north side, but not on the south side 1396 1445 i = 0 1397 !-- TO_DO: routine must be indentet by 1 space from here on, 1398 !-- and long lines must be wrapped 1399 IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND. & 1400 ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) ) THEN 1401 inc = -1 1402 wall_index = j + 1 1403 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1404 1405 ! 1406 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1407 logc_u_l(k,j,2) = inc * lc 1408 logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1409 lcr(0:ncorr-1) = 1.0_wp 1410 ENDIF 1411 1412 ! 1413 !-- Wall for w on the south side, but not on the north side. 1414 i = -1 1415 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) .AND. & 1416 ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) ) THEN 1417 inc = 1 1418 wall_index = j 1419 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1420 1421 ! 1422 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1423 logc_w_l(k,j,2) = inc * lc 1424 logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1425 lcr(0:ncorr-1) = 1.0_wp 1426 ENDIF 1427 1428 ! 1429 !-- Wall for w on the north side, but not on the south side. 1430 i = -1 1431 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND. & 1432 ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) ) THEN 1433 inc = -1 1434 wall_index = j+1 1435 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1436 1437 ! 1438 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1439 logc_w_l(k,j,2) = inc * lc 1440 logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1441 lcr(0:ncorr-1) = 1.0_wp 1442 ENDIF 1443 ENDDO 1444 ENDDO 1445 ENDIF ! IF ( nest_bound_l ) 1446 IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND. & 1447 ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) ) THEN 1448 inc = -1 1449 wall_index = j + 1 1450 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1451 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1452 ! 1453 !-- The direction of the wall-normal index is stored as the 1454 !-- sign of the logc-element. 1455 logc_u_l(k,j,2) = inc * lc 1456 logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1457 lcr(0:ncorr-1) = 1.0_wp 1458 ENDIF 1459 1460 ! 1461 !-- Wall for w on the south side, but not on the north side. 1462 i = -1 1463 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) .AND. & 1464 ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) ) THEN 1465 inc = 1 1466 wall_index = j 1467 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1468 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1469 ! 1470 !-- The direction of the wall-normal index is stored as the 1471 !-- sign of the logc-element. 1472 logc_w_l(k,j,2) = inc * lc 1473 logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1474 lcr(0:ncorr-1) = 1.0_wp 1475 ENDIF 1476 1477 ! 1478 !-- Wall for w on the north side, but not on the south side. 1479 i = -1 1480 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND. & 1481 ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) ) THEN 1482 inc = -1 1483 wall_index = j+1 1484 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1485 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1486 ! 1487 !-- The direction of the wall-normal index is stored as the 1488 !-- sign of the logc-element. 1489 logc_w_l(k,j,2) = inc * lc 1490 logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1491 lcr(0:ncorr-1) = 1.0_wp 1492 ENDIF 1493 1494 ENDDO 1495 ENDDO 1496 1497 ENDIF ! IF ( nest_bound_l ) 1446 1498 1447 1499 ! 1448 !-- Right boundary. 1449 IF ( nest_bound_r ) THEN 1450 ALLOCATE( logc_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) ) 1451 ALLOCATE( logc_ratio_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) ) 1452 logc_w_r = 0 1453 logc_ratio_w_r = 1.0_wp 1454 direction = 2 1455 i = nxr + 1 1456 DO j = nys, nyn 1457 DO k = nzb, nzt_topo_nestbc_r 1458 1459 ! 1460 !-- Wall for u on the south side, but not on the north side. 1461 IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND. & 1462 ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) ) THEN 1463 inc = 1 1464 wall_index = j 1465 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1466 1467 ! 1468 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1469 logc_u_r(k,j,2) = inc * lc 1470 logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1471 lcr(0:ncorr-1) = 1.0_wp 1472 ENDIF 1473 1474 ! 1475 !-- Wall for u on the north side, but not on the south side. 1476 IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND. & 1477 ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) ) THEN 1478 inc = -1 1479 wall_index = j+1 1480 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1481 1482 ! 1483 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1484 logc_u_r(k,j,2) = inc * lc 1485 logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1486 lcr(0:ncorr-1) = 1.0_wp 1487 ENDIF 1488 1489 ! 1490 !-- Wall for w on the south side, but not on the north side. 1491 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) .AND. & 1492 ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) ) THEN 1493 inc = 1 1494 wall_index = j 1495 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1496 1497 ! 1498 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1499 logc_w_r(k,j,2) = inc * lc 1500 logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1501 lcr(0:ncorr-1) = 1.0_wp 1502 ENDIF 1503 ! 1504 !-- Wall for w on the north side, but not on the south side. 1505 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND. & 1506 ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) ) THEN 1507 inc = -1 1508 wall_index = j+1 1509 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1510 1511 ! 1512 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1513 logc_w_r(k,j,2) = inc * lc 1514 logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1515 lcr(0:ncorr-1) = 1.0_wp 1516 ENDIF 1517 ENDDO 1518 ENDDO 1519 ENDIF ! IF ( nest_bound_r ) 1500 !-- Right boundary 1501 IF ( nest_bound_r ) THEN 1502 1503 ALLOCATE( logc_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) ) 1504 ALLOCATE( logc_ratio_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2, & 1505 0:ncorr-1) ) 1506 logc_w_r = 0 1507 logc_ratio_w_r = 1.0_wp 1508 direction = 2 1509 i = nxr + 1 1510 1511 DO j = nys, nyn 1512 DO k = nzb, nzt_topo_nestbc_r 1513 ! 1514 !-- Wall for u on the south side, but not on the north side 1515 IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND. & 1516 ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) ) THEN 1517 inc = 1 1518 wall_index = j 1519 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1520 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1521 ! 1522 !-- The direction of the wall-normal index is stored as the 1523 !-- sign of the logc-element. 1524 logc_u_r(k,j,2) = inc * lc 1525 logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1526 lcr(0:ncorr-1) = 1.0_wp 1527 ENDIF 1528 1529 ! 1530 !-- Wall for u on the north side, but not on the south side 1531 IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND. & 1532 ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) ) THEN 1533 inc = -1 1534 wall_index = j+1 1535 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1536 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1537 ! 1538 !-- The direction of the wall-normal index is stored as the 1539 !-- sign of the logc-element. 1540 logc_u_r(k,j,2) = inc * lc 1541 logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1542 lcr(0:ncorr-1) = 1.0_wp 1543 ENDIF 1544 1545 ! 1546 !-- Wall for w on the south side, but not on the north side 1547 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) .AND. & 1548 ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) ) THEN 1549 inc = 1 1550 wall_index = j 1551 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1552 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1553 ! 1554 !-- The direction of the wall-normal index is stored as the 1555 !-- sign of the logc-element. 1556 logc_w_r(k,j,2) = inc * lc 1557 logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1558 lcr(0:ncorr-1) = 1.0_wp 1559 ENDIF 1560 1561 ! 1562 !-- Wall for w on the north side, but not on the south side 1563 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND. & 1564 ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) ) THEN 1565 inc = -1 1566 wall_index = j+1 1567 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1568 k, j, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1569 1570 ! 1571 !-- The direction of the wall-normal index is stored as the 1572 !-- sign of the logc-element. 1573 logc_w_r(k,j,2) = inc * lc 1574 logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1) 1575 lcr(0:ncorr-1) = 1.0_wp 1576 ENDIF 1577 1578 ENDDO 1579 ENDDO 1580 1581 ENDIF ! IF ( nest_bound_r ) 1520 1582 1521 1583 ! 1522 !-- South boundary. 1523 IF ( nest_bound_s ) THEN 1524 ALLOCATE( logc_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2) ) 1525 ALLOCATE( logc_ratio_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2, 0:ncorr-1) ) 1526 logc_w_s = 0 1527 logc_ratio_w_s = 1.0_wp 1528 direction = 3 1529 DO i = nxl, nxr 1530 DO k = nzb, nzt_topo_nestbc_s 1531 1532 ! 1533 !-- Wall for v on the left side, but not on the right side. 1534 j = 0 1535 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) .AND. & 1536 ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) ) THEN 1537 inc = 1 1538 wall_index = i 1539 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1540 1541 ! 1542 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1543 logc_v_s(k,i,2) = inc * lc 1544 logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1545 lcr(0:ncorr-1) = 1.0_wp 1546 ENDIF 1547 ! 1548 !-- Wall for v on the right side, but not on the left side. 1549 j = 0 1550 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) .AND. & 1551 ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) ) THEN 1552 inc = -1 1553 wall_index = i+1 1554 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1555 1556 ! 1557 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1558 logc_v_s(k,i,2) = inc * lc 1559 logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1560 lcr(0:ncorr-1) = 1.0_wp 1561 ENDIF 1562 1563 ! 1564 !-- Wall for w on the left side, but not on the right side. 1565 j = -1 1566 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) .AND. & 1567 ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) ) THEN 1568 inc = 1 1569 wall_index = i 1570 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1571 1572 ! 1573 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1574 logc_w_s(k,i,2) = inc * lc 1575 logc_ratio_w_s(k,i,2,0:ncorr - 1) = lcr(0:ncorr-1) 1576 lcr(0:ncorr-1) = 1.0_wp 1577 ENDIF 1578 1579 ! 1580 !-- Wall for w on the right side, but not on the left side. 1581 j = -1 1582 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) .AND. & 1583 ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) ) THEN 1584 inc = -1 1585 wall_index = i+1 1586 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1587 1588 ! 1589 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1590 logc_w_s(k,i,2) = inc * lc 1591 logc_ratio_w_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1592 lcr(0:ncorr-1) = 1.0_wp 1593 ENDIF 1594 ENDDO 1595 ENDDO 1596 ENDIF ! IF (nest_bound_s ) 1584 !-- South boundary 1585 IF ( nest_bound_s ) THEN 1586 1587 ALLOCATE( logc_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2) ) 1588 ALLOCATE( logc_ratio_w_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2, & 1589 0:ncorr-1) ) 1590 logc_w_s = 0 1591 logc_ratio_w_s = 1.0_wp 1592 direction = 3 1593 1594 DO i = nxl, nxr 1595 DO k = nzb, nzt_topo_nestbc_s 1596 ! 1597 !-- Wall for v on the left side, but not on the right side 1598 j = 0 1599 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) .AND. & 1600 ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) ) THEN 1601 inc = 1 1602 wall_index = i 1603 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1604 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1605 ! 1606 !-- The direction of the wall-normal index is stored as the 1607 !-- sign of the logc-element. 1608 logc_v_s(k,i,2) = inc * lc 1609 logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1610 lcr(0:ncorr-1) = 1.0_wp 1611 ENDIF 1612 1613 ! 1614 !-- Wall for v on the right side, but not on the left side 1615 j = 0 1616 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) .AND. & 1617 ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) ) THEN 1618 inc = -1 1619 wall_index = i+1 1620 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1621 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1622 ! 1623 !-- The direction of the wall-normal index is stored as the 1624 !-- sign of the logc-element. 1625 logc_v_s(k,i,2) = inc * lc 1626 logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1627 lcr(0:ncorr-1) = 1.0_wp 1628 ENDIF 1629 1630 ! 1631 !-- Wall for w on the left side, but not on the right side 1632 j = -1 1633 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) .AND. & 1634 ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) ) THEN 1635 inc = 1 1636 wall_index = i 1637 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1638 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1639 ! 1640 !-- The direction of the wall-normal index is stored as the 1641 !-- sign of the logc-element. 1642 logc_w_s(k,i,2) = inc * lc 1643 logc_ratio_w_s(k,i,2,0:ncorr - 1) = lcr(0:ncorr-1) 1644 lcr(0:ncorr-1) = 1.0_wp 1645 ENDIF 1646 1647 ! 1648 !-- Wall for w on the right side, but not on the left side 1649 j = -1 1650 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) .AND. & 1651 ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) ) THEN 1652 inc = -1 1653 wall_index = i+1 1654 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1655 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1656 ! 1657 !-- The direction of the wall-normal index is stored as the 1658 !-- sign of the logc-element. 1659 logc_w_s(k,i,2) = inc * lc 1660 logc_ratio_w_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1661 lcr(0:ncorr-1) = 1.0_wp 1662 ENDIF 1663 1664 ENDDO 1665 ENDDO 1666 1667 ENDIF ! IF (nest_bound_s ) 1597 1668 1598 1669 ! 1599 !-- North boundary. 1600 IF ( nest_bound_n ) THEN 1601 ALLOCATE( logc_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2) ) 1602 ALLOCATE( logc_ratio_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2, 0:ncorr-1) ) 1603 logc_w_n = 0 1604 logc_ratio_w_n = 1.0_wp 1605 direction = 3 1606 j = nyn + 1 1607 DO i = nxl, nxr 1608 DO k = nzb, nzt_topo_nestbc_n 1609 1610 ! 1611 !-- Wall for v on the left side, but not on the right side. 1612 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) .AND. & 1613 ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) ) THEN 1614 inc = 1 1615 wall_index = i 1616 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1617 1618 ! 1619 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1620 logc_v_n(k,i,2) = inc * lc 1621 logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1622 lcr(0:ncorr-1) = 1.0_wp 1623 ENDIF 1624 1625 ! 1626 !-- Wall for v on the right side, but not on the left side. 1627 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) .AND. & 1628 ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) ) THEN 1629 inc = -1 1630 wall_index = i + 1 1631 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1632 1633 ! 1634 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1635 logc_v_n(k,i,2) = inc * lc 1636 logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1637 lcr(0:ncorr-1) = 1.0_wp 1638 ENDIF 1639 1640 ! 1641 !-- Wall for w on the left side, but not on the right side. 1642 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) .AND. & 1643 ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) ) THEN 1644 inc = 1 1645 wall_index = i 1646 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1647 1648 ! 1649 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1650 logc_w_n(k,i,2) = inc * lc 1651 logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1652 lcr(0:ncorr-1) = 1.0_wp 1653 ENDIF 1654 1655 ! 1656 !-- Wall for w on the right side, but not on the left side. 1657 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) .AND. & 1658 ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) ) THEN 1659 inc = -1 1660 wall_index = i+1 1661 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1662 1663 ! 1664 !-- The direction of the wall-normal index is stored as the sign of the logc-element. 1665 logc_w_n(k,i,2) = inc * lc 1666 logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1667 lcr(0:ncorr-1) = 1.0_wp 1668 ENDIF 1669 ENDDO 1670 ENDDO 1671 ENDIF ! IF ( nest_bound_n ) 1672 ENDIF ! IF ( topography /= 'flat' ) 1673 1674 END SUBROUTINE pmci_init_loglaw_correction 1670 !-- North boundary 1671 IF ( nest_bound_n ) THEN 1672 1673 ALLOCATE( logc_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2) ) 1674 ALLOCATE( logc_ratio_w_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2, & 1675 0:ncorr-1) ) 1676 logc_w_n = 0 1677 logc_ratio_w_n = 1.0_wp 1678 direction = 3 1679 j = nyn + 1 1680 1681 DO i = nxl, nxr 1682 DO k = nzb, nzt_topo_nestbc_n 1683 ! 1684 !-- Wall for v on the left side, but not on the right side 1685 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) .AND. & 1686 ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) ) THEN 1687 inc = 1 1688 wall_index = i 1689 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1690 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1691 ! 1692 !-- The direction of the wall-normal index is stored as the 1693 !-- sign of the logc-element. 1694 logc_v_n(k,i,2) = inc * lc 1695 logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1696 lcr(0:ncorr-1) = 1.0_wp 1697 ENDIF 1698 1699 ! 1700 !-- Wall for v on the right side, but not on the left side 1701 IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) .AND. & 1702 ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) ) THEN 1703 inc = -1 1704 wall_index = i + 1 1705 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1706 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1707 ! 1708 !-- The direction of the wall-normal index is stored as the 1709 !-- sign of the logc-element. 1710 logc_v_n(k,i,2) = inc * lc 1711 logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1712 lcr(0:ncorr-1) = 1.0_wp 1713 ENDIF 1714 1715 ! 1716 !-- Wall for w on the left side, but not on the right side 1717 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) .AND. & 1718 ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) ) THEN 1719 inc = 1 1720 wall_index = i 1721 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1722 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1723 ! 1724 !-- The direction of the wall-normal index is stored as the 1725 !-- sign of the logc-element. 1726 logc_w_n(k,i,2) = inc * lc 1727 logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1728 lcr(0:ncorr-1) = 1.0_wp 1729 ENDIF 1730 1731 ! 1732 !-- Wall for w on the right side, but not on the left side 1733 IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) .AND. & 1734 ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) ) THEN 1735 inc = -1 1736 wall_index = i+1 1737 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1738 k, i, inc, wall_index, z0(j,i), kb, direction, ncorr ) 1739 ! 1740 !-- The direction of the wall-normal index is stored as the 1741 !-- sign of the logc-element. 1742 logc_w_n(k,i,2) = inc * lc 1743 logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1) 1744 lcr(0:ncorr-1) = 1.0_wp 1745 ENDIF 1746 1747 ENDDO 1748 ENDDO 1749 1750 ENDIF ! IF ( nest_bound_n ) 1751 1752 ENDIF ! IF ( topography /= 'flat' ) 1753 1754 END SUBROUTINE pmci_init_loglaw_correction 1675 1755 1676 1756 … … 1901 1981 1902 1982 1903 !-- TO_DO: indentation and wrap long lines from here on to the end of the file 1904 SUBROUTINE pmci_init_anterp_tophat 1905 ! 1906 !-- Precomputation of the client-array indices for 1907 !-- corresponding coarse-grid array index and the 1908 !-- Under-relaxation coefficients to be used by anterp_tophat. 1909 1910 IMPLICIT NONE 1911 1912 INTEGER(iwp) :: i !: Fine-grid index 1913 INTEGER(iwp) :: ii !: Coarse-grid index 1914 INTEGER(iwp) :: istart !: 1915 INTEGER(iwp) :: j !: Fine-grid index 1916 INTEGER(iwp) :: jj !: Coarse-grid index 1917 INTEGER(iwp) :: jstart !: 1918 INTEGER(iwp) :: k !: Fine-grid index 1919 INTEGER(iwp) :: kk !: Coarse-grid index 1920 INTEGER(iwp) :: kstart !: 1921 REAL(wp) :: xi !: 1922 REAL(wp) :: eta !: 1923 REAL(wp) :: zeta !: 1983 1984 1985 SUBROUTINE pmci_init_anterp_tophat 1986 ! 1987 !-- Precomputation of the client-array indices for 1988 !-- corresponding coarse-grid array index and the 1989 !-- Under-relaxation coefficients to be used by anterp_tophat. 1990 1991 IMPLICIT NONE 1992 1993 INTEGER(iwp) :: i !: Fine-grid index 1994 INTEGER(iwp) :: ii !: Coarse-grid index 1995 INTEGER(iwp) :: istart !: 1996 INTEGER(iwp) :: j !: Fine-grid index 1997 INTEGER(iwp) :: jj !: Coarse-grid index 1998 INTEGER(iwp) :: jstart !: 1999 INTEGER(iwp) :: k !: Fine-grid index 2000 INTEGER(iwp) :: kk !: Coarse-grid index 2001 INTEGER(iwp) :: kstart !: 2002 REAL(wp) :: xi !: 2003 REAL(wp) :: eta !: 2004 REAL(wp) :: zeta !: 1924 2005 1925 1926 ! 1927 !-- Default values: 1928 IF ( anterp_relax_length_l < 0.0_wp ) THEN 1929 anterp_relax_length_l = 0.1_wp * ( nx + 1 ) * dx 1930 ENDIF 1931 IF ( anterp_relax_length_r < 0.0_wp ) THEN 1932 anterp_relax_length_r = 0.1_wp * ( nx + 1 ) * dx 1933 ENDIF 1934 IF ( anterp_relax_length_s < 0.0_wp ) THEN 1935 anterp_relax_length_s = 0.1_wp * ( ny + 1 ) * dy 1936 ENDIF 1937 IF ( anterp_relax_length_n < 0.0_wp ) THEN 1938 anterp_relax_length_n = 0.1_wp * ( ny + 1 ) * dy 1939 ENDIF 1940 IF ( anterp_relax_length_t < 0.0_wp ) THEN 1941 anterp_relax_length_t = 0.1_wp * zu(nzt) 1942 ENDIF 1943 1944 ! 1945 !-- First determine kctu and kctw that are the coarse-grid upper bounds for index k. 1946 kk = 0 1947 DO WHILE ( cg%zu(kk) < zu(nzt) ) 1948 kk = kk + 1 1949 ENDDO 1950 kctu = kk - 1 1951 1952 kk = 0 1953 DO WHILE ( cg%zw(kk) < zw(nzt-1) ) 1954 kk = kk + 1 1955 ENDDO 1956 kctw = kk - 1 1957 1958 ALLOCATE( iflu(icl:icr) ) 1959 ALLOCATE( iflo(icl:icr) ) 1960 ALLOCATE( ifuu(icl:icr) ) 1961 ALLOCATE( ifuo(icl:icr) ) 1962 ALLOCATE( jflv(jcs:jcn) ) 1963 ALLOCATE( jflo(jcs:jcn) ) 1964 ALLOCATE( jfuv(jcs:jcn) ) 1965 ALLOCATE( jfuo(jcs:jcn) ) 1966 ALLOCATE( kflw(0:kctw) ) 1967 ALLOCATE( kflo(0:kctu) ) 1968 ALLOCATE( kfuw(0:kctw) ) 1969 ALLOCATE( kfuo(0:kctu) ) 1970 1971 ! 1972 !-- i-indices of u for each l-index value. 1973 istart = nxlg 1974 DO ii = icl, icr 1975 i = istart 1976 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx ) .AND. ( i < nxrg ) ) 1977 i = i + 1 1978 ENDDO 1979 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 1980 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) + 0.5_wp * cg%dx ) .AND. ( i < nxrg ) ) 1981 i = i + 1 1982 ENDDO 1983 ifuu(ii) = MIN( MAX( i, nxlg ), nxrg ) 1984 istart = iflu(ii) 1985 ENDDO 1986 1987 ! 1988 !-- i-indices of others for each l-index value. 1989 istart = nxlg 1990 DO ii = icl, icr 1991 i = istart 1992 DO WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) ) .AND. ( i < nxrg ) ) 1993 i = i + 1 1994 ENDDO 1995 iflo(ii) = MIN( MAX( i, nxlg ), nxrg ) 1996 DO WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) + cg%dx ) .AND. ( i < nxrg ) ) 1997 i = i + 1 1998 ENDDO 1999 ifuo(ii) = MIN(MAX(i,nxlg),nxrg) 2000 istart = iflo(ii) 2001 ENDDO 2002 2003 ! 2004 !-- j-indices of v for each m-index value. 2005 jstart = nysg 2006 DO jj = jcs, jcn 2007 j = jstart 2008 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy ) .AND. ( j < nyng ) ) 2009 j = j + 1 2010 ENDDO 2011 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 2012 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) + 0.5_wp * cg%dy ) .AND. ( j < nyng ) ) 2013 j = j + 1 2014 ENDDO 2015 jfuv(jj) = MIN( MAX( j, nysg ), nyng ) 2016 jstart = jflv(jj) 2017 ENDDO 2018 2019 ! 2020 !-- j-indices of others for each m-index value. 2021 jstart = nysg 2022 DO jj = jcs, jcn 2023 j = jstart 2024 DO WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) ) .AND. ( j < nyng ) ) 2025 j = j + 1 2026 ENDDO 2027 jflo(jj) = MIN( MAX( j, nysg ), nyng ) 2028 DO WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) + cg%dy ) .AND. ( j < nyng ) ) 2029 j = j + 1 2030 ENDDO 2031 jfuo(jj) = MIN( MAX( j, nysg ), nyng ) 2032 jstart = jflv(jj) 2033 ENDDO 2034 2035 ! 2036 !-- k-indices of w for each n-index value. 2037 kstart = 0 2038 kflw(0) = 0 2039 kfuw(0) = 0 2040 DO kk = 1, kctw 2041 k = kstart 2042 DO WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) ) .AND. ( k < nzt ) ) 2043 k = k + 1 2044 ENDDO 2045 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2046 DO WHILE ( ( zw(k) < cg%zw(kk) + 0.5_wp * cg%dzw(kk+1) ) .AND. ( k < nzt ) ) 2047 k = k + 1 2048 ENDDO 2049 kfuw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2050 kstart = kflw(kk) 2051 ENDDO 2052 2053 ! 2054 !-- k-indices of others for each n-index value. 2055 kstart = 0 2056 kflo(0) = 0 2057 kfuo(0) = 0 2058 DO kk = 1, kctu 2059 k = kstart 2060 DO WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) ) .AND. ( k < nzt ) ) 2061 k = k + 1 2062 ENDDO 2063 kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2064 DO WHILE ( ( zu(k) < cg%zu(kk) + 0.5_wp * cg%dzu(kk+1) ) .AND. ( k < nzt ) ) 2065 k = k + 1 2066 ENDDO 2067 kfuo(kk) = MIN( MAX( k-1, 1 ), nzt + 1 ) 2068 kstart = kflo(kk) 2069 ENDDO 2006 ! 2007 !-- Default values: 2008 IF ( anterp_relax_length_l < 0.0_wp ) THEN 2009 anterp_relax_length_l = 0.1_wp * ( nx + 1 ) * dx 2010 ENDIF 2011 IF ( anterp_relax_length_r < 0.0_wp ) THEN 2012 anterp_relax_length_r = 0.1_wp * ( nx + 1 ) * dx 2013 ENDIF 2014 IF ( anterp_relax_length_s < 0.0_wp ) THEN 2015 anterp_relax_length_s = 0.1_wp * ( ny + 1 ) * dy 2016 ENDIF 2017 IF ( anterp_relax_length_n < 0.0_wp ) THEN 2018 anterp_relax_length_n = 0.1_wp * ( ny + 1 ) * dy 2019 ENDIF 2020 IF ( anterp_relax_length_t < 0.0_wp ) THEN 2021 anterp_relax_length_t = 0.1_wp * zu(nzt) 2022 ENDIF 2023 2024 ! 2025 !-- First determine kctu and kctw that are the coarse-grid upper bounds for 2026 !-- index k 2027 kk = 0 2028 DO WHILE ( cg%zu(kk) < zu(nzt) ) 2029 kk = kk + 1 2030 ENDDO 2031 kctu = kk - 1 2032 2033 kk = 0 2034 DO WHILE ( cg%zw(kk) < zw(nzt-1) ) 2035 kk = kk + 1 2036 ENDDO 2037 kctw = kk - 1 2038 2039 ALLOCATE( iflu(icl:icr) ) 2040 ALLOCATE( iflo(icl:icr) ) 2041 ALLOCATE( ifuu(icl:icr) ) 2042 ALLOCATE( ifuo(icl:icr) ) 2043 ALLOCATE( jflv(jcs:jcn) ) 2044 ALLOCATE( jflo(jcs:jcn) ) 2045 ALLOCATE( jfuv(jcs:jcn) ) 2046 ALLOCATE( jfuo(jcs:jcn) ) 2047 ALLOCATE( kflw(0:kctw) ) 2048 ALLOCATE( kflo(0:kctu) ) 2049 ALLOCATE( kfuw(0:kctw) ) 2050 ALLOCATE( kfuo(0:kctu) ) 2051 2052 ! 2053 !-- i-indices of u for each ii-index value 2054 istart = nxlg 2055 DO ii = icl, icr 2056 i = istart 2057 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx ) .AND. & 2058 ( i < nxrg ) ) 2059 i = i + 1 2060 ENDDO 2061 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 2062 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) + 0.5_wp * cg%dx ) .AND. & 2063 ( i < nxrg ) ) 2064 i = i + 1 2065 ENDDO 2066 ifuu(ii) = MIN( MAX( i, nxlg ), nxrg ) 2067 istart = iflu(ii) 2068 ENDDO 2069 2070 ! 2071 !-- i-indices of others for each ii-index value 2072 istart = nxlg 2073 DO ii = icl, icr 2074 i = istart 2075 DO WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) ) .AND. & 2076 ( i < nxrg ) ) 2077 i = i + 1 2078 ENDDO 2079 iflo(ii) = MIN( MAX( i, nxlg ), nxrg ) 2080 DO WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) + cg%dx ) & 2081 .AND. ( i < nxrg ) ) 2082 i = i + 1 2083 ENDDO 2084 ifuo(ii) = MIN(MAX(i,nxlg),nxrg) 2085 istart = iflo(ii) 2086 ENDDO 2087 2088 ! 2089 !-- j-indices of v for each jj-index value 2090 jstart = nysg 2091 DO jj = jcs, jcn 2092 j = jstart 2093 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy ) .AND. & 2094 ( j < nyng ) ) 2095 j = j + 1 2096 ENDDO 2097 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 2098 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) + 0.5_wp * cg%dy ) .AND. & 2099 ( j < nyng ) ) 2100 j = j + 1 2101 ENDDO 2102 jfuv(jj) = MIN( MAX( j, nysg ), nyng ) 2103 jstart = jflv(jj) 2104 ENDDO 2105 2106 ! 2107 !-- j-indices of others for each jj-index value 2108 jstart = nysg 2109 DO jj = jcs, jcn 2110 j = jstart 2111 DO WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) ) .AND. & 2112 ( j < nyng ) ) 2113 j = j + 1 2114 ENDDO 2115 jflo(jj) = MIN( MAX( j, nysg ), nyng ) 2116 DO WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) + cg%dy ) & 2117 .AND. ( j < nyng ) ) 2118 j = j + 1 2119 ENDDO 2120 jfuo(jj) = MIN( MAX( j, nysg ), nyng ) 2121 jstart = jflv(jj) 2122 ENDDO 2123 2124 ! 2125 !-- k-indices of w for each kk-index value 2126 kstart = 0 2127 kflw(0) = 0 2128 kfuw(0) = 0 2129 DO kk = 1, kctw 2130 k = kstart 2131 DO WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) ) .AND. & 2132 ( k < nzt ) ) 2133 k = k + 1 2134 ENDDO 2135 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2136 DO WHILE ( ( zw(k) < cg%zw(kk) + 0.5_wp * cg%dzw(kk+1) ) .AND. & 2137 ( k < nzt ) ) 2138 k = k + 1 2139 ENDDO 2140 kfuw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2141 kstart = kflw(kk) 2142 ENDDO 2143 2144 ! 2145 !-- k-indices of others for each kk-index value 2146 kstart = 0 2147 kflo(0) = 0 2148 kfuo(0) = 0 2149 DO kk = 1, kctu 2150 k = kstart 2151 DO WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) ) .AND. & 2152 ( k < nzt ) ) 2153 k = k + 1 2154 ENDDO 2155 kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 2156 DO WHILE ( ( zu(k) < cg%zu(kk) + 0.5_wp * cg%dzu(kk+1) ) .AND. & 2157 ( k < nzt ) ) 2158 k = k + 1 2159 ENDDO 2160 kfuo(kk) = MIN( MAX( k-1, 1 ), nzt + 1 ) 2161 kstart = kflo(kk) 2162 ENDDO 2070 2163 2071 2164 ! 2072 !-- Spatial under-relaxation coefficients 2073 ALLOCATE( frax(icl:icr) ) 2074 2075 DO ii = icl, icr 2076 IF ( nest_bound_l ) THEN 2077 xi = ( MAX( 0.0_wp, ( cg%coord_x(ii) - lower_left_coord_x ) ) / anterp_relax_length_l )**4 2078 ELSEIF ( nest_bound_r ) THEN 2079 xi = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx - cg%coord_x(ii) ) ) / anterp_relax_length_r )**4 2080 ELSE 2081 xi = 999999.9_wp 2082 ENDIF 2083 frax(ii) = xi / ( 1.0_wp + xi ) 2084 ENDDO 2085 2086 ALLOCATE( fray(jcs:jcn) ) 2087 2088 DO jj = jcs, jcn 2089 IF ( nest_bound_s ) THEN 2090 eta = ( MAX( 0.0_wp, ( cg%coord_y(jj) - lower_left_coord_y ) ) / anterp_relax_length_s )**4 2091 ELSEIF ( nest_bound_n ) THEN 2092 eta = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy - cg%coord_y(jj)) ) / anterp_relax_length_n )**4 2093 ELSE 2094 eta = 999999.9_wp 2095 ENDIF 2096 fray(jj) = eta / ( 1.0_wp + eta ) 2097 ENDDO 2165 !-- Spatial under-relaxation coefficients 2166 ALLOCATE( frax(icl:icr) ) 2167 2168 DO ii = icl, icr 2169 IF ( nest_bound_l ) THEN 2170 xi = ( MAX( 0.0_wp, ( cg%coord_x(ii) - lower_left_coord_x ) ) / & 2171 anterp_relax_length_l )**4 2172 ELSEIF ( nest_bound_r ) THEN 2173 xi = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx - & 2174 cg%coord_x(ii) ) ) / & 2175 anterp_relax_length_r )**4 2176 ELSE 2177 xi = 999999.9_wp 2178 ENDIF 2179 frax(ii) = xi / ( 1.0_wp + xi ) 2180 ENDDO 2181 2182 ALLOCATE( fray(jcs:jcn) ) 2183 2184 DO jj = jcs, jcn 2185 IF ( nest_bound_s ) THEN 2186 eta = ( MAX( 0.0_wp, ( cg%coord_y(jj) - lower_left_coord_y ) ) / & 2187 anterp_relax_length_s )**4 2188 ELSEIF ( nest_bound_n ) THEN 2189 eta = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy - & 2190 cg%coord_y(jj)) ) / & 2191 anterp_relax_length_n )**4 2192 ELSE 2193 eta = 999999.9_wp 2194 ENDIF 2195 fray(jj) = eta / ( 1.0_wp + eta ) 2196 ENDDO 2098 2197 2099 ALLOCATE( fraz(0:kctu) ) 2100 DO kk = 0, kctu 2101 zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4 2102 fraz(kk) = zeta / ( 1.0_wp + zeta ) 2103 ENDDO 2104 2105 END SUBROUTINE pmci_init_anterp_tophat 2106 2107 2108 2109 SUBROUTINE pmci_init_tkefactor 2110 2111 ! 2112 !-- Computes the scaling factor for the SGS TKE from coarse grid to be used 2113 !-- as BC for the fine grid. Based on the Kolmogorov energy spectrum 2114 !-- for the inertial subrange and assumption of sharp cut-off of the resolved 2115 !-- energy spectrum. Near the surface, the reduction of TKE is made 2116 !-- smaller than further away from the surface. 2117 ! 2118 ! Antti Hellsten 4.3.2015 2119 ! 2120 !-- Extended for non-flat topography and variable dz. 2121 ! 2122 ! Antti Hellsten 26.3.2015 2123 ! 2124 !-- The current near-wall adaptation can be replaced by a new one which 2125 !-- uses a step function [0,1] based on the logc-arrays. AH 30.12.2015 2126 IMPLICIT NONE 2127 REAL(wp), PARAMETER :: cfw = 0.2_wp !: 2128 REAL(wp), PARAMETER :: c_tkef = 0.6_wp !: 2129 REAL(wp) :: fw !: 2130 REAL(wp), PARAMETER :: fw0 = 0.9_wp !: 2131 REAL(wp) :: glsf !: 2132 REAL(wp) :: glsc !: 2133 REAL(wp) :: height !: 2134 REAL(wp), PARAMETER :: p13 = 1.0_wp/3.0_wp !: 2135 REAL(wp), PARAMETER :: p23 = 2.0_wp/3.0_wp !: 2136 INTEGER(iwp) :: k !: 2137 INTEGER(iwp) :: kc !: 2198 ALLOCATE( fraz(0:kctu) ) 2199 DO kk = 0, kctu 2200 zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4 2201 fraz(kk) = zeta / ( 1.0_wp + zeta ) 2202 ENDDO 2203 2204 END SUBROUTINE pmci_init_anterp_tophat 2205 2206 2207 2208 SUBROUTINE pmci_init_tkefactor 2209 2210 ! 2211 !-- Computes the scaling factor for the SGS TKE from coarse grid to be used 2212 !-- as BC for the fine grid. Based on the Kolmogorov energy spectrum 2213 !-- for the inertial subrange and assumption of sharp cut-off of the resolved 2214 !-- energy spectrum. Near the surface, the reduction of TKE is made 2215 !-- smaller than further away from the surface. 2216 2217 IMPLICIT NONE 2218 REAL(wp), PARAMETER :: cfw = 0.2_wp !: 2219 REAL(wp), PARAMETER :: c_tkef = 0.6_wp !: 2220 REAL(wp) :: fw !: 2221 REAL(wp), PARAMETER :: fw0 = 0.9_wp !: 2222 REAL(wp) :: glsf !: 2223 REAL(wp) :: glsc !: 2224 REAL(wp) :: height !: 2225 REAL(wp), PARAMETER :: p13 = 1.0_wp/3.0_wp !: 2226 REAL(wp), PARAMETER :: p23 = 2.0_wp/3.0_wp !: 2227 INTEGER(iwp) :: k !: 2228 INTEGER(iwp) :: kc !: 2138 2229 2139 2230 2140 IF ( nest_bound_l ) THEN 2141 ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) ) 2142 tkefactor_l = 0.0_wp 2143 i = nxl - 1 2144 DO j = nysg, nyng 2145 DO k = nzb_s_inner(j,i) + 1, nzt 2146 kc = kco(k+1) 2147 glsf = ( dx * dy * dzu(k) )**p13 2148 glsc = ( cg%dx * cg%dy *cg%dzu(kc) )**p13 2149 height = zu(k) - zu(nzb_s_inner(j,i)) 2150 fw = EXP( -cfw * height / glsf ) 2151 tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 ) 2152 ENDDO 2153 tkefactor_l(nzb_s_inner(j,i),j) = c_tkef * fw0 2154 ENDDO 2155 ENDIF 2156 2157 IF ( nest_bound_r ) THEN 2158 ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) ) 2159 tkefactor_r = 0.0_wp 2160 i = nxr + 1 2161 DO j = nysg, nyng 2162 DO k = nzb_s_inner(j,i) + 1, nzt 2163 kc = kco(k+1) 2164 glsf = ( dx * dy * dzu(k) )**p13 2165 glsc = ( cg%dx * cg%dy * cg%dzu(kc) )**p13 2166 height = zu(k) - zu(nzb_s_inner(j,i)) 2167 fw = EXP( -cfw * height / glsf ) 2168 tkefactor_r(k,j) = c_tkef * (fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 ) 2169 ENDDO 2170 tkefactor_r(nzb_s_inner(j,i),j) = c_tkef * fw0 2171 ENDDO 2172 ENDIF 2231 IF ( nest_bound_l ) THEN 2232 ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) ) 2233 tkefactor_l = 0.0_wp 2234 i = nxl - 1 2235 DO j = nysg, nyng 2236 DO k = nzb_s_inner(j,i) + 1, nzt 2237 kc = kco(k+1) 2238 glsf = ( dx * dy * dzu(k) )**p13 2239 glsc = ( cg%dx * cg%dy *cg%dzu(kc) )**p13 2240 height = zu(k) - zu(nzb_s_inner(j,i)) 2241 fw = EXP( -cfw * height / glsf ) 2242 tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * & 2243 ( glsf / glsc )**p23 ) 2244 ENDDO 2245 tkefactor_l(nzb_s_inner(j,i),j) = c_tkef * fw0 2246 ENDDO 2247 ENDIF 2248 2249 IF ( nest_bound_r ) THEN 2250 ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) ) 2251 tkefactor_r = 0.0_wp 2252 i = nxr + 1 2253 DO j = nysg, nyng 2254 DO k = nzb_s_inner(j,i) + 1, nzt 2255 kc = kco(k+1) 2256 glsf = ( dx * dy * dzu(k) )**p13 2257 glsc = ( cg%dx * cg%dy * cg%dzu(kc) )**p13 2258 height = zu(k) - zu(nzb_s_inner(j,i)) 2259 fw = EXP( -cfw * height / glsf ) 2260 tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * & 2261 ( glsf / glsc )**p23 ) 2262 ENDDO 2263 tkefactor_r(nzb_s_inner(j,i),j) = c_tkef * fw0 2264 ENDDO 2265 ENDIF 2173 2266 2174 2267 IF ( nest_bound_s ) THEN 2175 ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) ) 2176 tkefactor_s = 0.0_wp 2177 j = nys - 1 2178 DO i = nxlg, nxrg 2179 DO k = nzb_s_inner(j,i) + 1, nzt 2180 kc = kco(k+1) 2181 glsf = ( dx * dy * dzu(k) )**p13 2182 glsc = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13 2183 height = zu(k) - zu(nzb_s_inner(j,i)) 2184 fw = EXP( -cfw*height / glsf ) 2185 tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 ) 2186 ENDDO 2187 tkefactor_s(nzb_s_inner(j,i),i) = c_tkef * fw0 2188 ENDDO 2189 ENDIF 2190 2191 IF ( nest_bound_n ) THEN 2192 ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) ) 2193 tkefactor_n = 0.0_wp 2194 j = nyn + 1 2195 DO i = nxlg, nxrg 2196 DO k = nzb_s_inner(j,i)+1, nzt 2197 kc = kco(k+1) 2198 glsf = ( dx * dy * dzu(k) )**p13 2199 glsc = ( cg%dx * cg%dy * cg%dzu(kc) )**p13 2200 height = zu(k) - zu(nzb_s_inner(j,i)) 2201 fw = EXP( -cfw * height / glsf ) 2202 tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 ) 2203 ENDDO 2204 tkefactor_n(nzb_s_inner(j,i),i) = c_tkef * fw0 2205 ENDDO 2206 ENDIF 2207 2208 ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) ) 2209 k = nzt 2210 DO i = nxlg, nxrg 2211 DO j = nysg, nyng 2212 kc = kco(k+1) 2213 glsf = ( dx * dy * dzu(k) )**p13 2214 glsc = ( cg%dx * cg%dy * cg%dzu(kc) )**p13 2215 height = zu(k) - zu(nzb_s_inner(j,i)) 2216 fw = EXP( -cfw * height / glsf ) 2217 tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 ) 2218 ENDDO 2219 ENDDO 2268 ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) ) 2269 tkefactor_s = 0.0_wp 2270 j = nys - 1 2271 DO i = nxlg, nxrg 2272 DO k = nzb_s_inner(j,i) + 1, nzt 2273 kc = kco(k+1) 2274 glsf = ( dx * dy * dzu(k) )**p13 2275 glsc = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13 2276 height = zu(k) - zu(nzb_s_inner(j,i)) 2277 fw = EXP( -cfw*height / glsf ) 2278 tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * & 2279 ( glsf / glsc )**p23 ) 2280 ENDDO 2281 tkefactor_s(nzb_s_inner(j,i),i) = c_tkef * fw0 2282 ENDDO 2283 ENDIF 2284 2285 IF ( nest_bound_n ) THEN 2286 ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) ) 2287 tkefactor_n = 0.0_wp 2288 j = nyn + 1 2289 DO i = nxlg, nxrg 2290 DO k = nzb_s_inner(j,i)+1, nzt 2291 kc = kco(k+1) 2292 glsf = ( dx * dy * dzu(k) )**p13 2293 glsc = ( cg%dx * cg%dy * cg%dzu(kc) )**p13 2294 height = zu(k) - zu(nzb_s_inner(j,i)) 2295 fw = EXP( -cfw * height / glsf ) 2296 tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * & 2297 ( glsf / glsc )**p23 ) 2298 ENDDO 2299 tkefactor_n(nzb_s_inner(j,i),i) = c_tkef * fw0 2300 ENDDO 2301 ENDIF 2302 2303 ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) ) 2304 k = nzt 2305 DO i = nxlg, nxrg 2306 DO j = nysg, nyng 2307 kc = kco(k+1) 2308 glsf = ( dx * dy * dzu(k) )**p13 2309 glsc = ( cg%dx * cg%dy * cg%dzu(kc) )**p13 2310 height = zu(k) - zu(nzb_s_inner(j,i)) 2311 fw = EXP( -cfw * height / glsf ) 2312 tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * & 2313 ( glsf / glsc )**p23 ) 2314 ENDDO 2315 ENDDO 2220 2316 2221 2317 END SUBROUTINE pmci_init_tkefactor … … 2366 2462 2367 2463 ! 2368 !-- List of array names, which can be coupled. 2369 !-- TO_DO: Antti: what is the meaning of the next line? 2370 !-- AH: Note that the k-range of the *c arrays is changed from 1:nz to 0:nz+1. 2464 !-- List of array names, which can be coupled 2371 2465 IF ( TRIM( name ) == "u" ) THEN 2372 2466 IF ( .NOT. ALLOCATED( uc ) ) ALLOCATE( uc(0:nzc+1, js:je, is:ie) ) … … 2406 2500 ELSE 2407 2501 ! 2408 !-- Avoid others to continue2502 !-- Prevent others from continuing 2409 2503 CALL MPI_BARRIER( comm2d, ierr ) 2410 2504 ENDIF … … 2650 2744 2651 2745 2746 SUBROUTINE pmci_check_setting_mismatches 2747 ! 2748 !-- Check for mismatches between settings of master and client variables 2749 !-- (e.g., all clients have to follow the end_time settings of the root model). 2750 !-- The root model overwrites variables in the other models, so these variables 2751 !-- only need to be set once in file PARIN. 2752 2753 #if defined( __parallel ) 2754 2755 USE control_parameters, & 2756 ONLY: dt_restart, end_time, message_string, restart_time, time_restart 2757 2758 IMPLICIT NONE 2759 2760 INTEGER :: ierr 2761 2762 REAL(wp) :: dt_restart_root 2763 REAL(wp) :: end_time_root 2764 REAL(wp) :: restart_time_root 2765 REAL(wp) :: time_restart_root 2766 2767 ! 2768 !-- Check the time to be simulated. 2769 !-- Here, and in the following, the root process communicates the respective 2770 !-- variable to all others, and its value will then be compared with the local 2771 !-- values. 2772 IF ( pmc_is_rootmodel() ) end_time_root = end_time 2773 CALL MPI_BCAST( end_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr ) 2774 2775 IF ( .NOT. pmc_is_rootmodel() ) THEN 2776 IF ( end_time /= end_time_root ) THEN 2777 WRITE( message_string, * ) 'mismatch between root model and ', & 2778 'client settings & end_time(root) = ', end_time_root, & 2779 ' & end_time(client) = ', end_time, ' & client value is set', & 2780 ' to root value' 2781 CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, & 2782 0 ) 2783 end_time = end_time_root 2784 ENDIF 2785 ENDIF 2786 2787 ! 2788 !-- Same for restart time 2789 IF ( pmc_is_rootmodel() ) restart_time_root = restart_time 2790 CALL MPI_BCAST( restart_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr ) 2791 2792 IF ( .NOT. pmc_is_rootmodel() ) THEN 2793 IF ( restart_time /= restart_time_root ) THEN 2794 WRITE( message_string, * ) 'mismatch between root model and ', & 2795 'client settings & restart_time(root) = ', restart_time_root, & 2796 ' & restart_time(client) = ', restart_time, ' & client ', & 2797 'value is set to root value' 2798 CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, & 2799 0 ) 2800 restart_time = restart_time_root 2801 ENDIF 2802 ENDIF 2803 2804 ! 2805 !-- Same for dt_restart 2806 IF ( pmc_is_rootmodel() ) dt_restart_root = dt_restart 2807 CALL MPI_BCAST( dt_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr ) 2808 2809 IF ( .NOT. pmc_is_rootmodel() ) THEN 2810 IF ( dt_restart /= dt_restart_root ) THEN 2811 WRITE( message_string, * ) 'mismatch between root model and ', & 2812 'client settings & dt_restart(root) = ', dt_restart_root, & 2813 ' & dt_restart(client) = ', dt_restart, ' & client ', & 2814 'value is set to root value' 2815 CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, & 2816 0 ) 2817 dt_restart = dt_restart_root 2818 ENDIF 2819 ENDIF 2820 2821 ! 2822 !-- Same for time_restart 2823 IF ( pmc_is_rootmodel() ) time_restart_root = time_restart 2824 CALL MPI_BCAST( time_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr ) 2825 2826 IF ( .NOT. pmc_is_rootmodel() ) THEN 2827 IF ( time_restart /= time_restart_root ) THEN 2828 WRITE( message_string, * ) 'mismatch between root model and ', & 2829 'client settings & time_restart(root) = ', time_restart_root, & 2830 ' & time_restart(client) = ', time_restart, ' & client ', & 2831 'value is set to root value' 2832 CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, & 2833 0 ) 2834 time_restart = time_restart_root 2835 ENDIF 2836 ENDIF 2837 2838 #endif 2839 2840 END SUBROUTINE pmci_check_setting_mismatches 2841 2842 2843 2652 2844 SUBROUTINE pmci_ensure_nest_mass_conservation 2653 2845 … … 2787 2979 REAL(wp), DIMENSION(1) :: dtl !: 2788 2980 2981 2982 CALL cpu_log( log_point_s(70), 'pmc sync', 'start' ) 2983 2789 2984 ! 2790 2985 !-- First find the smallest native time step of all the clients of the current … … 2813 3008 ENDDO 2814 3009 3010 CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' ) 3011 2815 3012 #endif 2816 3013 END SUBROUTINE pmci_server_synchronize … … 2836 3033 dtl(1) = dt_3d 2837 3034 IF ( cpl_id > 1 ) THEN 3035 3036 CALL cpu_log( log_point_s(70), 'pmc sync', 'start' ) 3037 2838 3038 IF ( myid==0 ) THEN 2839 3039 CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr ) … … 2845 3045 !-- Broadcast the unified time step to all server processes 2846 3046 CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr ) 3047 3048 CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' ) 3049 2847 3050 ENDIF 2848 3051 … … 2853 3056 2854 3057 SUBROUTINE pmci_set_swaplevel( swaplevel ) 3058 ! 3059 !-- After each Runge-Kutta sub-timestep, alternately set buffer one or buffer 3060 !-- two active 2855 3061 2856 3062 IMPLICIT NONE … … 2862 3068 INTEGER(iwp) :: m !: 2863 3069 2864 !2865 !-- After each timestep, alternately set buffer one or buffer two active2866 3070 DO m = 1, SIZE( pmc_server_for_client )-1 2867 3071 client_id = pmc_server_for_client(m) … … 2870 3074 2871 3075 END SUBROUTINE pmci_set_swaplevel 3076 3077 3078 3079 SUBROUTINE pmci_datatrans( local_nesting_mode ) 3080 ! 3081 !-- Althoug nesting_mode is a variable of this model, pass it as an argument to 3082 !-- allow for example to force one-way initialization phase 3083 3084 IMPLICIT NONE 3085 3086 INTEGER(iwp) :: ierr !: 3087 INTEGER(iwp) :: istat !: 3088 3089 CHARACTER(LEN=*),INTENT(IN) :: local_nesting_mode 3090 3091 IF ( local_nesting_mode == 'one-way' ) THEN 3092 3093 CALL pmci_client_datatrans( server_to_client ) 3094 CALL pmci_server_datatrans( server_to_client ) 3095 3096 ELSE 3097 3098 IF( nesting_datatransfer_mode == 'cascade' ) THEN 3099 3100 CALL pmci_client_datatrans( server_to_client ) 3101 CALL pmci_server_datatrans( server_to_client ) 3102 3103 CALL pmci_server_datatrans( client_to_server ) 3104 CALL pmci_client_datatrans( client_to_server ) 3105 3106 ELSEIF( nesting_datatransfer_mode == 'overlap') THEN 3107 3108 CALL pmci_server_datatrans( server_to_client ) 3109 CALL pmci_client_datatrans( server_to_client ) 3110 3111 CALL pmci_client_datatrans( client_to_server ) 3112 CALL pmci_server_datatrans( client_to_server ) 3113 3114 ELSEIF( TRIM( nesting_datatransfer_mode ) == 'mixed' ) THEN 3115 3116 CALL pmci_server_datatrans( server_to_client ) 3117 CALL pmci_client_datatrans( server_to_client ) 3118 3119 CALL pmci_server_datatrans( client_to_server ) 3120 CALL pmci_client_datatrans( client_to_server ) 3121 3122 ENDIF 3123 3124 ENDIF 3125 3126 END SUBROUTINE pmci_datatrans 3127 2872 3128 2873 3129 … … 2890 3146 REAL(wp), DIMENSION(1) :: dtl !: 2891 3147 2892 !2893 !-- First find the smallest native time step of all the clients of the current2894 !-- server.2895 dtl(1) = 999999.9_wp2896 DO m = 1, SIZE( pmc_server_for_client )-12897 client_id = pmc_server_for_client(m)2898 IF ( myid==0 ) THEN2899 CALL pmc_recv_from_client( client_id, dtc, SIZE( dtc ), 0, 101, ierr )2900 dtl(1) = MIN( dtl(1), dtc(1) )2901 dt_3d = dtl(1)2902 ENDIF2903 ENDDO2904 2905 !2906 !-- Broadcast the unified time step to all server processes2907 CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )2908 3148 2909 3149 DO m = 1, SIZE( PMC_Server_for_Client )-1 2910 3150 client_id = PMC_Server_for_Client(m) 2911 CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )2912 2913 !2914 !-- Send the new time step to all the clients of the current server2915 IF ( myid == 0 ) THEN2916 CALL pmc_send_to_client( client_id, dtl, SIZE( dtl ), 0, 102, ierr )2917 ENDIF2918 CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )2919 3151 2920 3152 IF ( direction == server_to_client ) THEN 2921 3153 CALL cpu_log( log_point_s(71), 'pmc server send', 'start' ) 2922 CALL pmc_s_fillbuffer( client_id , waittime=waittime)3154 CALL pmc_s_fillbuffer( client_id ) 2923 3155 CALL cpu_log( log_point_s(71), 'pmc server send', 'stop' ) 2924 3156 ELSE … … 2936 3168 ! 2937 3169 !-- Inside buildings/topography reset velocities and TKE back to zero. 2938 !-- TO_DO: at least temperature should be included here immediately2939 3170 !-- Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at 2940 3171 !-- present, maybe revise later. 2941 3172 DO i = nxlg, nxrg 2942 3173 DO j = nysg, nyng 2943 u(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp 2944 v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp 2945 w(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 2946 e(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp 3174 u(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp 3175 v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp 3176 w(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 3177 e(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp 3178 pt(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp 3179 IF ( humidity .OR. passive_scalar ) THEN 3180 q(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp 3181 ENDIF 2947 3182 ENDDO 2948 3183 ENDDO … … 2969 3204 INTEGER(iwp) :: jcn !: 2970 3205 2971 REAL(wp) :: waittime !:2972 2973 3206 REAL(wp), DIMENSION(1) :: dtl !: 2974 3207 REAL(wp), DIMENSION(1) :: dts !: … … 2977 3210 dtl = dt_3d 2978 3211 IF ( cpl_id > 1 ) THEN 2979 CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )2980 IF ( myid==0 ) THEN2981 CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )2982 CALL pmc_recv_from_server( dts, SIZE( dts ), 0, 102, ierr )2983 dt_3d = dts(1)2984 ENDIF2985 2986 !2987 !-- Broadcast the unified time step to all server processes.2988 CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )2989 CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )2990 2991 3212 ! 2992 3213 !-- Client domain boundaries in the server indice space. … … 2999 3220 3000 3221 CALL cpu_log( log_point_s(73), 'pmc client recv', 'start' ) 3001 CALL pmc_c_getbuffer( WaitTime = WaitTime)3222 CALL pmc_c_getbuffer( ) 3002 3223 CALL cpu_log( log_point_s(73), 'pmc client recv', 'stop' ) 3003 3224 3225 CALL cpu_log( log_point_s(75), 'pmc interpolation', 'start' ) 3004 3226 CALL pmci_interpolation 3227 CALL cpu_log( log_point_s(75), 'pmc interpolation', 'stop' ) 3005 3228 3006 3229 ELSE 3007 3230 ! 3008 !-- direction == server_to_client 3231 !-- direction == client_to_server 3232 CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'start' ) 3009 3233 CALL pmci_anterpolation 3234 CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'stop' ) 3010 3235 3011 3236 CALL cpu_log( log_point_s(74), 'pmc client send', 'start' ) 3012 CALL pmc_c_putbuffer( WaitTime = WaitTime)3237 CALL pmc_c_putbuffer( ) 3013 3238 CALL cpu_log( log_point_s(74), 'pmc client send', 'stop' ) 3014 3239 … … 3234 3459 3235 3460 CALL pmci_anterp_tophat( u, uc, kctu, iflu, ifuu, jflo, jfuo, kflo, & 3236 kfuo, nzb_u_inner,'u' )3461 kfuo, 'u' ) 3237 3462 CALL pmci_anterp_tophat( v, vc, kctu, iflo, ifuo, jflv, jfuv, kflo, & 3238 kfuo, nzb_v_inner,'v' )3463 kfuo, 'v' ) 3239 3464 CALL pmci_anterp_tophat( w, wc, kctw, iflo, ifuo, jflo, jfuo, kflw, & 3240 kfuw, nzb_w_inner,'w' )3465 kfuw, 'w' ) 3241 3466 CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo, & 3242 kfuo, nzb_s_inner,'s' )3467 kfuo, 's' ) 3243 3468 IF ( humidity .OR. passive_scalar ) THEN 3244 3469 CALL pmci_anterp_tophat( q, qc, kctu, iflo, ifuo, jflo, jfuo, kflo, & 3245 kfuo, nzb_s_inner,'s' )3470 kfuo, 's' ) 3246 3471 ENDIF 3247 3472 … … 3256 3481 !-- Interpolation of ghost-node values used as the client-domain boundary 3257 3482 !-- conditions. This subroutine handles the left and right boundaries. It is 3258 !-- based on trilinear interpolation. Constant dz is still assumed. 3259 !-- TO_DO: constant dz is an important restriction and should be checked 3260 !-- somewhere in order to let users not run into this trap 3483 !-- based on trilinear interpolation. 3484 3261 3485 IMPLICIT NONE 3262 3486 3263 !-- TO_DO: wrap long lines in this and the remaining interp_tril routines 3264 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !: 3265 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc !: 3266 REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2,0:ncorr-1), INTENT(IN) :: logc_ratio !: 3267 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !: 3268 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !: 3269 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !: 3270 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !: 3271 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !: 3272 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !: 3487 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 3488 INTENT(INOUT) :: f !: 3489 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), & 3490 INTENT(IN) :: fc !: 3491 REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2,0:ncorr-1), & 3492 INTENT(IN) :: logc_ratio !: 3493 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !: 3494 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !: 3495 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !: 3496 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !: 3497 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !: 3498 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !: 3273 3499 3274 INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN) :: ic !: 3275 INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) :: jc !: 3276 INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) :: kb !: 3277 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !: 3278 INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2), INTENT(IN) :: logc !: 3279 INTEGER(iwp) :: nzt_topo_nestbc !: 3500 INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN) :: ic !: 3501 INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) :: jc !: 3502 INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) :: kb !: 3503 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !: 3504 INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2), & 3505 INTENT(IN) :: logc !: 3506 INTEGER(iwp) :: nzt_topo_nestbc !: 3280 3507 3281 3508 CHARACTER(LEN=1),INTENT(IN) :: edge !: … … 3312 3539 3313 3540 ! 3314 !-- Check which edge is to be handled: left or right. Note the assumption that the same PE never 3315 !-- holds both left and right nest boundaries. Should this be changed? 3541 !-- Check which edge is to be handled 3316 3542 IF ( edge == 'l' ) THEN 3317 IF ( var == 'u' ) THEN ! For u, nxl is a ghost node, but not for the other variables. 3543 ! 3544 !-- For u, nxl is a ghost node, but not for the other variables 3545 IF ( var == 'u' ) THEN 3318 3546 i = nxl 3319 3547 ib = nxl - 1 … … 3327 3555 ENDIF 3328 3556 3329 DO j = nys, nyn +13330 DO k = kb(j,i), nzt +13557 DO j = nys, nyn+1 3558 DO k = kb(j,i), nzt+1 3331 3559 l = ic(i) 3332 3560 m = jc(j) … … 3345 3573 !-- Generalized log-law-correction algorithm. 3346 3574 !-- Doubly two-dimensional index arrays logc(:,:,1:2) and log-ratio arrays 3347 !-- logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine pmci_init_loglaw_correction. 3575 !-- logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine 3576 !-- pmci_init_loglaw_correction. 3348 3577 ! 3349 3578 !-- Solid surface below the node 3350 3579 IF ( var == 'u' .OR. var == 'v' ) THEN 3351 3580 DO j = nys, nyn 3352 k = kb(j,i) +13581 k = kb(j,i)+1 3353 3582 IF ( ( logc(k,j,1) /= 0 ) .AND. ( logc(k,j,2) == 0 ) ) THEN 3354 3583 k1 = logc(k,j,1) … … 3362 3591 3363 3592 ! 3364 !-- In case of non-flat topography, also vertical walls and corners need to be treated. 3365 !-- Only single and double wall nodes are corrected. Triple and higher-multiple wall nodes 3366 !-- are not corrected as the log law would not be valid anyway in such locations. 3593 !-- In case of non-flat topography, also vertical walls and corners need to be 3594 !-- treated. Only single and double wall nodes are corrected. Triple and 3595 !-- higher-multiple wall nodes are not corrected as the log law would not be 3596 !-- valid anyway in such locations. 3367 3597 IF ( topography /= 'flat' ) THEN 3368 3598 IF ( var == 'u' .OR. var == 'w' ) THEN … … 3371 3601 !-- Solid surface only on south/north side of the node 3372 3602 DO j = nys, nyn 3373 DO k = kb(j,i) +1, nzt_topo_nestbc3603 DO k = kb(j,i)+1, nzt_topo_nestbc 3374 3604 IF ( ( logc(k,j,2) /= 0 ) .AND. ( logc(k,j,1) == 0 ) ) THEN 3375 3605 3376 3606 ! 3377 !-- Direction of the wall-normal index is carried in as the sign of logc. 3607 !-- Direction of the wall-normal index is carried in as the 3608 !-- sign of logc 3378 3609 jinc = SIGN( 1, logc(k,j,2) ) 3379 3610 j1 = ABS( logc(k,j,2) ) 3380 DO jcorr =0, ncorr -13611 DO jcorr = 0, ncorr-1 3381 3612 jco = j + jinc * jcorr 3382 3613 f(k,jco,i) = logc_ratio(k,j,2,jcorr) * f(k,j1,i) … … 3396 3627 jinc = SIGN( 1, logc(k,j,2) ) 3397 3628 j1 = ABS( logc(k,j,2) ) 3398 DO jcorr = 0, ncorr -13629 DO jcorr = 0, ncorr-1 3399 3630 jco = j + jinc * jcorr 3400 DO kcorr = 0, ncorr -13631 DO kcorr = 0, ncorr-1 3401 3632 kco = k + kcorr 3402 f(kco,jco,i) = 0.5_wp * ( logc_ratio(k,j,1,kcorr) * f(k1,j,i) & 3403 + logc_ratio(k,j,2,jcorr) * f(k,j1,i) ) 3633 f(kco,jco,i) = 0.5_wp * ( logc_ratio(k,j,1,kcorr) * & 3634 f(k1,j,i) & 3635 + logc_ratio(k,j,2,jcorr) * & 3636 f(k,j1,i) ) 3404 3637 ENDDO 3405 3638 ENDDO … … 3420 3653 ENDDO 3421 3654 ELSEIF ( edge == 'r' ) THEN 3422 DO j = nys, nyn +13423 DO k = kb(j,i), nzt +13655 DO j = nys, nyn+1 3656 DO k = kb(j,i), nzt+1 3424 3657 f(k,j,i) = tkefactor_r(k,j) * f(k,j,i) 3425 3658 ENDDO … … 3431 3664 !-- Store the boundary values also into the other redundant ghost node layers 3432 3665 IF ( edge == 'l' ) THEN 3433 DO ibgp = -nbgp, ib3666 DO ibgp = -nbgp, ib 3434 3667 f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i) 3435 3668 ENDDO 3436 ELSE 3437 DO ibgp = ib, nx +nbgp3669 ELSEIF ( edge == 'r' ) THEN 3670 DO ibgp = ib, nx+nbgp 3438 3671 f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i) 3439 3672 ENDDO … … 3444 3677 3445 3678 3446 SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, logc, logc_ratio, & 3447 nzt_topo_nestbc, edge, var ) 3679 SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, & 3680 r2z, kb, logc, logc_ratio, & 3681 nzt_topo_nestbc, edge, var ) 3448 3682 3449 3683 ! … … 3451 3685 !-- conditions. This subroutine handles the south and north boundaries. 3452 3686 !-- This subroutine is based on trilinear interpolation. 3453 !-- Constant dz is still assumed. 3454 ! 3455 !-- Antti Hellsten 22.2.2015. 3456 ! 3457 !-- Rewritten so that all the coefficients and client-array indices are 3458 !-- precomputed in the initialization phase by pmci_init_interp_tril. 3459 ! 3460 !-- Antti Hellsten 3.3.2015. 3461 ! 3462 !-- Constant dz no more assumed. 3463 !-- Antti Hellsten 23.3.2015. 3464 ! 3465 !-- Adapted for non-flat topography. However, the near-wall velocities 3466 !-- are log-corrected only over horifontal surfaces, not yet near vertical 3467 !-- walls. 3468 !-- Antti Hellsten 26.3.2015. 3469 ! 3470 !-- Indexing in the principal direction (j) is changed. Now, the nest-boundary 3471 !-- values are interpolated only into the first ghost-node layers on each later 3472 !-- boundary. These values are then simply copied to the second ghost-node layer. 3473 ! 3474 !-- Antti Hellsten 6.10.2015. 3687 3475 3688 IMPLICIT NONE 3476 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !: 3477 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc !: 3478 REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2,0:ncorr-1), INTENT(IN) :: logc_ratio !: 3479 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !: 3480 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !: 3481 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !: 3482 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !: 3483 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !: 3484 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !: 3689 3690 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 3691 INTENT(INOUT) :: f !: 3692 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), & 3693 INTENT(IN) :: fc !: 3694 REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2,0:ncorr-1), & 3695 INTENT(IN) :: logc_ratio !: 3696 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !: 3697 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !: 3698 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !: 3699 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !: 3700 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !: 3701 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !: 3485 3702 3486 INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN) :: ic !: 3487 INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) :: jc !: 3488 INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) :: kb !: 3489 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !: 3490 INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2), INTENT(IN) :: logc !: 3491 INTEGER(iwp) :: nzt_topo_nestbc !: 3703 INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN) :: ic !: 3704 INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) :: jc !: 3705 INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) :: kb !: 3706 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !: 3707 INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2), & 3708 INTENT(IN) :: logc !: 3709 INTEGER(iwp) :: nzt_topo_nestbc !: 3492 3710 3493 3711 CHARACTER(LEN=1), INTENT(IN) :: edge !: … … 3521 3739 3522 3740 ! 3523 !-- Check which edge is to be handled: south or north. Note the assumption that the same PE never 3524 !-- holds both south and north nest boundaries. Should this be changed? 3741 !-- Check which edge is to be handled: south or north 3525 3742 IF ( edge == 's' ) THEN 3526 IF ( var == 'v' ) THEN ! For v, nys is a ghost node, but not for the other variables. 3743 ! 3744 !-- For v, nys is a ghost node, but not for the other variables 3745 IF ( var == 'v' ) THEN 3527 3746 j = nys 3528 3747 jb = nys - 1 … … 3536 3755 ENDIF 3537 3756 3538 DO i = nxl, nxr +13539 DO k = kb(j,i), nzt +13757 DO i = nxl, nxr+1 3758 DO k = kb(j,i), nzt+1 3540 3759 l = ic(i) 3541 3760 m = jc(j) … … 3554 3773 !-- Generalized log-law-correction algorithm. 3555 3774 !-- Multiply two-dimensional index arrays logc(:,:,1:2) and log-ratio arrays 3556 !-- logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine pmci_init_loglaw_correction. 3775 !-- logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine 3776 !-- pmci_init_loglaw_correction. 3557 3777 ! 3558 3778 !-- Solid surface below the node … … 3562 3782 IF ( ( logc(k,i,1) /= 0 ) .AND. ( logc(k,i,2) == 0 ) ) THEN 3563 3783 k1 = logc(k,i,1) 3564 DO kcorr = 0, ncorr -13784 DO kcorr = 0, ncorr-1 3565 3785 kco = k + kcorr 3566 3786 f(kco,j,i) = logc_ratio(k,i,1,kcorr) * f(k1,j,i) … … 3571 3791 3572 3792 ! 3573 !-- In case of non-flat topography, also vertical walls and corners need to be treated. 3574 !-- Only single and double wall nodes are corrected. 3575 !-- Triple and higher-multiple wall nodes are not corrected as it would be extremely complicated 3576 !-- and the log law would not be valid anyway in such locations. 3793 !-- In case of non-flat topography, also vertical walls and corners need to be 3794 !-- treated. Only single and double wall nodes are corrected. 3795 !-- Triple and higher-multiple wall nodes are not corrected as it would be 3796 !-- extremely complicated and the log law would not be valid anyway in such 3797 !-- locations. 3577 3798 IF ( topography /= 'flat' ) THEN 3578 3799 IF ( var == 'v' .OR. var == 'w' ) THEN … … 3585 3806 3586 3807 ! 3587 !-- Direction of the wall-normal index is carried in as the sign of logc. 3808 !-- Direction of the wall-normal index is carried in as the 3809 !-- sign of logc 3588 3810 iinc = SIGN( 1, logc(k,i,2) ) 3589 3811 i1 = ABS( logc(k,i,2) ) 3590 DO icorr = 0, ncorr -13812 DO icorr = 0, ncorr-1 3591 3813 ico = i + iinc * icorr 3592 3814 f(k,j,ico) = logc_ratio(k,i,2,icorr) * f(k,j,i1) … … 3606 3828 iinc = SIGN( 1, logc(k,i,2) ) 3607 3829 i1 = ABS( logc(k,i,2) ) 3608 DO icorr = 0, ncorr -13830 DO icorr = 0, ncorr-1 3609 3831 ico = i + iinc * icorr 3610 DO kcorr = 0, ncorr -13832 DO kcorr = 0, ncorr-1 3611 3833 kco = k + kcorr 3612 f(kco,i,ico) = 0.5_wp * ( logc_ratio(k,i,1,kcorr) * f(k1,j,i) & 3613 + logc_ratio(k,i,2,icorr) * f(k,j,i1) ) 3834 f(kco,i,ico) = 0.5_wp * ( logc_ratio(k,i,1,kcorr) * & 3835 f(k1,j,i) & 3836 + logc_ratio(k,i,2,icorr) * & 3837 f(k,j,i1) ) 3614 3838 ENDDO 3615 3839 ENDDO … … 3625 3849 IF ( edge == 's' ) THEN 3626 3850 DO i = nxl, nxr + 1 3627 DO k = kb(j,i), nzt +13851 DO k = kb(j,i), nzt+1 3628 3852 f(k,j,i) = tkefactor_s(k,i) * f(k,j,i) 3629 3853 ENDDO … … 3631 3855 ELSEIF ( edge == 'n' ) THEN 3632 3856 DO i = nxl, nxr + 1 3633 DO k = kb(j,i), nzt +13857 DO k = kb(j,i), nzt+1 3634 3858 f(k,j,i) = tkefactor_n(k,i) * f(k,j,i) 3635 3859 ENDDO … … 3641 3865 !-- Store the boundary values also into the other redundant ghost node layers 3642 3866 IF ( edge == 's' ) THEN 3643 DO jbgp = -nbgp, jb3867 DO jbgp = -nbgp, jb 3644 3868 f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg) 3645 3869 ENDDO 3646 ELSE 3647 DO jbgp = jb, ny +nbgp3870 ELSEIF ( edge == 'n' ) THEN 3871 DO jbgp = jb, ny+nbgp 3648 3872 f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg) 3649 3873 ENDDO … … 3654 3878 3655 3879 3656 SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, var ) 3880 SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, & 3881 r2z, var ) 3657 3882 3658 3883 ! … … 3660 3885 !-- conditions. This subroutine handles the top boundary. 3661 3886 !-- This subroutine is based on trilinear interpolation. 3662 !-- Constant dz is still assumed. 3887 3663 3888 IMPLICIT NONE 3664 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !: 3665 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc !: 3666 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !: 3667 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !: 3668 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !: 3669 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !: 3670 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !: 3671 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !: 3889 3890 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 3891 INTENT(INOUT) :: f !: 3892 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), & 3893 INTENT(IN) :: fc !: 3894 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !: 3895 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !: 3896 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !: 3897 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !: 3898 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !: 3899 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !: 3672 3900 3673 INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN) 3674 INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) 3675 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) 3901 INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN) :: ic !: 3902 INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) :: jc !: 3903 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !: 3676 3904 3677 3905 CHARACTER(LEN=1), INTENT(IN) :: var !: … … 3701 3929 ENDIF 3702 3930 3703 DO i = nxl - 1, nxr +13704 DO j = nys - 1, nyn +13931 DO i = nxl-1, nxr+1 3932 DO j = nys-1, nyn+1 3705 3933 l = ic(i) 3706 3934 m = jc(j) … … 3724 3952 ! 3725 3953 !-- Rescale if f is the TKE. 3726 !-- It is assumed that the bottom surface never reaches the top 3727 !-- - boundary of a nest domain.3728 IF ( var == 'e' ) THEN3954 !-- It is assumed that the bottom surface never reaches the top boundary of a 3955 !-- nest domain. 3956 IF ( var == 'e' ) THEN 3729 3957 DO i = nxl, nxr 3730 3958 DO j = nys, nyn … … 3746 3974 !-- subroutine handles the left and right boundaries. However, this operation 3747 3975 !-- is only needed in case of one-way coupling. 3976 3748 3977 IMPLICIT NONE 3749 3978 … … 3785 4014 ENDIF 3786 4015 3787 DO j = nys, nyn +13788 DO k = kb(j,i), nzt 4016 DO j = nys, nyn+1 4017 DO k = kb(j,i), nzt+1 3789 4018 vdotnor = outnor * u(k,j,ied) 3790 4019 ! … … 3806 4035 ENDDO 3807 4036 ELSEIF ( edge == 'r' ) THEN 3808 DO ibgp = ib, nx +nbgp4037 DO ibgp = ib, nx+nbgp 3809 4038 f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i) 3810 4039 ENDDO … … 3822 4051 !-- interpolated values by values extrapolated from the domain. This 3823 4052 !-- subroutine handles the south and north boundaries. 4053 3824 4054 IMPLICIT NONE 3825 4055 … … 3861 4091 ENDIF 3862 4092 3863 DO i = nxl, nxr +13864 DO k = kb(j,i), nzt +14093 DO i = nxl, nxr+1 4094 DO k = kb(j,i), nzt+1 3865 4095 vdotnor = outnor * v(k,jed,i) 3866 4096 ! … … 3878 4108 !-- Store the boundary values also into the redundant ghost node layers. 3879 4109 IF ( edge == 's' ) THEN 3880 DO jbgp = -nbgp, jb4110 DO jbgp = -nbgp, jb 3881 4111 f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg) 3882 4112 ENDDO 3883 4113 ELSEIF ( edge == 'n' ) THEN 3884 DO jbgp = jb, ny +nbgp4114 DO jbgp = jb, ny+nbgp 3885 4115 f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg) 3886 4116 ENDDO … … 3896 4126 !-- conditions. This subroutine handles the top boundary. It is based on 3897 4127 !-- trilinear interpolation. 4128 3898 4129 IMPLICIT NONE 3899 4130 … … 3941 4172 3942 4173 SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu, & 3943 kb,var )4174 var ) 3944 4175 ! 3945 4176 !-- Anterpolation of internal-node values to be used as the server-domain … … 3947 4178 !-- integration of the fine-grid values contained within the coarse-grid 3948 4179 !-- cell. 4180 3949 4181 IMPLICIT NONE 3950 4182 … … 3968 4200 INTEGER(iwp), INTENT(IN) :: kct !: 3969 4201 3970 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifl !: 3971 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifu !: 3972 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !: 3973 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !: 3974 !-- TO_DO: is the next line really unnecessary? 3975 INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) :: kb !: may be unnecessary 3976 INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !: 3977 INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !: 4202 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifl !: 4203 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifu !: 4204 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !: 4205 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !: 4206 INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !: 4207 INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !: 3978 4208 3979 4209 … … 4049 4279 ENDIF 4050 4280 ENDIF 4051 4052 fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) + & 4281 ! 4282 !-- TO DO: introduce 3-d coarse grid array for precomputed 4283 !-- 1/REAL(nfc) values 4284 fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) + & 4053 4285 fra * cellsum / REAL( nfc, KIND = wp ) 4054 4286
Note: See TracChangeset
for help on using the changeset viewer.
