Changeset 2602
- Timestamp:
- Nov 3, 2017 11:06:41 AM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pmc_interface_mod.f90
r2599 r2602 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Some cleanup and commenting improvements only. 28 ! Index-limit related bug (occurred with nesting_mode='vertical') fixed in 29 ! pmci_interp_tril_t. Check for too high nest domain added in pmci_setup_parent. 30 ! Some cleaning up made. 29 31 ! 30 32 ! 2582 2017-10-26 13:19:46Z hellstea … … 159 161 ! only the total number of PEs is given for the domains, npe_x and npe_y 160 162 ! replaced by npe_total, two unused elements removed from array 161 ! define_coarse_grid_real,163 ! parent_grid_info_real, 162 164 ! array management changed from linked list to sequential loop 163 165 ! … … 258 260 259 261 PRIVATE 260 261 262 ! 262 263 !-- Constants 263 264 INTEGER(iwp), PARAMETER :: child_to_parent = 2 !: 264 265 INTEGER(iwp), PARAMETER :: parent_to_child = 1 !: 265 266 266 ! 267 267 !-- Coupler setup … … 271 271 INTEGER(iwp), SAVE :: cpl_npe_total !: 272 272 INTEGER(iwp), SAVE :: cpl_parent_id !: 273 274 273 ! 275 274 !-- Control parameters, will be made input parameters later … … 287 286 REAL(wp), SAVE :: anterp_relax_length_n = -1.0_wp !: 288 287 REAL(wp), SAVE :: anterp_relax_length_t = -1.0_wp !: 289 290 288 ! 291 289 !-- Geometry … … 295 293 REAL(wp), SAVE :: lower_left_coord_x !: 296 294 REAL(wp), SAVE :: lower_left_coord_y !: 297 298 295 ! 299 296 !-- Child coarse data arrays … … 321 318 REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ncc !: 322 319 REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: sc !: 323 324 320 ! 325 321 !-- Child interpolation coefficients and child-array indices to be … … 343 339 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: r1zw !: 344 340 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: r2zw !: 345 346 341 ! 347 342 !-- Child index arrays and log-ratio arrays for the log-law near-wall … … 372 367 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: logc_ratio_w_r !: 373 368 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: logc_ratio_w_s !: 374 375 369 ! 376 370 !-- Upper bounds for k in anterpolation. 377 371 INTEGER(iwp), SAVE :: kctu !: 378 372 INTEGER(iwp), SAVE :: kctw !: 379 380 373 ! 381 374 !-- Upper bound for k in log-law correction in interpolation. … … 384 377 INTEGER(iwp), SAVE :: nzt_topo_nestbc_r !: 385 378 INTEGER(iwp), SAVE :: nzt_topo_nestbc_s !: 386 387 379 ! 388 380 !-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation. … … 391 383 INTEGER(iwp), SAVE :: nhls !: 392 384 INTEGER(iwp), SAVE :: nhln !: 393 394 385 ! 395 386 !-- Spatial under-relaxation coefficients for anterpolation. … … 397 388 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: fray !: 398 389 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: fraz !: 399 400 390 ! 401 391 !-- Child-array indices to be precomputed and stored for anterpolation. … … 412 402 INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) :: kflo !: 413 403 INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) :: kfuo !: 414 415 404 ! 416 405 !-- Number of fine-grid nodes inside coarse-grid ij-faces … … 422 411 INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) :: kfc_s !: 423 412 424 INTEGER(iwp), DIMENSION(3) :: define_coarse_grid_int !:425 REAL(wp), DIMENSION(7) :: define_coarse_grid_real !:413 INTEGER(iwp), DIMENSION(3) :: parent_grid_info_int !: 414 REAL(wp), DIMENSION(7) :: parent_grid_info_real !: 426 415 427 416 TYPE coarsegrid_def … … 535 524 536 525 ENDIF 537 538 526 ! 539 527 !-- Check steering parameter values … … 554 542 CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 ) 555 543 ENDIF 556 557 544 ! 558 545 !-- Set the general steering switch which tells PALM that its a nested run 559 546 nested_run = .TRUE. 560 561 547 ! 562 548 !-- Get some variables required by the pmc-interface (and in some cases in the … … 650 636 REAL(wp), DIMENSION(:), ALLOCATABLE :: ch_ys !: 651 637 REAL(wp), DIMENSION(:), ALLOCATABLE :: ch_yn !: 638 REAL(wp) :: cl_height !: 652 639 REAL(wp) :: dx_cl !: 653 640 REAL(wp) :: dy_cl !: … … 663 650 REAL(wp), DIMENSION(:), ALLOCATABLE :: cl_coord_x !: 664 651 REAL(wp), DIMENSION(:), ALLOCATABLE :: cl_coord_y !: 665 666 652 667 653 ! 668 654 ! Initialize the pmc parent 669 655 CALL pmc_parentinit 670 671 656 ! 672 657 !-- Corners of all children of the present parent … … 677 662 ALLOCATE( ch_yn(1:SIZE( pmc_parent_for_child ) - 1) ) 678 663 ENDIF 679 680 664 ! 681 665 !-- Get coordinates from all children … … 692 676 dx_cl = val(4) 693 677 dy_cl = val(5) 694 678 cl_height = fval(1) 695 679 nz_cl = nz 696 697 680 ! 698 681 !-- Find the highest nest level in the coarse grid for the reduced z … … 704 687 ENDIF 705 688 ENDDO 706 707 689 ! 708 690 !-- Get absolute coordinates from the child … … 711 693 712 694 CALL pmc_recv_from_child( child_id, cl_coord_x, SIZE( cl_coord_x ), & 713 695 0, 11, ierr ) 714 696 CALL pmc_recv_from_child( child_id, cl_coord_y, SIZE( cl_coord_y ), & 715 0, 12, ierr ) 716 ! WRITE ( 0, * ) 'receive from pmc child ', child_id, nx_cl, ny_cl 697 0, 12, ierr ) 717 698 718 define_coarse_grid_real(1) = lower_left_coord_x 719 define_coarse_grid_real(2) = lower_left_coord_y 720 define_coarse_grid_real(3) = dx 721 define_coarse_grid_real(4) = dy 722 define_coarse_grid_real(5) = lower_left_coord_x + ( nx + 1 ) * dx 723 define_coarse_grid_real(6) = lower_left_coord_y + ( ny + 1 ) * dy 724 define_coarse_grid_real(7) = dz 725 726 define_coarse_grid_int(1) = nx 727 define_coarse_grid_int(2) = ny 728 define_coarse_grid_int(3) = nz_cl 729 699 parent_grid_info_real(1) = lower_left_coord_x 700 parent_grid_info_real(2) = lower_left_coord_y 701 parent_grid_info_real(3) = dx 702 parent_grid_info_real(4) = dy 703 parent_grid_info_real(5) = lower_left_coord_x + ( nx + 1 ) * dx 704 parent_grid_info_real(6) = lower_left_coord_y + ( ny + 1 ) * dy 705 parent_grid_info_real(7) = dz 706 707 parent_grid_info_int(1) = nx 708 parent_grid_info_int(2) = ny 709 parent_grid_info_int(3) = nz_cl 730 710 ! 731 711 !-- Check that the child domain matches parent domain. 732 712 nomatch = 0 733 713 IF ( nesting_mode == 'vertical' ) THEN 734 right_limit = define_coarse_grid_real(5)735 north_limit = define_coarse_grid_real(6)714 right_limit = parent_grid_info_real(5) 715 north_limit = parent_grid_info_real(6) 736 716 IF ( ( cl_coord_x(nx_cl+1) /= right_limit ) .OR. & 737 717 ( cl_coord_y(ny_cl+1) /= north_limit ) ) THEN 738 718 nomatch = 1 739 719 ENDIF 740 ELSE 741 720 ELSE 742 721 ! 743 722 !-- Check that the child domain is completely inside the parent domain. … … 745 724 yez = ( nbgp + 1 ) * dy 746 725 left_limit = lower_left_coord_x + xez 747 right_limit = define_coarse_grid_real(5) - xez726 right_limit = parent_grid_info_real(5) - xez 748 727 south_limit = lower_left_coord_y + yez 749 north_limit = define_coarse_grid_real(6) - yez728 north_limit = parent_grid_info_real(6) - yez 750 729 IF ( ( cl_coord_x(0) < left_limit ) .OR. & 751 730 ( cl_coord_x(nx_cl+1) > right_limit ) .OR. & … … 755 734 ENDIF 756 735 ENDIF 757 736 ! 737 !-- Child domain must be lower than the parent domain such 738 !-- that the top ghost layer of the child grid does not exceed 739 !-- the parent domain top boundary. 740 IF ( cl_height > zw(nz) ) THEN 741 nomatch = 1 742 ENDIF 758 743 ! 759 744 !-- Check that parallel nest domains, if any, do not overlap. … … 779 764 DEALLOCATE( cl_coord_x ) 780 765 DEALLOCATE( cl_coord_y ) 781 782 766 ! 783 767 !-- Send coarse grid information to child 784 CALL pmc_send_to_child( child_id, define_coarse_grid_real,&785 SIZE( define_coarse_grid_real ), 0, 21,&768 CALL pmc_send_to_child( child_id, parent_grid_info_real, & 769 SIZE( parent_grid_info_real ), 0, 21, & 786 770 ierr ) 787 CALL pmc_send_to_child( child_id, define_coarse_grid_int, 3, 0,&771 CALL pmc_send_to_child( child_id, parent_grid_info_int, 3, 0, & 788 772 22, ierr ) 789 790 773 ! 791 774 !-- Send local grid to child … … 794 777 CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25, & 795 778 ierr ) 796 797 779 ! 798 780 !-- Also send the dzu-, dzw-, zu- and zw-arrays here … … 818 800 819 801 CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr ) 820 821 802 ! 822 803 !-- TO_DO: Klaus: please give a comment what is done here 823 804 CALL pmci_create_index_list 824 825 805 ! 826 806 !-- Include couple arrays into parent content … … 870 850 871 851 IF ( myid == 0 ) THEN 852 ! 872 853 !-- TO_DO: Klaus: give more specific comment what size_of_array stands for 873 854 CALL pmc_recv_from_child( child_id, size_of_array, 2, 0, 40, ierr ) … … 875 856 CALL pmc_recv_from_child( child_id, coarse_bound_all, & 876 857 SIZE( coarse_bound_all ), 0, 41, ierr ) 877 878 858 ! 879 859 !-- Compute size of index_list. … … 1023 1003 1024 1004 CALL pmc_set_dataarray_name( lastentry = .TRUE. ) 1025 1026 1005 ! 1027 1006 !-- Send grid to parent … … 1039 1018 CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr ) 1040 1019 CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr ) 1041 1042 1020 ! 1043 1021 !-- Receive Coarse grid information. 1044 !-- TO_DO: find shorter and more meaningful name for define_coarse_grid_real 1045 CALL pmc_recv_from_parent( define_coarse_grid_real, & 1046 SIZE(define_coarse_grid_real), 0, 21, ierr ) 1047 CALL pmc_recv_from_parent( define_coarse_grid_int, 3, 0, 22, ierr ) 1022 CALL pmc_recv_from_parent( parent_grid_info_real, & 1023 SIZE(parent_grid_info_real), 0, 21, ierr ) 1024 CALL pmc_recv_from_parent( parent_grid_info_int, 3, 0, 22, ierr ) 1048 1025 ! 1049 1026 !-- Debug-printouts - keep them 1050 1027 ! WRITE(0,*) 'Coarse grid from parent ' 1051 ! WRITE(0,*) 'startx_tot = ', define_coarse_grid_real(1)1052 ! WRITE(0,*) 'starty_tot = ', define_coarse_grid_real(2)1053 ! WRITE(0,*) 'endx_tot = ', define_coarse_grid_real(5)1054 ! WRITE(0,*) 'endy_tot = ', define_coarse_grid_real(6)1055 ! WRITE(0,*) 'dx = ', define_coarse_grid_real(3)1056 ! WRITE(0,*) 'dy = ', define_coarse_grid_real(4)1057 ! WRITE(0,*) 'dz = ', define_coarse_grid_real(7)1058 ! WRITE(0,*) 'nx_coarse = ', define_coarse_grid_int(1)1059 ! WRITE(0,*) 'ny_coarse = ', define_coarse_grid_int(2)1060 ! WRITE(0,*) 'nz_coarse = ', define_coarse_grid_int(3)1061 ENDIF 1062 1063 CALL MPI_BCAST( define_coarse_grid_real, SIZE(define_coarse_grid_real),&1028 ! WRITE(0,*) 'startx_tot = ',parent_grid_info_real(1) 1029 ! WRITE(0,*) 'starty_tot = ',parent_grid_info_real(2) 1030 ! WRITE(0,*) 'endx_tot = ',parent_grid_info_real(5) 1031 ! WRITE(0,*) 'endy_tot = ',parent_grid_info_real(6) 1032 ! WRITE(0,*) 'dx = ',parent_grid_info_real(3) 1033 ! WRITE(0,*) 'dy = ',parent_grid_info_real(4) 1034 ! WRITE(0,*) 'dz = ',parent_grid_info_real(7) 1035 ! WRITE(0,*) 'nx_coarse = ',parent_grid_info_int(1) 1036 ! WRITE(0,*) 'ny_coarse = ',parent_grid_info_int(2) 1037 ! WRITE(0,*) 'nz_coarse = ',parent_grid_info_int(3) 1038 ENDIF 1039 1040 CALL MPI_BCAST( parent_grid_info_real, SIZE(parent_grid_info_real), & 1064 1041 MPI_REAL, 0, comm2d, ierr ) 1065 CALL MPI_BCAST( define_coarse_grid_int, 3, MPI_INTEGER, 0, comm2d, ierr ) 1066 1067 cg%dx = define_coarse_grid_real(3) 1068 cg%dy = define_coarse_grid_real(4) 1069 cg%dz = define_coarse_grid_real(7) 1070 cg%nx = define_coarse_grid_int(1) 1071 cg%ny = define_coarse_grid_int(2) 1072 cg%nz = define_coarse_grid_int(3) 1073 1042 CALL MPI_BCAST( parent_grid_info_int, 3, MPI_INTEGER, 0, comm2d, ierr ) 1043 1044 cg%dx = parent_grid_info_real(3) 1045 cg%dy = parent_grid_info_real(4) 1046 cg%dz = parent_grid_info_real(7) 1047 cg%nx = parent_grid_info_int(1) 1048 cg%ny = parent_grid_info_int(2) 1049 cg%nz = parent_grid_info_int(3) 1074 1050 ! 1075 1051 !-- Get parent coordinates on coarse grid … … 1081 1057 ALLOCATE( cg%zu(0:cg%nz+1) ) 1082 1058 ALLOCATE( cg%zw(0:cg%nz+1) ) 1083 1084 1059 ! 1085 1060 !-- Get coarse grid coordinates and values of the z-direction from the parent … … 1094 1069 1095 1070 ENDIF 1096 1097 1071 ! 1098 1072 !-- Broadcast this information … … 1110 1084 !-- TO_DO: Klaus give a comment what is happening here 1111 1085 CALL pmc_c_get_2d_index_list 1112 1113 1086 ! 1114 1087 !-- Include couple arrays into child content … … 1116 1089 CALL pmc_c_clear_next_array_list 1117 1090 DO WHILE ( pmc_c_getnextarray( myname ) ) 1118 !-- TO_DO: Klaus, why the child-arrays are still up to cg%nz?? 1091 !-- Note that cg%nz is not th eoriginal nz of parent, but the highest 1092 !-- parent-grid level needed for nesting. 1119 1093 CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz ) 1120 1094 ENDDO 1121 1095 CALL pmc_c_setind_and_allocmem 1122 1123 1096 ! 1124 1097 !-- Precompute interpolation coefficients and child-array indices 1125 1098 CALL pmci_init_interp_tril 1126 1127 1099 ! 1128 1100 !-- Precompute the log-law correction index- and ratio-arrays 1129 1101 CALL pmci_init_loglaw_correction 1130 1131 1102 ! 1132 1103 !-- Define the SGS-TKE scaling factor based on the grid-spacing ratio 1133 1104 CALL pmci_init_tkefactor 1134 1135 1105 ! 1136 1106 !-- Two-way coupling for general and vertical nesting. … … 1141 1111 CALL pmci_init_anterp_tophat 1142 1112 ENDIF 1143 1144 1113 ! 1145 1114 !-- Finally, compute the total area of the top-boundary face of the domain. … … 1150 1119 1151 1120 CONTAINS 1121 1152 1122 1153 1123 SUBROUTINE pmci_map_fine_to_coarse_grid … … 1263 1233 INTEGER(iwp) :: k !: 1264 1234 INTEGER(iwp) :: kc !: 1235 INTEGER(iwp) :: kdzo !: 1236 INTEGER(iwp) :: kdzw !: 1265 1237 1266 1238 REAL(wp) :: xb !: … … 1301 1273 ALLOCATE( r1zo(nzb:nzt+1) ) 1302 1274 ALLOCATE( r2zo(nzb:nzt+1) ) 1303 1304 1275 ! 1305 1276 !-- Note that the node coordinates xfs... and xcs... are relative to the … … 1336 1307 zfso = zu(k) 1337 1308 1338 kc = 0 1339 DO WHILE ( cg%zw(kc) <= zfsw ) 1340 kc = kc + 1 1309 DO kc = 0, cg%nz+1 1310 IF ( cg%zw(kc) > zfsw ) EXIT 1341 1311 ENDDO 1342 1312 kcw(k) = kc - 1 1343 1344 kc = 0 1345 DO WHILE ( cg%zu(kc) <= zfso ) 1346 kc = kc + 1 1313 1314 !if ( myid == 0 .and. nx==191 ) then 1315 ! write(162,*)nx, nzt+1, k, zw(k), cg%nz+1, kcw(k), cg%zw(kcw(k)) 1316 !endif 1317 !if ( myid == 0 .and. nx==383 ) then 1318 ! write(163,*)nx, nzt+1, k, zw(k), cg%nz+1, kcw(k), cg%zw(kcw(k)) 1319 !endif 1320 1321 DO kc = 0, cg%nz+1 1322 IF ( cg%zu(kc) > zfso ) EXIT 1347 1323 ENDDO 1348 1324 kco(k) = kc - 1 … … 1350 1326 zcsw = cg%zw(kcw(k)) 1351 1327 zcso = cg%zu(kco(k)) 1352 r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k)+1) 1353 r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k)+1) 1328 kdzw = MIN( kcw(k)+1, cg%nz+1 ) 1329 kdzo = MIN( kco(k)+1, cg%nz+1 ) 1330 r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kdzw) 1331 r2zo(k) = ( zfso - zcso ) / cg%dzu(kdzo) 1354 1332 r1zw(k) = 1.0_wp - r2zw(k) 1355 1333 r1zo(k) = 1.0_wp - r2zo(k) … … 1507 1485 nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1 1508 1486 ENDIF 1509 1510 1487 ! 1511 1488 !-- Then determine the maximum number of near-wall nodes per wall point based … … 1513 1490 nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r, & 1514 1491 nzt_topo_nestbc_s, nzt_topo_nestbc_n ) 1515 1516 1492 ! 1517 1493 !-- Note that the outer division must be integer division. … … 1529 1505 1530 1506 z0_topo = roughness_length 1531 1532 1507 ! 1533 1508 !-- First horizontal walls. Note that also logc_w_? and logc_ratio_w_? need to … … 1590 1565 1591 1566 ENDIF 1592 1593 1567 ! 1594 1568 !-- Right boundary … … 1650 1624 1651 1625 ENDIF 1652 1653 1626 ! 1654 1627 !-- South boundary … … 1707 1680 1708 1681 ENDIF 1709 1710 1682 ! 1711 1683 !-- North boundary … … 1763 1735 1764 1736 ENDIF 1765 1766 1737 ! 1767 1738 !-- Then vertical walls and corners if necessary … … 1817 1788 lcr(0:ncorr-1) = 1.0_wp 1818 1789 ENDIF 1819 1820 1790 ! 1821 1791 !-- Wall for u on the north side, but not on the south side … … 1851 1821 lcr(0:ncorr-1) = 1.0_wp 1852 1822 ENDIF 1853 1854 1823 ! 1855 1824 !-- Wall for w on the north side, but not on the south side. … … 1872 1841 1873 1842 ENDIF ! IF ( nest_bound_l ) 1874 1875 1843 ! 1876 1844 !-- Right boundary … … 1907 1875 lcr(0:ncorr-1) = 1.0_wp 1908 1876 ENDIF 1909 1910 1877 ! 1911 1878 !-- Wall for u on the north side, but not on the south side … … 1923 1890 lcr(0:ncorr-1) = 1.0_wp 1924 1891 ENDIF 1925 1926 1892 ! 1927 1893 !-- Wall for w on the south side, but not on the north side … … 1939 1905 lcr(0:ncorr-1) = 1.0_wp 1940 1906 ENDIF 1941 1942 1907 ! 1943 1908 !-- Wall for w on the north side, but not on the south side … … 1948 1913 CALL pmci_define_loglaw_correction_parameters( lc, lcr, & 1949 1914 k, j, inc, wall_index, z0_topo, kb, direction, ncorr ) 1950 1951 1915 ! 1952 1916 !-- The direction of the wall-normal index is stored as the … … 1961 1925 1962 1926 ENDIF ! IF ( nest_bound_r ) 1963 1964 1927 ! 1965 1928 !-- South boundary … … 1996 1959 lcr(0:ncorr-1) = 1.0_wp 1997 1960 ENDIF 1998 1999 1961 ! 2000 1962 !-- Wall for v on the right side, but not on the left side … … 2013 1975 lcr(0:ncorr-1) = 1.0_wp 2014 1976 ENDIF 2015 2016 1977 ! 2017 1978 !-- Wall for w on the left side, but not on the right side … … 2052 2013 2053 2014 ENDIF ! IF (nest_bound_s ) 2054 2055 2015 ! 2056 2016 !-- North boundary … … 2087 2047 lcr(0:ncorr-1) = 1.0_wp 2088 2048 ENDIF 2089 2090 2049 ! 2091 2050 !-- Wall for v on the right side, but not on the left side … … 2103 2062 lcr(0:ncorr-1) = 1.0_wp 2104 2063 ENDIF 2105 2106 2064 ! 2107 2065 !-- Wall for w on the left side, but not on the right side … … 2119 2077 lcr(0:ncorr-1) = 1.0_wp 2120 2078 ENDIF 2121 2122 2079 ! 2123 2080 !-- Wall for w on the right side, but not on the left side … … 2204 2161 z0_l, inc ) 2205 2162 ENDIF 2206 2207 2163 ! 2208 2164 !-- The role of inc here is to make the comparison operation "<" … … 2311 2267 REAL(wp) :: logyc1 !: 2312 2268 REAL(wp) :: yc1 !: 2313 2269 2314 2270 ! 2315 2271 !-- yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from … … 2416 2372 anterp_relax_length_t = 0.1_wp * zu(nzt) 2417 2373 ENDIF 2418 2419 2374 ! 2420 2375 !-- First determine kctu and kctw that are the coarse-grid upper bounds for … … 2472 2427 iflu(icr) = nxrg 2473 2428 ifuu(icr) = nxrg 2474 2475 2429 ! 2476 2430 !-- i-indices of others for each ii-index value … … 2495 2449 iflo(icr) = nxrg 2496 2450 ifuo(icr) = nxrg 2497 2498 2451 ! 2499 2452 !-- j-indices of v for each jj-index value … … 2540 2493 jflo(jcn) = nyng 2541 2494 jfuo(jcn) = nyng 2542 2543 2495 ! 2544 2496 !-- k-indices of w for each kk-index value … … 2558 2510 kstart = kflw(kk) 2559 2511 ENDDO 2560 2561 2512 ! 2562 2513 !-- k-indices of others for each kk-index value … … 2576 2527 kstart = kflo(kk) 2577 2528 ENDDO 2578 2579 2529 ! 2580 2530 !-- Precomputation of number of fine-grid nodes inside coarse-grid ij-faces. … … 2645 2595 2646 2596 2647 SUBROUTINE pmci_init_tkefactor2597 SUBROUTINE pmci_init_tkefactor 2648 2598 2649 2599 ! … … 2712 2662 ENDIF 2713 2663 2714 IF ( nest_bound_s ) THEN2664 IF ( nest_bound_s ) THEN 2715 2665 ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) ) 2716 2666 tkefactor_s = 0.0_wp … … 2718 2668 DO i = nxlg, nxrg 2719 2669 k_wall = get_topography_top_index( j, i, 's' ) 2720 2670 2721 2671 DO k = k_wall + 1, nzt 2722 2672 2723 2673 kc = kco(k+1) 2724 2674 glsf = ( dx * dy * dzu(k) )**p13 … … 2727 2677 fw = EXP( -cfw*height / glsf ) 2728 2678 tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * & 2729 2679 ( glsf / glsc )**p23 ) 2730 2680 ENDDO 2731 2681 tkefactor_s(k_wall,i) = c_tkef * fw0 … … 2827 2777 NULLIFY( p_3d ) 2828 2778 NULLIFY( p_2d ) 2829 2830 2779 ! 2831 2780 !-- List of array names, which can be coupled. … … 2928 2877 NULLIFY( p_3d ) 2929 2878 NULLIFY( p_2d ) 2930 2931 2879 ! 2932 2880 !-- List of array names, which can be coupled … … 3036 2984 3037 2985 ! 3038 !-- Root id is never achild2986 !-- Root model is never anyone's child 3039 2987 IF ( cpl_id > 1 ) THEN 3040 3041 2988 ! 3042 2989 !-- Child domain boundaries in the parent index space … … 3045 2992 jcs = coarse_bound(3) 3046 2993 jcn = coarse_bound(4) 3047 3048 2994 ! 3049 2995 !-- Get data from the parent 3050 2996 CALL pmc_c_getbuffer( waittime = waittime ) 3051 3052 2997 ! 3053 2998 !-- The interpolation. … … 3203 3148 ENDIF 3204 3149 ENDIF 3205 3206 3150 ! 3207 3151 !-- Trilinear interpolation. … … 3222 3166 ENDDO 3223 3167 ENDDO 3224 3225 3168 ! 3226 3169 !-- Correct the interpolated values of u and v in near-wall nodes, i.e. in … … 3321 3264 ENDIF 3322 3265 ENDIF 3323 3324 3266 ! 3325 3267 !-- Same for restart time … … 3338 3280 ENDIF 3339 3281 ENDIF 3340 3341 3282 ! 3342 3283 !-- Same for dt_restart … … 3355 3296 ENDIF 3356 3297 ENDIF 3357 3358 3298 ! 3359 3299 !-- Same for time_restart … … 3432 3372 volume_flow(1) = volume_flow_l(1) 3433 3373 #endif 3434 3374 3435 3375 ! 3436 3376 !-- Sum up the volume flow through the south/north boundaries … … 3512 3452 !-- MPI_ALLREDUCE with the MPI_MIN operator over all processes using 3513 3453 !-- the global communicator MPI_COMM_WORLD. 3454 3514 3455 IMPLICIT NONE 3515 3456 … … 3605 3546 3606 3547 END SUBROUTINE pmci_datatrans 3607 3608 3548 3609 3549 … … 3642 3582 CALL pmc_s_getdata_from_buffer( child_id ) 3643 3583 CALL cpu_log( log_point_s(72), 'pmc parent recv', 'stop' ) 3644 3645 3584 ! 3646 3585 !-- The anterpolated data is now available in u etc 3647 3586 IF ( topography /= 'flat' ) THEN 3648 3649 3587 ! 3650 3588 !-- Inside buildings/topography reset velocities back to zero. … … 3732 3670 ENDIF 3733 3671 3734 CONTAINS 3735 3672 CONTAINS 3673 3674 3736 3675 SUBROUTINE pmci_interpolation 3737 3676 3738 3677 ! 3739 3678 !-- A wrapper routine for all interpolation and extrapolation actions 3679 3740 3680 IMPLICIT NONE 3741 3681 3742 3682 ! 3743 3683 !-- In case of vertical nesting no interpolation is needed for the … … 3746 3686 3747 3687 ! 3748 !-- Left border pe:3688 !-- Left border pe: 3749 3689 IF ( nest_bound_l ) THEN 3750 3690 CALL pmci_interp_tril_lr( u, uc, icu, jco, kco, r1xu, r2xu, & … … 3855 3795 3856 3796 ENDIF 3857 3858 ! 3859 !-- Right border pe 3797 ! 3798 !-- Right border pe 3860 3799 IF ( nest_bound_r ) THEN 3861 3800 CALL pmci_interp_tril_lr( u, uc, icu, jco, kco, r1xu, r2xu, & … … 3969 3908 3970 3909 ENDIF 3971 3972 ! 3973 !-- South border pe 3910 ! 3911 !-- South border pe 3974 3912 IF ( nest_bound_s ) THEN 3975 3913 CALL pmci_interp_tril_sn( u, uc, icu, jco, kco, r1xu, r2xu, & … … 4077 4015 4078 4016 ENDIF 4079 4080 ! 4081 !-- North border pe 4017 ! 4018 !-- North border pe 4082 4019 IF ( nest_bound_n ) THEN 4083 4020 CALL pmci_interp_tril_sn( u, uc, icu, jco, kco, r1xu, r2xu, & … … 4189 4126 4190 4127 ENDIF ! IF ( nesting_mode /= 'vertical' ) 4191 4192 4128 ! 4193 4129 !-- All PEs are top-border PEs … … 4413 4349 4414 4350 DO j = nys, nyn+1 4415 DO k = nzb, nzt+1 4351 ! 4352 !-- Determine vertical index of topography top at grid point (j,i) 4353 k_wall = get_topography_top_index( j, i, TRIM( var ) ) 4354 4355 DO k = k_wall, nzt+1 4416 4356 l = ic(i) 4417 4357 m = jc(j) … … 4426 4366 ENDDO 4427 4367 ENDDO 4428 4429 4368 ! 4430 4369 !-- Generalized log-law-correction algorithm. … … 4450 4389 ENDDO 4451 4390 ENDIF 4452 4453 4391 ! 4454 4392 !-- In case of non-flat topography, also vertical walls and corners need to be … … 4513 4451 4514 4452 ENDIF ! ( topography /= 'flat' ) 4515 4516 4453 ! 4517 4454 !-- Rescale if f is the TKE. … … 4539 4476 ENDIF 4540 4477 ENDIF 4541 4542 4478 ! 4543 4479 !-- Store the boundary values also into the other redundant ghost node layers … … 4634 4570 ENDIF 4635 4571 4636 4637 4572 DO i = nxl, nxr+1 4638 4573 ! … … 4653 4588 ENDDO 4654 4589 ENDDO 4655 4656 4590 ! 4657 4591 !-- Generalized log-law-correction algorithm. … … 4677 4611 ENDDO 4678 4612 ENDIF 4679 4680 4613 ! 4681 4614 !-- In case of non-flat topography, also vertical walls and corners need to be … … 4741 4674 4742 4675 ENDIF ! ( topography /= 'flat' ) 4743 4744 4676 ! 4745 4677 !-- Rescale if f is the TKE. … … 4765 4697 ENDIF 4766 4698 ENDIF 4767 4768 4699 ! 4769 4700 !-- Store the boundary values also into the other redundant ghost node layers … … 4810 4741 4811 4742 INTEGER(iwp) :: i !: 4743 INTEGER(iwp) :: ib !: 4744 INTEGER(iwp) :: ie !: 4812 4745 INTEGER(iwp) :: j !: 4746 INTEGER(iwp) :: jb !: 4747 INTEGER(iwp) :: je !: 4813 4748 INTEGER(iwp) :: k !: 4814 4749 INTEGER(iwp) :: l !: … … 4832 4767 k = nzt + 1 4833 4768 ENDIF 4834 4835 DO i = nxl-1, nxr+1 4836 DO j = nys-1, nyn+1 4769 ! 4770 !-- These exceedings by one are needed only to avoid stripes 4771 !-- and spots in visualization. They have no effect on the 4772 !-- actual solution. 4773 ib = nxl-1 4774 ie = nxr+1 4775 jb = nys-1 4776 je = nyn+1 4777 ! 4778 !-- The exceedings must not be made past the outer edges in 4779 !-- case of pure vertical nesting. 4780 IF ( nesting_mode == 'vertical' ) THEN 4781 IF ( nxl == 0 ) ib = nxl 4782 IF ( nxr == nx ) ie = nxr 4783 IF ( nys == 0 ) jb = nys 4784 IF ( nyn == ny ) je = nyn 4785 ENDIF 4786 4787 DO i = ib, ie 4788 DO j = jb, je 4837 4789 l = ic(i) 4838 4790 m = jc(j) 4839 n = kc(k) 4791 n = kc(k) 4840 4792 fkj = r1x(i) * fc(n,m,l) + r2x(i) * fc(n,m,l+1) 4841 4793 fkjp = r1x(i) * fc(n,m+1,l) + r2x(i) * fc(n,m+1,l+1) … … 4847 4799 ENDDO 4848 4800 ENDDO 4849 4850 4801 ! 4851 4802 !-- Just fill up the second ghost-node layer for w. … … 4853 4804 f(nzt+1,:,:) = f(nzt,:,:) 4854 4805 ENDIF 4855 4856 4806 ! 4857 4807 !-- Rescale if f is the TKE. … … 4896 4846 4897 4847 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !: 4898 4899 4848 ! 4900 4849 !-- Check which edge is to be handled: left or right … … 4917 4866 ENDIF 4918 4867 4919 4920 4868 DO j = nys, nyn+1 4921 4869 ! … … 4935 4883 ENDIF 4936 4884 ENDDO 4937 4938 4885 ! 4939 4886 !-- Store the boundary values also into the redundant ghost node layers. … … 4998 4945 ENDIF 4999 4946 5000 5001 4947 DO i = nxl, nxr+1 5002 4948 ! … … 5016 4962 ENDIF 5017 4963 ENDDO 5018 5019 4964 ! 5020 4965 !-- Store the boundary values also into the redundant ghost node layers. … … 5072 5017 ENDDO 5073 5018 ENDDO 5074 5075 5019 ! 5076 5020 !-- Just fill up the second ghost-node layer for w … … 5126 5070 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT) :: fc !< Treated variable - parent domain 5127 5071 5128 5129 5072 ! 5130 5073 !-- Initialize the index bounds for anterpolation … … 5179 5122 jcnm = jcn - nhln - 1 5180 5123 ENDIF 5181 ENDIF 5182 5124 ENDIF 5183 5125 ! 5184 5126 !-- Note that ii, jj, and kk are coarse-grid indices and i,j, and k … … 5211 5153 ENDDO 5212 5154 ENDDO 5213 5214 5155 5215 5156 END SUBROUTINE pmci_anterp_tophat … … 5298 5239 ENDIF 5299 5240 ENDIF 5300 5301 5241 ! 5302 5242 !-- Set Neumann boundary conditions for humidity and cloud-physical quantities
Note: See TracChangeset
for help on using the changeset viewer.