Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4248)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4249)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Several grid-line matching tests changed to a round-off-error tolerant form
+! in pmci_setup_parent, pmci_define_index_mapping and pmci_check_grid_matching.
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -261,4 +265,6 @@
     INTEGER(iwp), PARAMETER ::  interpolation_scheme_lrsn  = 2  !< Interpolation scheme to be used on lateral boundaries
     INTEGER(iwp), PARAMETER ::  interpolation_scheme_t     = 3  !< Interpolation scheme to be used on top boundary
+
+    REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp                 !< Relative tolerence for grid-line matching tests and comparisons
 !
 !-- Coupler setup
@@ -682,4 +688,8 @@
     REAL(wp) ::  yez                  !< Minimum separation in the y-direction required between the child and
                                       !< parent boundaries (south or north)
+    REAL(wp)     ::  tolex            !< Tolerance for grid-line matching in x-direction    
+    REAL(wp)     ::  toley            !< Tolerance for grid-line matching in y-direction    
+    REAL(wp)     ::  tolez            !< Tolerance for grid-line matching in z-direction    
+
     CHARACTER(LEN=32) ::  myname      !< String for variable name such as 'u'
 
@@ -700,4 +710,11 @@
     LOGICAL :: msib_north_in_m        !< Logical auxiliary parameter for the overlap test: true if the north border
                                       !< of the child msib is within the y-range of the child m
+
+
+!
+!-- Grid-line tolerances.
+    tolex = tolefac * dx
+    toley = tolefac * dy
+    tolez = tolefac * dz(1)    
 !
 !-- Initialize the current pmc parent.
@@ -738,5 +755,6 @@
 !--       transfer
           DO  kp = 1, nzt                 
-             IF ( zw(kp) > child_height )  THEN
+!             IF ( zw(kp) > child_height )  THEN
+             IF ( zw(kp) - child_height > tolez )  THEN                   
                 nz_child = kp
                 EXIT
@@ -772,6 +790,6 @@
              right_limit = upper_right_coord_x
              north_limit = upper_right_coord_y
-             IF ( ( child_coord_x(nx_child+1) /= right_limit ) .OR.                                 &
-                  ( child_coord_y(ny_child+1) /= north_limit ) )  THEN
+             IF ( ( child_coord_x(nx_child+1) /= right_limit ) .OR.                                 &  ! Change this IF test to a round-off-error tolerant form
+                  ( child_coord_y(ny_child+1) /= north_limit ) )  THEN                
                 nomatch = 1
              ENDIF
@@ -786,8 +804,8 @@
              south_limit = lower_left_coord_y + yez
              north_limit = upper_right_coord_y - yez
-             IF ( ( child_coord_x(0) < left_limit )  .OR.                                           &
-                  ( child_coord_x(nx_child+1) > right_limit )  .OR.                                 &
-                  ( child_coord_y(0) < south_limit )  .OR.                                          &
-                  ( child_coord_y(ny_child+1) > north_limit ) )  THEN
+             IF ( ( left_limit - child_coord_x(0) > tolex )  .OR.                                   &
+                  ( child_coord_x(nx_child+1) - right_limit > tolex )  .OR.                         &
+                  ( south_limit - child_coord_y(0) > toley )  .OR.                                  &
+                  ( child_coord_y(ny_child+1) - north_limit > toley ) )  THEN
                 nomatch = 1
              ENDIF
@@ -797,5 +815,5 @@
 !--       that the top ghost layer of the child grid does not exceed 
 !--       the parent domain top boundary.
-          IF ( child_height > zw(nzt) ) THEN
+          IF ( child_height > zw(nzt) ) THEN   ! Consider changing also this IF-test although it is not critical.
              nomatch = 1
           ENDIF
@@ -1350,5 +1368,5 @@
        xpl     = coord_x(nxl) - xexl
        DO  ip = 0, pg%nx
-          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl )  THEN
+          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl )  THEN  ! Consider changing xpl to xpl - tolex
              ipl = MAX( 0, ip )
              EXIT
@@ -1367,5 +1385,5 @@
        xpr  = coord_x(nxr+1) + xexr
        DO  ip = pg%nx, 0 , -1
-          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr )  THEN
+          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr )  THEN  ! Consider changing xpr to xpr + tolex
              ipr = MIN( pg%nx, MAX( ipl, ip ) )
              EXIT
