Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2995)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2997)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Mask topography grid points in anterpolation
+! 
+! 2984 2018-04-18 11:51:30Z hellstea
 ! Bugfix in the log-law correction initialization. Pivot node cannot any more be 
 ! selected from outside the subdomain in order to prevent array under/overflows.
@@ -511,24 +514,23 @@
 !
 !-- Child-array indices to be precomputed and stored for anterpolation.
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !<
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !< child index indicating left bound of parent grid box on u-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !< child index indicating right bound of parent grid box on u-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !< child index indicating left bound of parent grid box on scalar-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !< child index indicating right bound of parent grid box on scalar-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !< child index indicating south bound of parent grid box on v-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !< child index indicating north bound of parent grid box on v-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !< child index indicating south bound of parent grid box on scalar-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !< child index indicating north bound of parent grid box on scalar-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !< child index indicating lower bound of parent grid box on w-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !< child index indicating upper bound of parent grid box on w-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !< child index indicating lower bound of parent grid box on scalar-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !< child index indicating upper bound of parent grid box on scalar-grid
 !
 !-- Number of fine-grid nodes inside coarse-grid ij-faces
 !-- to be precomputed for anterpolation.
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_u        !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_v        !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_s        !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_w         !<
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_s         !<
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_u  !< number of child grid boxes contribution to a parent grid box, u-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_v  !< number of child grid boxes contribution to a parent grid box, v-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_w  !< number of child grid boxes contribution to a parent grid box, w-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_s  !< number of child grid boxes contribution to a parent grid box, scalar-grid
     
     INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !<
@@ -2978,6 +2980,4 @@
 
        INTEGER(iwp) ::  i        !< Fine-grid index
-       INTEGER(iwp) ::  ifc_o    !<
-       INTEGER(iwp) ::  ifc_u    !<
        INTEGER(iwp) ::  ii       !< Coarse-grid index
        INTEGER(iwp) ::  istart   !<
@@ -3039,9 +3039,12 @@
        ALLOCATE( kfuo(0:kctu) )
 
-       ALLOCATE( ijfc_u(jcs:jcn,icl:icr) )
-       ALLOCATE( ijfc_v(jcs:jcn,icl:icr) )
-       ALLOCATE( ijfc_s(jcs:jcn,icl:icr) )
-       ALLOCATE( kfc_w(0:kctw) )
-       ALLOCATE( kfc_s(0:kctu) )
+       ALLOCATE( ijkfc_u(0:kctu,jcs:jcn,icl:icr) )
+       ALLOCATE( ijkfc_v(0:kctu,jcs:jcn,icl:icr) )
+       ALLOCATE( ijkfc_w(0:kctw,jcs:jcn,icl:icr) )
+       ALLOCATE( ijkfc_s(0:kctu,jcs:jcn,icl:icr) )
+       ijkfc_u = 0
+       ijkfc_v = 0
+       ijkfc_w = 0
+       ijkfc_s = 0
 !
 !--    i-indices of u for each ii-index value
@@ -3171,17 +3174,52 @@
 !--    This information is needed in anterpolation.
        DO  ii = icl, icr
-          ifc_u = ifuu(ii) - iflu(ii) + 1
-          ifc_o = ifuo(ii) - iflo(ii) + 1
           DO  jj = jcs, jcn
