Changeset 3299 for palm/trunk/SOURCE
- Timestamp:
- Oct 2, 2018 2:02:54 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/turbulence_closure_mod.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/turbulence_closure_mod.f90
r3294 r3299 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall 28 ! calculation 29 ! - removed maxloc call as this produced different results for different 30 ! compiler options 31 ! 32 ! 3294 2018-10-01 02:37:10Z raasch 27 33 ! changes concerning modularization of ocean option 28 34 ! … … 1320 1326 ! -----------------------------------------------------------------------------! 1321 1327 !> Pre-computation of grid-dependent and near-wall mixing length. 1328 !> @todo consider walls in horizontal direction at a distance further than a 1329 !> single grid point (RANS mode) 1322 1330 !------------------------------------------------------------------------------! 1323 1331 SUBROUTINE tcm_init_mixing_length … … 1362 1370 INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k) 1363 1371 1364 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_0_global !< wall_flags_0 of whole domain1365 1372 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_dummy !< dummy array required for MPI_ALLREDUCE command 1366 1373 … … 1510 1517 1511 1518 ! 1512 !-- Gather topography information of whole domain 1513 !> @todo reduce amount of data sent by MPI call 1514 !> By now, a whole global 3D-array is sent and received with 1515 !> MPI_ALLREDUCE although most of the array is 0. This can be 1516 !> drastically reduced if only the local subarray is sent and stored 1517 !> in a global array. For that, an MPI data type or subarray must be 1518 !> defined. 1519 !> 2018-03-19, gronemeier 1520 ALLOCATE( wall_flags_0_global(nzb:nzt+1,0:ny,0:nx) ) 1521 1522 #if defined ( __parallel ) 1523 ALLOCATE( wall_flags_dummy(nzb:nzt+1,0:ny,0:nx) ) 1524 wall_flags_dummy = 0 1525 wall_flags_dummy(nzb:nzt+1,nys:nyn,nxl:nxr) = & 1526 wall_flags_0(nzb:nzt+1,nys:nyn,nxl:nxr) 1527 1528 CALL MPI_ALLREDUCE( wall_flags_dummy, & 1529 wall_flags_0_global, & 1530 (nzt-nzb+2)*(ny+1)*(nx+1), & 1531 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 1532 DEALLOCATE( wall_flags_dummy ) 1533 #else 1534 wall_flags_0_global(nzb:nzt+1,nys:nyn,nxl:nxr) = & 1535 wall_flags_0(nzb:nzt+1,nys:nyn,nxl:nxr) 1536 #endif 1537 ! 1538 !-- Get height level of highest topography 1539 DO i = 0, nx 1540 DO j = 0, ny 1519 !-- Get height level of highest topography within local subdomain 1520 DO i = nxlg, nxrg 1521 DO j = nysg, nyng 1541 1522 DO k = nzb+1, nzt-1 1542 IF ( .NOT. BTEST( wall_flags_0 _global(k,j,i), 0 ) .AND. &1523 IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND. & 1543 1524 k > k_max_topo ) & 1544 1525 k_max_topo = k … … 1614 1595 !-- First, limit size of data copied to vicinity by the domain 1615 1596 !-- border 1616 rad_i_l = MIN( rad_i, i ) 1617 rad_i_r = MIN( rad_i, nx-i ) 1618 1619 rad_j_s = MIN( rad_j, j ) 1620 rad_j_n = MIN( rad_j, ny-j ) 1597 !> @note limit copied area to 1 grid point in hor. dir. 1598 !> Ignore walls in horizontal direction which are 1599 !> further away than a single grid point. This allows to 1600 !> only search within local subdomain without the need 1601 !> of global topography information. 1602 !> The error made by this assumption are acceptable at 1603 !> the moment. 1604 !> 2018-10-01, gronemeier 1605 rad_i_l = MIN( 1, rad_i, i ) 1606 rad_i_r = MIN( 1, rad_i, nx-i ) 1607 1608 rad_j_s = MIN( 1, rad_j, j ) 1609 rad_j_n = MIN( 1, rad_j, ny-j ) 1621 1610 1622 1611 CALL copy_into_vicinity( k, j, i, & … … 1624 1613 -rad_j_s, rad_j_n, & 1625 1614 -rad_i_l, rad_i_r ) 1626 ! 1627 !-- In case of cyclic boundaries, copy parts into vicinity 1628 !-- where vicinity reaches over the domain borders. 1629 IF ( bc_lr_cyc ) THEN 1630 ! 1631 !-- Vicinity reaches over left domain boundary 1632 IF ( rad_i > rad_i_l ) THEN 1633 CALL copy_into_vicinity( k, j, nx+rad_i_l+1, & 1634 -rad_k_b, rad_k_t, & 1635 -rad_j_s, rad_j_n, & 1636 -rad_i, -rad_i_l-1 ) 1637 ! 1638 !-- ...and over southern domain boundary 1639 IF ( bc_ns_cyc .AND. rad_j > rad_j_s ) & 1640 CALL copy_into_vicinity( k, ny+rad_j_s+1, & 1641 nx+rad_i_l+1, & 1642 -rad_k_b, rad_k_t, & 1643 -rad_j, -rad_j_s-1, & 1644 -rad_i, -rad_i_l-1 ) 1645 ! 1646 !-- ...and over northern domain boundary 1647 IF ( bc_ns_cyc .AND. rad_j > rad_j_n ) & 1648 CALL copy_into_vicinity( k, 0-rad_j_n-1, & 1649 nx+rad_i_l+1, & 1650 -rad_k_b, rad_k_t, & 1651 rad_j_n+1, rad_j, & 1652 -rad_i, -rad_i_l-1 ) 1653 ENDIF 1654 ! 1655 !-- Vicinity reaches over right domain boundary 1656 IF ( rad_i > rad_i_r ) THEN 1657 CALL copy_into_vicinity( k, j, 0-rad_i_r-1, & 1658 -rad_k_b, rad_k_t, & 1659 -rad_j_s, rad_j_n, & 1660 rad_i_r+1, rad_i ) 1661 ! 1662 !-- ...and over southern domain boundary 1663 IF ( bc_ns_cyc .AND. rad_j > rad_j_s ) & 1664 CALL copy_into_vicinity( k, ny+rad_j_s+1, & 1665 0-rad_i_r-1, & 1666 -rad_k_b, rad_k_t, & 1667 -rad_j, -rad_j_s-1, & 1668 rad_i_r+1, rad_i ) 1669 ! 1670 !-- ...and over northern domain boundary 1671 IF ( bc_ns_cyc .AND. rad_j > rad_j_n ) & 1672 CALL copy_into_vicinity( k, 0-rad_j_n-1, & 1673 0-rad_i_r-1, & 1674 -rad_k_b, rad_k_t, & 1675 rad_j_n+1, rad_j, & 1676 rad_i_r+1, rad_i ) 1677 ENDIF 1678 ENDIF 1679 1680 IF ( bc_ns_cyc ) THEN 1681 ! 1682 !-- Vicinity reaches over southern domain boundary 1683 IF ( rad_j > rad_j_s ) & 1684 CALL copy_into_vicinity( k, ny+rad_j_s+1, i, & 1685 -rad_k_b, rad_k_t, & 1686 -rad_j, -rad_j_s-1, & 1687 -rad_i_l, rad_i_r ) 1688 ! 1689 !-- Vicinity reaches over northern domain boundary 1690 IF ( rad_j > rad_j_n ) & 1691 CALL copy_into_vicinity( k, 0-rad_j_n-1, i, & 1692 -rad_k_b, rad_k_t, & 1693 rad_j_n+1, rad_j, & 1694 rad_i_l, rad_i_r ) 1695 ENDIF 1615 !> @note in case of cyclic boundaries, those parts of the 1616 !> topography which lies beyond the domain borders but 1617 !> still within the search radius still needs to be 1618 !> copied into 'vicinity'. As the effective search 1619 !> radius is limited to 1 at the moment, no further 1620 !> copying is needed. Old implementation (prior to 1621 !> 2018-10-01) had this covered but used a global array. 1622 !> 2018-10-01, gronemeier 1623 1696 1624 ! 1697 1625 !-- Search for walls only if there is any within vicinity … … 1791 1719 ENDDO !k loop 1792 1720 1793 DEALLOCATE( wall_flags_0_global )1794 1795 1721 ENDIF !LES or RANS mode 1796 1722 … … 1806 1732 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array' 1807 1733 !> closest to the origin (0/0) of 'array'. 1808 !> @todo this part of PALM does not reproduce the same results for optimized1809 !> and debug options for the compiler. This should be fixed1810 1734 !------------------------------------------------------------------------------! 1811 1735 REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i ) … … 1817 1741 INTEGER(iwp), INTENT(IN) :: pos_i !< x position of the yz-plane 'array' 1818 1742 1743 INTEGER(iwp) :: a !< loop index 1744 INTEGER(iwp) :: b !< loop index 1819 1745 INTEGER(iwp) :: jj !< loop index 1820 1746 1747 INTEGER(KIND=1) :: maximum !< maximum of array along z dimension 1748 1821 1749 INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension 1822 1750 … … 1825 1753 ! 1826 1754 !-- Get coordinate of first maximum along vertical dimension 1827 !-- at each y position of array. 1828 !-- Substract 1 because indices count from 1 instead of 0 by MAXLOC 1829 loc_k = MAXLOC( array, DIM = 1) - 1 1830 1755 !-- at each y position of array (similar to function maxloc but more stable). 1756 DO a = 0, rad_j 1757 loc_k(a) = rad_k+1 1758 maximum = MAXVAL( array(:,a) ) 1759 DO b = 0, rad_k+1 1760 IF ( array(b,a) == maxi ) THEN 1761 loc_k(a) = b 1762 EXIT 1763 ENDIF 1764 ENDDO 1765 ENDDO 1831 1766 ! 1832 1767 !-- Set distance to the default maximum value (=search radius) … … 1888 1823 INTEGER(iwp) :: k !< loop index 1889 1824 1890 1891 1825 DO i = il, ir 1892 1826 DO j = js, jn 1893 1827 DO k = kb, kt 1894 1828 vicinity(k,j,i) = MERGE( 0, 1, & 1895 BTEST( wall_flags_0 _global(kp+k,jp+j,ip+i), 0 ) )1829 BTEST( wall_flags_0(kp+k,jp+j,ip+i), 0 ) ) 1896 1830 ENDDO 1897 1831 ENDDO … … 2536 2470 IF ( intermediate_timestep_count == 2 ) diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i) 2537 2471 IF ( intermediate_timestep_count == 3 ) diss_diff3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) - diss_prod3(:,j,i) 2538 IF ( intermediate_timestep_count == 3 ) dummy3(:,j,i) = km(:,j,i)2472 ! IF ( intermediate_timestep_count == 3 ) dummy3(:,j,i) = km(:,j,i) 2539 2473 2540 2474 dum_dif = tend(:,j,i) - dum_adv - dum_pro !> @todo remove later
Note: See TracChangeset
for help on using the changeset viewer.
