Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3987)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4010)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Mass (volume) flux correction included to ensure global mass conservation for child domains.
+! 
+! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
 ! 
@@ -648,7 +651,10 @@
     INTEGER(iwp) :: workarr_t_exchange_type_y
  
-    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !<
-
-    REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !<
+    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !< Array for communicating the parent-grid dimensions
+                                                                    !< to its children.
+
+    REAL(wp), DIMENSION(6)              ::  face_area               !< Surface area of each boundary face 
+    REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !< Array for communicating the real-type parent-grid 
+                                                                    !< parameters to its children.
 
     TYPE parentgrid_def
@@ -714,4 +720,12 @@
     END INTERFACE pmci_datatrans
 
+    INTERFACE pmci_ensure_nest_mass_conservation
+       MODULE PROCEDURE pmci_ensure_nest_mass_conservation
+    END INTERFACE pmci_ensure_nest_mass_conservation
+
+    INTERFACE pmci_ensure_nest_mass_conservation_vertical
+       MODULE PROCEDURE pmci_ensure_nest_mass_conservation_vertical
+    END INTERFACE pmci_ensure_nest_mass_conservation_vertical
+
     INTERFACE pmci_init
        MODULE PROCEDURE pmci_init
@@ -758,6 +772,7 @@
     PUBLIC pmci_set_swaplevel
     PUBLIC get_number_of_children, get_childid, get_child_edges, get_child_gridspacing
-
-
+    PUBLIC pmci_ensure_nest_mass_conservation
+    PUBLIC pmci_ensure_nest_mass_conservation_vertical
+    
  CONTAINS
 
@@ -1549,6 +1564,9 @@
 !
 !--    Check that the child and parent grid lines do match 
-       CALL pmci_check_grid_matching 
-
+       CALL pmci_check_grid_matching
+!       
+!--    Compute surface areas of the nest-boundary faces 
+       CALL pmci_compute_face_areas
+       
     ENDIF
 
@@ -2330,6 +2348,124 @@
        
     END SUBROUTINE pmci_check_grid_matching
