Index: /palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4259)
+++ /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4260)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Rest of the possibly round-off-error sensitive grid-line matching tests
+! changed to round-off-error tolerant forms throughout the module.
+! 
+! 4249 2019-10-01 12:27:47Z hellstea
 ! 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.
@@ -509,5 +513,5 @@
  SUBROUTINE pmci_init( world_comm )
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                         &
         ONLY:  message_string
 
@@ -521,5 +525,5 @@
 
 
-    CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,  &
+    CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,                       &
                          anterpolation_buffer_width, pmc_status )
 
@@ -536,7 +540,7 @@
 !
 !-- Check steering parameter values
-    IF ( TRIM( nesting_mode ) /= 'one-way'  .AND.                              &
-         TRIM( nesting_mode ) /= 'two-way'  .AND.                              &
-         TRIM( nesting_mode ) /= 'vertical' )                                  &
+    IF ( TRIM( nesting_mode ) /= 'one-way'  .AND.                                                   &
+         TRIM( nesting_mode ) /= 'two-way'  .AND.                                                   &
+         TRIM( nesting_mode ) /= 'vertical' )                                                       &
     THEN
        message_string = 'illegal nesting mode: ' // TRIM( nesting_mode )
@@ -544,10 +548,9 @@
     ENDIF
 
-    IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade'  .AND.                 &
-         TRIM( nesting_datatransfer_mode ) /= 'mixed'    .AND.                 &
-         TRIM( nesting_datatransfer_mode ) /= 'overlap' )                      &
+    IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade'  .AND.                                      &
+         TRIM( nesting_datatransfer_mode ) /= 'mixed'    .AND.                                      &
+         TRIM( nesting_datatransfer_mode ) /= 'overlap' )                                           &
     THEN
-       message_string = 'illegal nesting datatransfer mode: '                  &
-                        // TRIM( nesting_datatransfer_mode )
+       message_string = 'illegal nesting datatransfer mode: ' // TRIM( nesting_datatransfer_mode )
        CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 )
     ENDIF
@@ -558,8 +561,8 @@
 !-- Get some variables required by the pmc-interface (and in some cases in the
 !-- PALM code out of the pmci) out of the pmc-core
-    CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting,          &
-                             cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,   &
-                             cpl_name = cpl_name, npe_total = cpl_npe_total,   &
-                             lower_left_x = lower_left_coord_x,                &
+    CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting,                               &
+                             cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,                        &
+                             cpl_name = cpl_name, npe_total = cpl_npe_total,                        &
+                             lower_left_x = lower_left_coord_x,                                     &
                              lower_left_y = lower_left_coord_y )
 !
@@ -688,6 +691,6 @@
     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)     ::  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    
 
@@ -752,8 +755,6 @@
           child_height = child_grid_info(1)
 !
-!--       Find the highest child-domain level in the parent grid for the reduced z
-!--       transfer
+!--       Find the highest child-domain level in the parent grid for the reduced z transfer
           DO  kp = 1, nzt                 
-!             IF ( zw(kp) > child_height )  THEN
              IF ( zw(kp) - child_height > tolez )  THEN                   
                 nz_child = kp
@@ -790,6 +791,6 @@
              right_limit = upper_right_coord_x
              north_limit = upper_right_coord_y
-             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                
+             IF ( ( ABS( child_coord_x(nx_child+1) - right_limit ) > tolex )  .OR.                  &
+                  ( ABS( child_coord_y(ny_child+1) - north_limit ) > toley ) )  THEN                
                 nomatch = 1
              ENDIF
@@ -812,8 +813,7 @@
           ENDIF
 !             
-!--       Child domain must be lower than the parent domain such
-!--       that the top ghost layer of the child grid does not exceed 
-!--       the parent domain top boundary.
-          IF ( child_height > zw(nzt) ) THEN   ! Consider changing also this IF-test although it is not critical.
+!--       Child domain must be lower than the parent domain such that the top ghost
+!--       layer of the child grid does not exceed the parent domain top boundary.
+          IF ( child_height - zw(nzt) > tolez ) THEN
              nomatch = 1
           ENDIF
@@ -833,21 +833,21 @@
              DO  msib = 1, m - 1
 !