-             ijfc_u(jj,ii) = ifc_u * ( jfuo(jj) - jflo(jj) + 1 )
-             ijfc_v(jj,ii) = ifc_o * ( jfuv(jj) - jflv(jj) + 1 )
-             ijfc_s(jj,ii) = ifc_o * ( jfuo(jj) - jflo(jj) + 1 )             
+             DO kk = 0, kctu
+!
+!--             u-component
+                DO  i = iflu(ii), ifuu(ii)
+                   DO  j = jflo(jj), jfuo(jj)
+                      DO  k = kflo(kk), kfuo(kk)
+                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii) + MERGE( 1, 0,  &
+                                            BTEST( wall_flags_0(k,j,i), 1 ) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+!
+!--             v-component 
+                DO  i = iflo(ii), ifuo(ii)
+                   DO  j = jflv(jj), jfuv(jj)
+                      DO  k = kflo(kk), kfuo(kk)
+                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii) + MERGE( 1, 0,  &
+                                            BTEST( wall_flags_0(k,j,i), 2 ) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+!
+!--             scalars
+                DO  i = iflo(ii), ifuo(ii)
+                   DO  j = jflo(jj), jfuo(jj)
+                      DO  k = kflo(kk), kfuo(kk)
+                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii) + MERGE( 1, 0,  &
+                                            BTEST( wall_flags_0(k,j,i), 0 ) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+
+             DO kk = 0, kctw
+!
+!--             w-component 
+                DO  i = iflo(ii), ifuo(ii)
+                   DO  j = jflo(jj), jfuo(jj)
+                      DO  k = kflw(kk), kfuw(kk)
+                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii) + MERGE( 1, 0,  &
+                                            BTEST( wall_flags_0(k,j,i), 3 ) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+        
           ENDDO
-       ENDDO
-       DO kk = 0, kctw
-          kfc_w(kk) = kfuw(kk) - kflw(kk) + 1
-       ENDDO
-       DO kk = 0, kctu
-          kfc_s(kk) = kfuo(kk) - kflo(kk) + 1
        ENDDO
 !
@@ -5012,9 +5050,9 @@
 
       CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
-                               kfuo, ijfc_u, kfc_s, 'u' )
+                               kfuo, ijkfc_u, 'u' )
       CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
-                               kfuo, ijfc_v, kfc_s, 'v' )
+                               kfuo, ijkfc_v, 'v' )
       CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
-                               kfuw, ijfc_s, kfc_w, 'w' )
+                               kfuw, ijkfc_w, 'w' )
 !
 !--   Anterpolation of TKE and dissipation rate if parent and child are in 
@@ -5022,10 +5060,10 @@
       IF ( rans_mode_parent  .AND.  rans_mode )  THEN
          CALL pmci_anterp_tophat( e, ec, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
-                                  kfuo, ijfc_s, kfc_s, 'e' )
+                                  kfuo, ijkfc_s, 'e' )
 !
 !--      Anterpolation of dissipation rate only if TKE-e closure is applied.
          IF ( rans_tke_e )  THEN
             CALL pmci_anterp_tophat( diss, dissc, kctu, iflo, ifuo, jflo, jfuo,&
-                                     kflo, kfuo, ijfc_s, kfc_s, 'diss' )
+                                     kflo, kfuo, ijkfc_s, 'diss' )
          ENDIF
 
@@ -5034,5 +5072,5 @@
       IF ( .NOT. neutral )  THEN
          CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo, &
-                                  kfuo, ijfc_s, kfc_s, 'pt' )
+                                  kfuo, ijkfc_s, 'pt' )
       ENDIF
 
@@ -5040,13 +5078,13 @@
 
          CALL pmci_anterp_tophat( q, q_c, kctu, iflo, ifuo, jflo, jfuo, kflo,  &
-                                  kfuo, ijfc_s, kfc_s, 'q' )
+                                  kfuo, ijkfc_s, 'q' )
 
          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
 
             CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
-                                     kflo, kfuo, ijfc_s, kfc_s, 'qc' )
+                                     kflo, kfuo, ijkfc_s, 'qc' )
 
             CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
-                                     kflo, kfuo, ijfc_s, kfc_s,  'nc' )
+                                     kflo, kfuo, ijkfc_s, 'nc' )
 
          ENDIF
@@ -5055,8 +5093,8 @@
 
             CALL pmci_anterp_tophat( qr, qrc, kctu, iflo, ifuo, jflo, jfuo,    &
-                                     kflo, kfuo, ijfc_s, kfc_s, 'qr' )
+                                     kflo, kfuo, ijkfc_s, 'qr' )
 
             CALL pmci_anterp_tophat( nr, nrc, kctu, iflo, ifuo, jflo, jfuo,    &
-                                     kflo, kfuo, ijfc_s, kfc_s, 'nr' )
+                                     kflo, kfuo, ijkfc_s, 'nr' )
 
          ENDIF
@@ -5066,5 +5104,5 @@
       IF ( passive_scalar )  THEN
          CALL pmci_anterp_tophat( s, sc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
-                                  kfuo, ijfc_s, kfc_s, 's' )
+                                  kfuo, ijkfc_s, 's' )
       ENDIF
 