@@ -1384,5 +1402,5 @@
        yps  = coord_y(nys) - yexs
        DO  jp = 0, pg%ny
-          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps )  THEN
+          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps )  THEN  ! Consider changing yps to yps - toley
              jps = MAX( 0, jp )
              EXIT
@@ -1401,5 +1419,5 @@
        ypn  = coord_y(nyn+1) + yexn
        DO  jp = pg%ny, 0 , -1
-          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn )  THEN
+          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn )  THEN  ! Consider changing ypn to ypn + toley
              jpn = MIN( pg%ny, MAX( jps, jp ) )
              EXIT
@@ -1509,5 +1527,14 @@
        INTEGER(iwp) ::  kstart    !<
        INTEGER(iwp) ::  kw        !< Child-grid index limited to kw <= nzt+1 for wall_flags_0
-     
+
+       REAL(wp)     ::  tolex     !< Tolerance for grid-line matching in x-direction    
+       REAL(wp)     ::  toley     !< Tolerance for grid-line matching in y-direction    
+       REAL(wp)     ::  tolez     !< Tolerance for grid-line matching in z-direction    
+
+!
+!--    Grid-line tolerances.
+       tolex = tolefac * dx
+       toley = tolefac * dy
+       tolez = tolefac * dz(1)
 !
 !--    Allocate child-grid work arrays for interpolation.
@@ -1582,5 +1609,5 @@
 !--       subroutines.
           i = istart
-          DO WHILE ( coord_x(i) < pg%coord_x(ii) .AND. i < nxrg )
+          DO WHILE ( pg%coord_x(ii) - coord_x(i) > tolex  .AND. i < nxrg )
              i = i + 1
           ENDDO
@@ -1600,12 +1627,10 @@
        DO  ii = ipla, ipra
           i = istart
-          DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.     &
-                      ( i < nxrg ) )
+          DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.  ( i < nxrg ) )
              i  = i + 1
           ENDDO
           iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
           ir = i
-          DO WHILE ( ( coord_x(ir) + 0.5_wp * dx <= pg%coord_x(ii) + pg%dx )  &
-                      .AND.  ( i < nxrg+1 ) )
+          DO WHILE ( ( coord_x(ir) + 0.5_wp * dx < pg%coord_x(ii) + pg%dx )  .AND.  ( i < nxrg+1 ) )
              i  = i + 1
              ir = MIN( i, nxrg )
@@ -1631,5 +1656,5 @@
 !--       subroutines.
           j = jstart
-          DO WHILE ( coord_y(j) < pg%coord_y(jj) .AND. j < nyng )
+          DO WHILE ( pg%coord_y(jj) - coord_y(j) > toley  .AND. j < nyng )
              j = j + 1
           ENDDO
@@ -1654,5 +1679,5 @@
           jflo(jj) = MIN( MAX( j, nysg ), nyng )
           jr = j