-!--             Set some logical auxiliary parameters to simplify the IF-condition.  
-                m_left_in_msib  = ( child_x_left(m)  >= child_x_left(msib) )  .AND.                 &
-                                  ( child_x_left(m)  <= child_x_right(msib) )
-                m_right_in_msib = ( child_x_right(m) >= child_x_left(msib) )  .AND.                 &
-                                  ( child_x_right(m) <= child_x_right(msib) )
-                msib_left_in_m  = ( child_x_left(msib)  >= child_x_left(m) )  .AND.                 &
-                                  ( child_x_left(msib)  <= child_x_right(m) )
-                msib_right_in_m = ( child_x_right(msib) >= child_x_left(m) )  .AND.                 &
-                                  ( child_x_right(msib) <= child_x_right(m) )
-                m_south_in_msib = ( child_y_south(m) >= child_y_south(msib) )  .AND.                &
-                                  ( child_y_south(m) <= child_y_north(msib) )
-                m_north_in_msib = ( child_y_north(m) >= child_y_south(msib) )  .AND.                &
-                                  ( child_y_north(m) <= child_y_north(msib) )
-                msib_south_in_m = ( child_y_south(msib) >= child_y_south(m) )  .AND.                &
-                                  ( child_y_south(msib) <= child_y_north(m) )
-                msib_north_in_m = ( child_y_north(msib) >= child_y_south(m) )  .AND.                &
-                                  ( child_y_north(msib) <= child_y_north(m) )
+!--             Set some logical auxiliary parameters to simplify the IF-condition.                  
+                m_left_in_msib  = ( child_x_left(m)  >= child_x_left(msib)  - tolex )  .AND.        &
+                                  ( child_x_left(m)  <= child_x_right(msib) + tolex )
+                m_right_in_msib = ( child_x_right(m) >= child_x_left(msib)  - tolex )  .AND.        &
+                                  ( child_x_right(m) <= child_x_right(msib) + tolex )
+                msib_left_in_m  = ( child_x_left(msib)  >= child_x_left(m)  - tolex )  .AND.        &
+                                  ( child_x_left(msib)  <= child_x_right(m) + tolex )
+                msib_right_in_m = ( child_x_right(msib) >= child_x_left(m)  - tolex )  .AND.        &
+                                  ( child_x_right(msib) <= child_x_right(m) + tolex )
+                m_south_in_msib = ( child_y_south(m) >= child_y_south(msib) - toley )  .AND.        &
+                                  ( child_y_south(m) <= child_y_north(msib) + toley )
+                m_north_in_msib = ( child_y_north(m) >= child_y_south(msib) - toley )  .AND.        &
+                                  ( child_y_north(m) <= child_y_north(msib) + toley )
+                msib_south_in_m = ( child_y_south(msib) >= child_y_south(m) - toley )  .AND.        &
+                                  ( child_y_south(msib) <= child_y_north(m) + toley )
+                msib_north_in_m = ( child_y_north(msib) >= child_y_south(m) - toley )  .AND.        &
+                                  ( child_y_north(msib) <= child_y_north(m) + toley )
                 
                 IF ( ( m_left_in_msib  .OR.  m_right_in_msib  .OR.                                  &
@@ -1333,4 +1333,6 @@
        INTEGER(iwp) ::  jauxn   !< Offset between the index bound jpn and the auxiliary index bound jpna
 
+       REAL(wp) ::  tolex       !< Tolerance for grid-line matching in x-direction    
+       REAL(wp) ::  toley       !< Tolerance for grid-line matching in y-direction    
        REAL(wp) ::  xexl        !< Parent-grid array exceedance behind the left edge of the child PE subdomain
        REAL(wp) ::  xexr        !< Parent-grid array exceedance behind the right edge of the child PE subdomain 
@@ -1355,5 +1357,9 @@
 !--    included in the neighbouring subdomain's parent-grid array, or not included at all if 
 !--    we are at the outer edge of the child domain. This may occur especially when a large 
-!--    grid-spacing ratio is used.       
+!--    grid-spacing ratio is used.
+!
+!--    Tolerances for grid-line matching.
+       tolex = tolefac * dx
+       toley = tolefac * dy
 !
 !--    Left boundary.
@@ -1368,5 +1374,5 @@
        xpl     = coord_x(nxl) - xexl
        DO  ip = 0, pg%nx
-          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl )  THEN  ! Consider changing xpl to xpl - tolex
+          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl - tolex )  THEN
              ipl = MAX( 0, ip )
              EXIT
