Ignore:
Timestamp:
Oct 1, 2019 12:27:47 PM (22 months ago)
Author:
hellstea
Message:

Several grid-line matching tests changed to a round-off-error tolerant form

File:
1 edited

Legend:

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

    r4182 r4249  
    2525! -----------------
    2626! $Id$
     27! Several grid-line matching tests changed to a round-off-error tolerant form
     28! in pmci_setup_parent, pmci_define_index_mapping and pmci_check_grid_matching.
     29!
     30! 4182 2019-08-22 15:20:23Z scharf
    2731! Corrected "Former revisions" section
    2832!
     
    261265    INTEGER(iwp), PARAMETER ::  interpolation_scheme_lrsn  = 2  !< Interpolation scheme to be used on lateral boundaries
    262266    INTEGER(iwp), PARAMETER ::  interpolation_scheme_t     = 3  !< Interpolation scheme to be used on top boundary
     267
     268    REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp                 !< Relative tolerence for grid-line matching tests and comparisons
    263269!
    264270!-- Coupler setup
     
    682688    REAL(wp) ::  yez                  !< Minimum separation in the y-direction required between the child and
    683689                                      !< parent boundaries (south or north)
     690    REAL(wp)     ::  tolex            !< Tolerance for grid-line matching in x-direction   
     691    REAL(wp)     ::  toley            !< Tolerance for grid-line matching in y-direction   
     692    REAL(wp)     ::  tolez            !< Tolerance for grid-line matching in z-direction   
     693
    684694    CHARACTER(LEN=32) ::  myname      !< String for variable name such as 'u'
    685695
     
    700710    LOGICAL :: msib_north_in_m        !< Logical auxiliary parameter for the overlap test: true if the north border
    701711                                      !< of the child msib is within the y-range of the child m
     712
     713
     714!
     715!-- Grid-line tolerances.
     716    tolex = tolefac * dx
     717    toley = tolefac * dy
     718    tolez = tolefac * dz(1)   
    702719!
    703720!-- Initialize the current pmc parent.
     
    738755!--       transfer
    739756          DO  kp = 1, nzt                 
    740              IF ( zw(kp) > child_height )  THEN
     757!             IF ( zw(kp) > child_height )  THEN
     758             IF ( zw(kp) - child_height > tolez )  THEN                   
    741759                nz_child = kp
    742760                EXIT
     
    772790             right_limit = upper_right_coord_x
    773791             north_limit = upper_right_coord_y
    774              IF ( ( child_coord_x(nx_child+1) /= right_limit ) .OR.                                 &
    775                   ( child_coord_y(ny_child+1) /= north_limit ) )  THEN
     792             IF ( ( child_coord_x(nx_child+1) /= right_limit ) .OR.                                 &  ! Change this IF test to a round-off-error tolerant form
     793                  ( child_coord_y(ny_child+1) /= north_limit ) )  THEN               
    776794                nomatch = 1
    777795             ENDIF
     
    786804             south_limit = lower_left_coord_y + yez
    787805             north_limit = upper_right_coord_y - yez
    788              IF ( ( child_coord_x(0) < left_limit )  .OR.                                           &
    789                   ( child_coord_x(nx_child+1) > right_limit )  .OR.                                 &
    790                   ( child_coord_y(0) < south_limit )  .OR.                                          &
    791                   ( child_coord_y(ny_child+1) > north_limit ) )  THEN
     806             IF ( ( left_limit - child_coord_x(0) > tolex )  .OR.                                   &
     807                  ( child_coord_x(nx_child+1) - right_limit > tolex )  .OR.                         &
     808                  ( south_limit - child_coord_y(0) > toley )  .OR.                                  &
     809                  ( child_coord_y(ny_child+1) - north_limit > toley ) )  THEN
    792810                nomatch = 1
    793811             ENDIF
     
    797815!--       that the top ghost layer of the child grid does not exceed
    798816!--       the parent domain top boundary.
    799           IF ( child_height > zw(nzt) ) THEN
     817          IF ( child_height > zw(nzt) ) THEN   ! Consider changing also this IF-test although it is not critical.
    800818             nomatch = 1
    801819          ENDIF
     
    13501368       xpl     = coord_x(nxl) - xexl
    13511369       DO  ip = 0, pg%nx
    1352           IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl )  THEN
     1370          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl )  THEN  ! Consider changing xpl to xpl - tolex
    13531371             ipl = MAX( 0, ip )
    13541372             EXIT
     
    13671385       xpr  = coord_x(nxr+1) + xexr
    13681386       DO  ip = pg%nx, 0 , -1
    1369           IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr )  THEN
     1387          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr )  THEN  ! Consider changing xpr to xpr + tolex
    13701388             ipr = MIN( pg%nx, MAX( ipl, ip ) )
    13711389             EXIT
     
    13841402       yps  = coord_y(nys) - yexs
    13851403       DO  jp = 0, pg%ny
    1386           IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps )  THEN
     1404          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps )  THEN  ! Consider changing yps to yps - toley
    13871405             jps = MAX( 0, jp )
    13881406             EXIT
     
    14011419       ypn  = coord_y(nyn+1) + yexn
    14021420       DO  jp = pg%ny, 0 , -1
    1403           IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn )  THEN
     1421          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn )  THEN  ! Consider changing ypn to ypn + toley
    14041422             jpn = MIN( pg%ny, MAX( jps, jp ) )
    14051423             EXIT
     
    15091527       INTEGER(iwp) ::  kstart    !<
    15101528       INTEGER(iwp) ::  kw        !< Child-grid index limited to kw <= nzt+1 for wall_flags_0
    1511      
     1529
     1530       REAL(wp)     ::  tolex     !< Tolerance for grid-line matching in x-direction   
     1531       REAL(wp)     ::  toley     !< Tolerance for grid-line matching in y-direction   
     1532       REAL(wp)     ::  tolez     !< Tolerance for grid-line matching in z-direction   
     1533
     1534!
     1535!--    Grid-line tolerances.
     1536       tolex = tolefac * dx
     1537       toley = tolefac * dy
     1538       tolez = tolefac * dz(1)
    15121539!
    15131540!--    Allocate child-grid work arrays for interpolation.
     
    15821609!--       subroutines.
    15831610          i = istart
    1584           DO WHILE ( coord_x(i) < pg%coord_x(ii) .AND. i < nxrg )
     1611          DO WHILE ( pg%coord_x(ii) - coord_x(i) > tolex .AND. i < nxrg )
    15851612             i = i + 1
    15861613          ENDDO
     
    16001627       DO  ii = ipla, ipra
    16011628          i = istart
    1602           DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.     &
    1603                       ( i < nxrg ) )
     1629          DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.  ( i < nxrg ) )
    16041630             i  = i + 1
    16051631          ENDDO
    16061632          iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
    16071633          ir = i
    1608           DO WHILE ( ( coord_x(ir) + 0.5_wp * dx <= pg%coord_x(ii) + pg%dx )  &
    1609                       .AND.  ( i < nxrg+1 ) )
     1634          DO WHILE ( ( coord_x(ir) + 0.5_wp * dx < pg%coord_x(ii) + pg%dx )  .AND.  ( i < nxrg+1 ) )
    16101635             i  = i + 1
    16111636             ir = MIN( i, nxrg )
     
    16311656!--       subroutines.
    16321657          j = jstart
    1633           DO WHILE ( coord_y(j) < pg%coord_y(jj) .AND. j < nyng )
     1658          DO WHILE ( pg%coord_y(jj) - coord_y(j) > toley .AND. j < nyng )
    16341659             j = j + 1
    16351660          ENDDO
     
    16541679          jflo(jj) = MIN( MAX( j, nysg ), nyng )
    16551680          jr = j
    1656           DO WHILE ( ( coord_y(jr) + 0.5_wp * dy <= pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
     1681          DO WHILE ( ( coord_y(jr) + 0.5_wp * dy < pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
    16571682             j  = j + 1
    16581683             jr = MIN( j, nyng )
     
    16821707!--       subroutines.
    16831708          k = kstart
    1684           DO WHILE ( ( zw(k) < pg%zw(kk) )  .AND.  ( k < nzt+1 ) )
     1709          DO WHILE ( ( pg%zw(kk) - zw(k) > tolez )  .AND.  ( k < nzt+1 ) )
    16851710             k = k + 1
    16861711          ENDDO
     
    17091734          ENDDO
    17101735          kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
    1711           DO WHILE ( ( zu(k) <= pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
     1736          DO WHILE ( ( zu(k) < pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
    17121737             k = k + 1
    17131738             IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
     
    19701995       INTEGER(iwp) ::  too_narrow_pesd_x = 0                !< Flag for too narrow pe-subdomain in x-direction
    19711996       INTEGER(iwp) ::  too_narrow_pesd_y = 0                !< Flag for too narrow pe-subdomain in y-direction
    1972                                                          
    1973        REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp           !< Relative tolerence for grid-line matching
    1974                                                          
     1997                                                                                                                 
    19751998       REAL(wp) ::  child_ngp_x_l                            !< Number of gridpoints in child subdomain in x-direction
    19761999                                                             !< converted to REAL(wp)
    19772000       REAL(wp) ::  child_ngp_y_l                            !< Number of gridpoints in child subdomain in y-direction
    19782001                                                             !< converted to REAL(wp)
     2002       REAL(wp) ::  gridline_mismatch_x                      !< Mismatch between the parent and child gridlines in the x-direction
     2003       REAL(wp) ::  gridline_mismatch_y                      !< Mismatch between the parent and child gridlines in the y-direction
     2004       REAL(wp) ::  gsr_mismatch_x                           !< Deviation of the grid-spacing ratio from the nearest integer value, the x-direction
     2005       REAL(wp) ::  gsr_mismatch_y                           !< Deviation of the grid-spacing ratio from the nearest integer value, the y-direction
     2006       REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
     2007       REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
    19792008       REAL(wp) ::  tolex                                    !< Tolerance for grid-line matching in x-direction
    19802009       REAL(wp) ::  toley                                    !< Tolerance for grid-line matching in y-direction
    19812010       REAL(wp) ::  tolez                                    !< Tolerance for grid-line matching in z-direction
    1982        REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
    1983        REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
    19842011
    19852012       
     
    19912018!
    19922019!--       First check that the child domain lower left corner matches the parent grid lines.
    1993           IF ( MOD( lower_left_coord_x, pg%dx ) > tolex ) non_matching_lower_left_corner = 1
    1994           IF ( MOD( lower_left_coord_y, pg%dy ) > toley ) non_matching_lower_left_corner = 1
     2020          gridline_mismatch_x = ABS( NINT( lower_left_coord_x / pg%dx ) * pg%dx - lower_left_coord_x )
     2021          gridline_mismatch_y = ABS( NINT( lower_left_coord_y / pg%dy ) * pg%dy - lower_left_coord_y )
     2022          IF ( gridline_mismatch_x > tolex ) non_matching_lower_left_corner = 1
     2023          IF ( gridline_mismatch_y > toley ) non_matching_lower_left_corner = 1
    19952024!
    19962025!--       Then check that the child doman upper right corner matches the parent grid lines.
    19972026          upper_right_coord_x = lower_left_coord_x + ( nx + 1 ) * dx
    19982027          upper_right_coord_y = lower_left_coord_y + ( ny + 1 ) * dy
    1999           IF ( MOD( upper_right_coord_x, pg%dx ) > tolex ) non_matching_upper_right_corner = 1
    2000           IF ( MOD( upper_right_coord_y, pg%dy ) > toley ) non_matching_upper_right_corner = 1
     2028          gridline_mismatch_x = ABS( NINT( upper_right_coord_x / pg%dx ) * pg%dx - upper_right_coord_x )
     2029          gridline_mismatch_y = ABS( NINT( upper_right_coord_y / pg%dy ) * pg%dy - upper_right_coord_y )
     2030          IF ( gridline_mismatch_x > tolex ) non_matching_upper_right_corner = 1
     2031          IF ( gridline_mismatch_y > toley ) non_matching_upper_right_corner = 1
    20012032!
    20022033!--       Also check that the cild domain height matches the parent grid lines.
     
    20042035!
    20052036!--       Check that the grid-spacing ratios in each direction are integer valued.   
    2006           IF ( MOD( pg%dx, dx ) > tolex )  non_int_gsr_x = 1
    2007           IF ( MOD( pg%dy, dy ) > toley )  non_int_gsr_y = 1
     2037          gsr_mismatch_x = ABS( NINT( pg%dx / dx ) * dx - pg%dx )
     2038          gsr_mismatch_y = ABS( NINT( pg%dy / dy ) * dy - pg%dy )
     2039          IF ( gsr_mismatch_x > tolex )  non_int_gsr_x = 1
     2040          IF ( gsr_mismatch_y > toley )  non_int_gsr_y = 1
    20082041!
    20092042!--       In the z-direction, all levels need to be checked separately against grid stretching 
     
    20552088
    20562089          IF ( too_narrow_pesd_x > 0 )  THEN
    2057              WRITE( message_string, * ) 'child subdomain width in x-direction must not be ',        &
     2090            WRITE( message_string, * ) 'child subdomain width in x-direction must not be ',        &
    20582091                                        'smaller than its parent grid dx. Change the PE-grid ',     &
    20592092                                        'setting (npex, npey) to satisfy this requirement.' 
     
    21852218       ENDDO
    21862219
    2187        write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
    2188        flush( 9 )
     2220!       write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
     2221!       flush( 9 )
    21892222       
    21902223    END SUBROUTINE pmci_compute_face_areas
     
    50785111    REAL(wp) ::  w_corr_top                   !< Correction added to the top boundary value of w
    50795112
    5080 !   
    5081 !-- The computation of the areas will be moved to initialization and made properly there.
     5113
    50825114    top_area = face_area(5)
    50835115!
Note: See TracChangeset for help on using the changeset viewer.