-          DO WHILE ( ( coord_y(jr) + 0.5_wp * dy <= pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
+          DO WHILE ( ( coord_y(jr) + 0.5_wp * dy < pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
              j  = j + 1
              jr = MIN( j, nyng )
@@ -1682,5 +1707,5 @@
 !--       subroutines.
           k = kstart
-          DO WHILE ( ( zw(k) < pg%zw(kk) )  .AND.  ( k < nzt+1 ) )
+          DO WHILE ( ( pg%zw(kk) - zw(k) > tolez )  .AND.  ( k < nzt+1 ) )
              k = k + 1
           ENDDO
@@ -1709,5 +1734,5 @@
           ENDDO
           kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-          DO WHILE ( ( zu(k) <= pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
+          DO WHILE ( ( zu(k) < pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
              k = k + 1
              IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
@@ -1970,16 +1995,18 @@
        INTEGER(iwp) ::  too_narrow_pesd_x = 0                !< Flag for too narrow pe-subdomain in x-direction
        INTEGER(iwp) ::  too_narrow_pesd_y = 0                !< Flag for too narrow pe-subdomain in y-direction
-                                                         
-       REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp           !< Relative tolerence for grid-line matching
-                                                         
+                                                                                                                  
        REAL(wp) ::  child_ngp_x_l                            !< Number of gridpoints in child subdomain in x-direction
                                                              !< converted to REAL(wp) 
        REAL(wp) ::  child_ngp_y_l                            !< Number of gridpoints in child subdomain in y-direction
                                                              !< converted to REAL(wp)
+       REAL(wp) ::  gridline_mismatch_x                      !< Mismatch between the parent and child gridlines in the x-direction
+       REAL(wp) ::  gridline_mismatch_y                      !< Mismatch between the parent and child gridlines in the y-direction
+       REAL(wp) ::  gsr_mismatch_x                           !< Deviation of the grid-spacing ratio from the nearest integer value, the x-direction
+       REAL(wp) ::  gsr_mismatch_y                           !< Deviation of the grid-spacing ratio from the nearest integer value, the y-direction
+       REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
+       REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
        REAL(wp) ::  tolex                                    !< Tolerance for grid-line matching in x-direction
        REAL(wp) ::  toley                                    !< Tolerance for grid-line matching in y-direction
        REAL(wp) ::  tolez                                    !< Tolerance for grid-line matching in z-direction
-       REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
-       REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
 
        
@@ -1991,12 +2018,16 @@
 !
 !--       First check that the child domain lower left corner matches the parent grid lines.
-          IF ( MOD( lower_left_coord_x, pg%dx ) > tolex ) non_matching_lower_left_corner = 1
-          IF ( MOD( lower_left_coord_y, pg%dy ) > toley ) non_matching_lower_left_corner = 1
+          gridline_mismatch_x = ABS( NINT( lower_left_coord_x / pg%dx ) * pg%dx - lower_left_coord_x )
+          gridline_mismatch_y = ABS( NINT( lower_left_coord_y / pg%dy ) * pg%dy - lower_left_coord_y )
+          IF ( gridline_mismatch_x > tolex ) non_matching_lower_left_corner = 1
+          IF ( gridline_mismatch_y > toley ) non_matching_lower_left_corner = 1
 !
 !--       Then check that the child doman upper right corner matches the parent grid lines.
           upper_right_coord_x = lower_left_coord_x + ( nx + 1 ) * dx
           upper_right_coord_y = lower_left_coord_y + ( ny + 1 ) * dy
-          IF ( MOD( upper_right_coord_x, pg%dx ) > tolex ) non_matching_upper_right_corner = 1
-          IF ( MOD( upper_right_coord_y, pg%dy ) > toley ) non_matching_upper_right_corner = 1
+          gridline_mismatch_x = ABS( NINT( upper_right_coord_x / pg%dx ) * pg%dx - upper_right_coord_x )
+          gridline_mismatch_y = ABS( NINT( upper_right_coord_y / pg%dy ) * pg%dy - upper_right_coord_y )
+          IF ( gridline_mismatch_x > tolex ) non_matching_upper_right_corner = 1
+          IF ( gridline_mismatch_y > toley ) non_matching_upper_right_corner = 1
 !
 !--       Also check that the cild domain height matches the parent grid lines.
@@ -2004,6 +2035,8 @@
 !
 !--       Check that the grid-spacing ratios in each direction are integer valued.    
-          IF ( MOD( pg%dx, dx ) > tolex )  non_int_gsr_x = 1
-          IF ( MOD( pg%dy, dy ) > toley )  non_int_gsr_y = 1
+          gsr_mismatch_x = ABS( NINT( pg%dx / dx ) * dx - pg%dx )
+          gsr_mismatch_y = ABS( NINT( pg%dy / dy ) * dy - pg%dy )
+          IF ( gsr_mismatch_x > tolex )  non_int_gsr_x = 1
+          IF ( gsr_mismatch_y > toley )  non_int_gsr_y = 1
 !
 !--       In the z-direction, all levels need to be checked separately against grid stretching  
@@ -2055,5 +2088,5 @@
 
           IF ( too_narrow_pesd_x > 0 )  THEN
-             WRITE( message_string, * ) 'child subdomain width in x-direction must not be ',        &
+            WRITE( message_string, * ) 'child subdomain width in x-direction must not be ',        &
                                         'smaller than its parent grid dx. Change the PE-grid ',     &
                                         'setting (npex, npey) to satisfy this requirement.'  
@@ -2185,6 +2218,6 @@
        ENDDO
 
-       write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
-       flush( 9 )
+!       write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
+!       flush( 9 )
        
     END SUBROUTINE pmci_compute_face_areas
@@ -5078,6 +5111,5 @@
     REAL(wp) ::  w_corr_top                   !< Correction added to the top boundary value of w 
 
-!    
-!-- The computation of the areas will be moved to initialization and made properly there.
+
     top_area = face_area(5)
 !