@@ -5074,5 +5112,5 @@
                                      chem_spec_c(:,:,:,n),                     &
                                      kctu, iflo, ifuo, jflo, jfuo, kflo,       &
-                                     kfuo, ijfc_s, kfc_s, 's' )
+                                     kfuo, ijkfc_s, 's' )
          ENDDO
       ENDIF
@@ -5625,5 +5663,5 @@
 
     SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu,   &
-                                   ijfc, kfc, var )
+                                   ijkfc, var )
 !
 !--    Anterpolation of internal-node values to be used as the parent-domain
@@ -5635,4 +5673,6 @@
 
        CHARACTER(LEN=*), INTENT(IN) ::  var   !< identifyer for treated variable
+
+       INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
 
        INTEGER(iwp) ::  i         !< Running index x-direction - fine-grid
@@ -5647,17 +5687,14 @@
        INTEGER(iwp) ::  kk        !< Running index z-direction - coarse grid
        INTEGER(iwp) ::  kcb = 0   !< Bottom boundary index for anterpolation along z
-       INTEGER(iwp) ::  m         !< Running index surface type
-       INTEGER(iwp) ::  nfc       !<
+       INTEGER(iwp) ::  var_flag  !< bit number used to flag topography on respective grid
 
        INTEGER(iwp), INTENT(IN) ::  kct   !< Top boundary index for anterpolation along z
 
-       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl         !< Indicates start index of child cells belonging to certain parent cell - x direction
-       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu         !< Indicates end index of child cells belonging to certain parent cell - x direction
-       INTEGER(iwp), DIMENSION(jcs:jcn,icl:icr), INTENT(IN) :: ijfc !<
-       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl         !< Indicates start index of child cells belonging to certain parent cell - y direction
-       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu         !< Indicates start index of child cells belonging to certain parent cell - y direction
-       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfc         !<
-       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl         !< Indicates start index of child cells belonging to certain parent cell - z direction
-       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu         !< Indicates start index of child cells belonging to certain parent cell - z direction
+       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
 
        REAL(wp) ::  cellsum   !< sum of respective child cells belonging to parent cell 
@@ -5719,5 +5756,17 @@
              jcnm = jcn - nhln - 1
           ENDIF
-       ENDIF     
+       ENDIF
+!
+!--    Set masking bit for topography flags 
+       IF ( var == 'u' )  THEN 
+          var_flag = 1 
+       ELSEIF ( var == 'v' )  THEN
+          var_flag = 2 
+       ELSEIF ( var == 'w' )  THEN
+          var_flag = 3
+       ELSE
+          var_flag = 0
+       ENDIF  
+
 !
 !--    Note that ii, jj, and kk are coarse-grid indices and i,j, and k 
@@ -5729,12 +5778,11 @@
 !--          terrain too
              DO  kk = kcb, kct - 1
-!
-!--             ijfc and kfc are precomputed in pmci_init_anterp_tophat
-                nfc =  ijfc(jj,ii) * kfc(kk) 
+
                 cellsum = 0.0_wp
                 DO  i = ifl(ii), ifu(ii)
                    DO  j = jfl(jj), jfu(jj)
                       DO  k = kfl(kk), kfu(kk)
-                         cellsum = cellsum + f(k,j,i)
+                         cellsum = cellsum + MERGE( f(k,j,i), 0.0_wp,          &
+                                        BTEST( wall_flags_0(k,j,i), var_flag ) )
                       ENDDO
                    ENDDO
@@ -5743,6 +5791,15 @@
 !--             Spatial under-relaxation.
                 fra  = frax(ii) * fray(jj) * fraz(kk)
-                fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +               &
-                               fra * cellsum / REAL( nfc, KIND = wp )
+!
+!--             In case all child grid points are inside topography, i.e. 
+!--             cellsum is zero, also parent solution would have zero values at
+!--             that point, which may cause problems in particular for the 
+!--             temperature. Therefore, in case cellsum is zero, keep the 
+!--             parent solution at this point.  
+                fc(kk,jj,ii) = MERGE( ( 1.0_wp - fra ) * fc(kk,jj,ii) +        &
+                                      fra * cellsum  /                         &
+                                      REAL( ijkfc(kk,jj,ii), KIND = wp ),      &
+                                      fc(kk,jj,ii),                            &
+                                      cellsum /= 0.0_wp )       
 
              ENDDO