@@ -1385,5 +1391,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  ! Consider changing xpr to xpr + tolex
+          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr + tolex )  THEN
              ipr = MIN( pg%nx, MAX( ipl, ip ) )
              EXIT
@@ -1402,5 +1408,5 @@
        yps  = coord_y(nys) - yexs
        DO  jp = 0, pg%ny
-          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps )  THEN  ! Consider changing yps to yps - toley
+          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps - toley )  THEN             
              jps = MAX( 0, jp )
              EXIT
@@ -1419,5 +1425,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  ! Consider changing ypn to ypn + toley
+          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn + toley )  THEN
              jpn = MIN( pg%ny, MAX( jps, jp ) )
              EXIT
@@ -1555,4 +1561,6 @@
 !--    First determine kcto and kctw which refer to the uppermost 
 !--    parent-grid levels below the child top-boundary level.
+!--    Note that these comparison tests are not round-off-error
+!--    sensitive and therefore tolerance buffering is not needed here.
        kk = 0
        DO WHILE ( pg%zu(kk) <= zu(nzt) )
@@ -1608,4 +1616,6 @@
 !--       are passed as arguments to the interpolation and anterpolation 
 !--       subroutines.
+!--       Note that this comparison test is round-off-error sensitive 
+!--       and therefore tolerance buffering is needed here.
           i = istart
           DO WHILE ( pg%coord_x(ii) - coord_x(i) > tolex  .AND. i < nxrg )
@@ -1623,5 +1633,7 @@
        WRITE( 9, * )
 !
-!--    i-indices of others for each ii-index value
+!--    i-indices of others for each ii-index value.
+!--    Note that these comparison tests are not round-off-error
+!--    sensitive and therefore tolerance buffering is not needed here.
        istart = nxlg
        DO  ii = ipla, ipra
@@ -1655,4 +1667,6 @@
 !--       are passed as arguments to the interpolation and anterpolation 
 !--       subroutines.
+!--       Note that this comparison test is round-off-error sensitive 
+!--       and therefore tolerance buffering is needed here.
           j = jstart
           DO WHILE ( pg%coord_y(jj) - coord_y(j) > toley  .AND. j < nyng )
@@ -1671,4 +1685,6 @@
 !
 !--    j-indices of others for each jj-index value
+!--    Note that these comparison tests are not round-off-error
+!--    sensitive and therefore tolerance buffering is not needed here.
        jstart = nysg
        DO  jj = jpsa, jpna
@@ -1706,4 +1722,6 @@
 !--       are passed as arguments to the interpolation and anterpolation 
 !--       subroutines.
+!--       Note that this comparison test is round-off-error sensitive 
+!--       and therefore tolerance buffering is needed here.
           k = kstart
           DO WHILE ( ( pg%zw(kk) - zw(k) > tolez )  .AND.  ( k < nzt+1 ) )
@@ -1728,4 +1746,6 @@
 !--    Note that anterpolation index limits are needed also for the top boundary 
 !--    ghost cell level because they are used also in the interpolation.
+!--    Note that these comparison tests are not round-off-error
+!--    sensitive and therefore tolerance buffering is not needed here.
        DO  kk = 1, pg%nz+1
           k = kstart
@@ -1741,12 +1761,4 @@
           kstart = kflo(kk)
        ENDDO
-!
-!--    Set the k-index bounds separately for the parent-grid cells pg%nz and pg%nz+1       
-!--    although they are not actually needed.
-!  WHY IS THIS LIKE THIS? REVISE WITH CARE.
-       kflo(pg%nz)   = nzt+1    
-       kfuo(pg%nz)   = nzt+kgsr 
-       kflo(pg%nz+1) = nzt+kgsr 
-       kfuo(pg%nz+1) = nzt+kgsr 
 !
 !--    Print out the index bounds for checking and debugging purposes
@@ -4521,5 +4533,5 @@
 !--   values. This subroutine is based on the first-order numerical
 !--   integration of the child-grid values contained within the anterpolation
-!--   cell.
+!--   cell (Clark & Farley, Journal of the Atmospheric Sciences 41(3), 1984).
 
       IMPLICIT NONE