-     
-#endif  
+
+
+
+    SUBROUTINE pmci_compute_face_areas
+
+       IMPLICIT NONE
+       REAL(wp)  :: face_area_local   !< Local (for the current pe) integral face area of the left boundary
+       REAL(wp)  :: sub_sum           !< Intermediate sum in order to improve the accuracy of the summation
+
+       INTEGER(iwp) :: i              !< Running index in the x-direction
+       INTEGER(iwp) :: j              !< Running index in the y-direction
+       INTEGER(iwp) :: k              !< Running index in the z-direction 
+       INTEGER(iwp) :: k_wall         !< Local topography top k-index
+       INTEGER(iwp) :: n              !< Running index over boundary faces
+
+       
+!
+!--    Sum up the volume flow through the left boundary
+       face_area(1) = 0.0_wp
+       face_area_local = 0.0_wp
+       IF ( bc_dirichlet_l )  THEN
+          i = 0
+          DO  j = nys, nyn
+             sub_sum = 0.0_wp
+             k_wall = get_topography_top_index_ji( j, i, 'u' )
+             DO   k = k_wall + 1, nzt
+                sub_sum = sub_sum + dzw(k)
+             ENDDO
+             face_area_local =  face_area_local + dy * sub_sum
+          ENDDO
+       ENDIF
+       
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( face_area_local, face_area(1), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+       face_area(1) = face_area_local
+#endif
+!
+!--    Sum up the volume flow through the right boundary
+       face_area(2) = 0.0_wp
+       face_area_local = 0.0_wp
+       IF ( bc_dirichlet_r )  THEN
+          i = nx
+          DO  j = nys, nyn
+             sub_sum = 0.0_wp
+             k_wall = get_topography_top_index_ji( j, i, 'u' )
+             DO   k = k_wall + 1, nzt
+                sub_sum = sub_sum + dzw(k)
+             ENDDO
+             face_area_local =  face_area_local + dy * sub_sum
+          ENDDO
+       ENDIF
+       
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( face_area_local, face_area(2), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+       face_area(2) = face_area_local
+#endif
+!
+!--    Sum up the volume flow through the south boundary
+       face_area(3) = 0.0_wp
+       face_area_local = 0.0_wp
+       IF ( bc_dirichlet_s )  THEN
+          j = 0
+          DO  i = nxl, nxr
+             sub_sum = 0.0_wp
+             k_wall = get_topography_top_index_ji( j, i, 'v' )
+             DO  k = k_wall + 1, nzt
+                sub_sum = sub_sum + dzw(k)
+             ENDDO
+             face_area_local = face_area_local + dx * sub_sum
+          ENDDO
+       ENDIF
+       
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( face_area_local, face_area(3), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+       face_area(3) = face_area_local
+#endif
+!
+!--    Sum up the volume flow through the north boundary
+       face_area(4) = 0.0_wp
+       face_area_local = 0.0_wp
+       IF ( bc_dirichlet_n )  THEN
+          j = ny
+          DO  i = nxl, nxr
+             sub_sum = 0.0_wp
+             k_wall = get_topography_top_index_ji( j, i, 'v' )
+             DO  k = k_wall + 1, nzt
+                sub_sum = sub_sum + dzw(k)
+             ENDDO
+             face_area_local = face_area_local + dx * sub_sum
+          ENDDO
+       ENDIF
+       
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( face_area_local, face_area(4), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+       face_area(4) = face_area_local
+#endif
+!
+!--    The top face area does not depend on the topography at all.       
+       face_area(5) = ( nx + 1 ) * ( ny + 1 ) * dx * dy
+!
+!--    The 6th element is used for the total area
+       face_area(6) = 0.0_wp
+       DO  n = 1, 5
+          face_area(6) = face_area(6) + face_area(n)
+       ENDDO
+
+       write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
+       flush( 9 )
+       
+    END SUBROUTINE pmci_compute_face_areas
+#endif
+    
  END SUBROUTINE pmci_setup_child
 
@@ -3932,5 +4068,4 @@
       REAL(wp) ::  c_interp_1   !< Value interpolated to the flux point in x direction from the parent-grid data
       REAL(wp) ::  c_interp_2   !< Auxiliary value interpolated  to the flux point in x direction from the parent-grid data
-
 ! 
 !--   Check which edge is to be handled
@@ -4154,6 +4289,4 @@
       INTEGER(iwp) ::  jcbc     !< Fixed child-grid index in the y-direction pointing to the boundary-value nodes
       INTEGER(iwp) ::  jcbgp    !< Index running over the redundant boundary ghost points in y-direction
-      INTEGER(iwp) ::  kc       !< Running child-grid index in the z-direction      
-      INTEGER(iwp) ::  kp       !< Running parent-grid index in the z-direction
       INTEGER(iwp) ::  jpbeg    !< Parent-grid index in the y-direction pointing to the starting point of workarr_sn
                                 !< in the parent-grid array
@@ -4162,8 +4295,11 @@
       INTEGER(iwp) ::  jpwp     !< Reduced parent-grid index in the y-direction for workarr_sn pointing to
                                 !< the first prognostic node
+      INTEGER(iwp) ::  kc       !< Running child-grid index in the z-direction      
+      INTEGER(iwp) ::  kp       !< Running parent-grid index in the z-direction
       REAL(wp) ::  cb           !< Interpolation coefficient for the boundary ghost node  
       REAL(wp) ::  cp           !< Interpolation coefficient for the first prognostic node
       REAL(wp) ::  c_interp_1   !< Value interpolated to the flux point in x direction from the parent-grid data
       REAL(wp) ::  c_interp_2   !< Auxiliary value interpolated  to the flux point in x direction from the parent-grid data
+
       
 !
@@ -4973,4 +5109,293 @@
 #endif
  END SUBROUTINE pmci_boundary_conds
-   
+
+
+ 
+ SUBROUTINE pmci_ensure_nest_mass_conservation
+
+!
+!-- Adjust the volume-flow rate through the nested boundaries so that the net volume
+!-- flow through all boundaries of the current nest domain becomes zero.
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i                        !< Running index in the x-direction
+    INTEGER(iwp) ::  ierr                     !< MPI error code
+    INTEGER(iwp) ::  j                        !< Running index in the y-direction
+    INTEGER(iwp) ::  k                        !< Running index in the z-direction
+    INTEGER(iwp) ::  n                        !< Running index over the boundary faces: l, r, s, n and t
+
+    REAL(wp) ::  dxdy                         !< Surface area of grid cell top face
+    REAL(wp) ::  innor                        !< Inner normal vector of the grid cell face
+    REAL(wp) ::  sub_sum                      !< Intermediate sum for reducing the loss of signifigant digits in 2-D summations
+    REAL(wp) ::  u_corr_left                  !< Correction added to the left boundary value of u
+    REAL(wp) ::  u_corr_right                 !< Correction added to the right boundary value of u
+    REAL(wp) ::  v_corr_south                 !< Correction added to the south boundary value of v
+    REAL(wp) ::  v_corr_north                 !< Correction added to the north boundary value of v
+    REAL(wp) ::  volume_flux_integral         !< Surface integral of volume flux over the domain boundaries
+    REAL(wp) ::  volume_flux_local            !< Surface integral of volume flux over the subdomain boundary face
+    REAL(wp) ::  w_corr_top                   !< Correction added to the top boundary value of w 
+
+    REAL(wp), DIMENSION(5) ::  volume_flux    !< Surface integral of volume flux over each boundary face of the domain
+
+    
+!
+!-- Sum up the volume flow through the left boundary
+    volume_flux(1) = 0.0_wp
+    volume_flux_local = 0.0_wp
+    IF ( bc_dirichlet_l )  THEN
+       i = 0
+       innor = dy
+       DO   j = nys, nyn
+          sub_sum = 0.0_wp
+          DO   k = nzb+1, nzt
+             sub_sum = sub_sum + innor * u(k,j,i) * dzw(k)                                          &
+                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1 ) )
+          ENDDO
+          volume_flux_local = volume_flux_local + sub_sum
+       ENDDO
+    ENDIF
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(1), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    volume_flux(1) = volume_flux_local
+#endif
+!
+!-- Sum up the volume flow through the right boundary
+    volume_flux(2) = 0.0_wp
+    volume_flux_local = 0.0_wp
+    IF ( bc_dirichlet_r )  THEN
+       i = nx + 1
+       innor = -dy
+       DO   j = nys, nyn
+          sub_sum = 0.0_wp
+          DO   k = nzb+1, nzt
+             sub_sum = sub_sum + innor * u(k,j,i) * dzw(k)                                          &
+                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1 ) )
+          ENDDO
+          volume_flux_local = volume_flux_local + sub_sum
+       ENDDO
+    ENDIF
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(2), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    volume_flux(2) = volume_flux_local
+#endif
+!
+!-- Sum up the volume flow through the south boundary
+    volume_flux(3) = 0.0_wp    
+    volume_flux_local = 0.0_wp
+    IF ( bc_dirichlet_s )  THEN
+       j = 0
+       innor = dx
+       DO   i = nxl, nxr
+          sub_sum = 0.0_wp
+          DO   k = nzb+1, nzt
+             sub_sum = sub_sum + innor * v(k,j,i) * dzw(k)                                          &
+                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) )
+          ENDDO
+          volume_flux_local = volume_flux_local + sub_sum
+       ENDDO
+    ENDIF
+    
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(3), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    volume_flux(3) = volume_flux_local
+#endif
+!    
+!-- Sum up the volume flow through the north boundary
+    volume_flux(4) = 0.0_wp
+    volume_flux_local = 0.0_wp
+    IF ( bc_dirichlet_n )  THEN
+       j = ny + 1
+       innor = -dx
+       DO  i = nxl, nxr
+          sub_sum = 0.0_wp
+          DO  k = nzb+1, nzt
+             sub_sum = sub_sum + innor * v(k,j,i) * dzw(k)                                          &
+                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) )
+          ENDDO
+          volume_flux_local = volume_flux_local + sub_sum
+       ENDDO
+    ENDIF
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(4), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    volume_flux(4) = volume_flux_local
+#endif
+!
+!-- Sum up the volume flow through the top boundary
+    volume_flux(5) = 0.0_wp
+    volume_flux_local = 0.0_wp
+    dxdy = dx * dy
+    k = nzt
+    DO  i = nxl, nxr
+       sub_sum = 0.0_wp
+       DO   j = nys, nyn
+          sub_sum = sub_sum - w(k,j,i) * dxdy  ! Minus, because the inner unit normal vector is (0,0,-1)
+       ENDDO
+       volume_flux_local = volume_flux_local + sub_sum
+    ENDDO
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(5), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    volume_flux(5) = volume_flux_local
+#endif
+
+    volume_flux_integral = 0.0_wp
+    DO  n = 1, 5
+       volume_flux_integral = volume_flux_integral + volume_flux(n)
+    ENDDO
+!    
+!-- Correction equally distributed to all nest boundaries, area_total must be used as area.
+!-- Note that face_area(6) is the total area (=sum from 1 to 5)
+    w_corr_top   = volume_flux_integral / face_area(6)
+    u_corr_left  =-w_corr_top
+    u_corr_right = w_corr_top
+    v_corr_south =-w_corr_top
+    v_corr_north = w_corr_top
+!!
+!!-- Just print out the net volume fluxes through each boundary. Only the root process prints.    
+!    if ( myid == 0 )  then       
+!       write( 9, "(5(e14.7,2x),4x,e14.7,4x,e12.5,4x,5(e14.7,2x))" )                                 &
+!            volume_flux(1), volume_flux(2), volume_flux(3), volume_flux(4), volume_flux(5),         &
+!            volume_flux_integral, c_correc,                                                         &
+!            u_corr_left, u_corr_right,  v_corr_south, v_corr_north, w_corr_top
+!       flush( 9 )
+!    endif    
+!
+!-- Correct the top-boundary value of w
+    DO   i = nxl, nxr
+       DO   j = nys, nyn
+          DO  k = nzt, nzt + 1
+             w(k,j,i) = w(k,j,i) + w_corr_top
+          ENDDO
+       ENDDO
+    ENDDO
+!
+!-- Correct the left-boundary value of u
+    IF ( bc_dirichlet_l )  THEN
+       DO  i = nxlg, nxl
+          DO  j = nys, nyn
+             DO  k = nzb + 1, nzt
+                u(k,j,i) = u(k,j,i) + u_corr_left
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+!
+!-- Correct the right-boundary value of u
+    IF ( bc_dirichlet_r )  THEN
+       DO  i = nxr+1, nxrg
+          DO  j = nys, nyn
+             DO  k = nzb + 1, nzt
+                u(k,j,i) = u(k,j,i) + u_corr_right
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+!
+!-- Correct the south-boundary value of v
+    IF ( bc_dirichlet_s )  THEN
+       DO  i = nxl, nxr
+          DO  j = nysg, nys
+             DO  k = nzb + 1, nzt
+                v(k,j,i) = v(k,j,i) + v_corr_south
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+!
+!-- Correct the north-boundary value of v
+    IF ( bc_dirichlet_n )  THEN
+       DO  i = nxl, nxr
+          DO  j = nyn+1, nyng
+             DO  k = nzb + 1, nzt
+                v(k,j,i) = v(k,j,i) + v_corr_north
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+    
+    
+ END SUBROUTINE pmci_ensure_nest_mass_conservation
+
+
+ 
+ SUBROUTINE pmci_ensure_nest_mass_conservation_vertical
+
+!
+!-- Adjust the volume-flow rate through the top boundary so that the net volume
+!-- flow through all boundaries of the current nest domain becomes zero.
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i                        !< Running index in the x-direction
+    INTEGER(iwp) ::  ierr                     !< MPI error code
+    INTEGER(iwp) ::  j                        !< Running index in the y-direction
+    INTEGER(iwp) ::  k                        !< Running index in the z-direction
+
+    REAL(wp) ::  dxdy                         !< Surface area of grid cell top face
+    REAL(wp) ::  innor                        !< Inner normal vector of the grid cell face
+    REAL(wp) ::  sub_sum                      !< Intermediate sum for reducing the loss of signifigant digits in 2-D summations
+    REAL(wp) ::  top_area                     !< Top boundary face area
+    REAL(wp) ::  volume_flux                  !< Surface integral of volume flux over the top boundary face
+    REAL(wp) ::  volume_flux_local            !< Surface integral of volume flux over the subdomain boundary face
+    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)
+!
+!-- Sum up the volume flow through the top boundary
+    volume_flux = 0.0_wp
+    volume_flux_local = 0.0_wp
+    dxdy = dx * dy
+    k = nzt
+    DO  i = nxl, nxr
+       sub_sum = 0.0_wp
+       DO   j = nys, nyn
+          sub_sum = sub_sum - w(k,j,i) * dxdy  ! Minus, because the inner unit normal vector is (0,0,-1)
+       ENDDO
+       volume_flux_local = volume_flux_local + sub_sum
+    ENDDO
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    volume_flux = volume_flux_local
+#endif
+
+    w_corr_top   = volume_flux / top_area
+!!
+!!-- Just print out the net volume fluxes through each boundary. Only the root process prints.    
+!    if ( myid == 0 )  then       
+!       write( 9, "(5(e14.7,2x),4x,e14.7,4x,e12.5,4x,5(e14.7,2x))" )                                 &
+!            volume_flux(1), volume_flux(2), volume_flux(3), volume_flux(4), volume_flux(5),         &
+!            volume_flux_integral, c_correc,                                                         &
+!            u_corr_left, u_corr_right,  v_corr_south, v_corr_north, w_corr_top
+!       flush( 9 )
+!    endif    
+!
+!-- Correct the top-boundary value of w
+    DO   i = nxl, nxr
+       DO   j = nys, nyn
+          DO  k = nzt, nzt + 1
+             w(k,j,i) = w(k,j,i) + w_corr_top
+          ENDDO
+       ENDDO
+    ENDDO
+    
+ END SUBROUTINE pmci_ensure_nest_mass_conservation_vertical
+
+ 
 END MODULE pmc_interface
