Ignore:
Timestamp:
Oct 2, 2018 2:02:54 PM (3 years ago)
Author:
gronemeier
Message:

removed global array and call of function MAXLOC; updated include paths for IMUK config

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r3294 r3299  
    2525! -----------------
    2626! $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
    2733! changes concerning modularization of ocean option
    2834!
     
    13201326! -----------------------------------------------------------------------------!
    13211327!> 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)
    13221330!------------------------------------------------------------------------------!
    13231331 SUBROUTINE tcm_init_mixing_length
     
    13621370    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
    13631371
    1364     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_0_global !< wall_flags_0 of whole domain
    13651372    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_dummy    !< dummy array required for MPI_ALLREDUCE command
    13661373
     
    15101517
    15111518!
    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
    15411522             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.  &
    15431524                     k > k_max_topo )  &
    15441525                   k_max_topo = k
     
    16141595!--                   First, limit size of data copied to vicinity by the domain
    16151596!--                   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 )
    16211610
    16221611                      CALL copy_into_vicinity( k, j, i,           &
     
    16241613                                               -rad_j_s, rad_j_n, &
    16251614                                               -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
    16961624!
    16971625!--                   Search for walls only if there is any within vicinity
     
    17911719       ENDDO  !k loop
    17921720
    1793        DEALLOCATE( wall_flags_0_global )
    1794 
    17951721    ENDIF  !LES or RANS mode
    17961722
     
    18061732!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
    18071733!> closest to the origin (0/0) of 'array'.
    1808 !> @todo this part of PALM does not reproduce the same results for optimized
    1809 !>   and debug options for the compiler. This should be fixed
    18101734!------------------------------------------------------------------------------!
    18111735    REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
     
    18171741       INTEGER(iwp), INTENT(IN) :: pos_i     !< x position of the yz-plane 'array'
    18181742
     1743       INTEGER(iwp) :: a                     !< loop index
     1744       INTEGER(iwp) :: b                     !< loop index
    18191745       INTEGER(iwp) :: jj                    !< loop index
    18201746
     1747       INTEGER(KIND=1) :: maximum            !< maximum of array along z dimension
     1748
    18211749       INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
    18221750
     
    18251753!
    18261754!--    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
    18311766!
    18321767!--    Set distance to the default maximum value (=search radius)
     
    18881823       INTEGER(iwp) :: k   !< loop index
    18891824
    1890 
    18911825       DO  i = il, ir
    18921826          DO  j = js, jn
    18931827             DO  k = kb, kt
    18941828                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 ) )
    18961830             ENDDO
    18971831          ENDDO
     
    25362470       IF ( intermediate_timestep_count == 2 )  diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i)
    25372471       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)
    25392473
    25402474       dum_dif = tend(:,j,i) - dum_adv - dum_pro                                !> @todo remove later
Note: See TracChangeset for help on using the changeset viewer.