Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3655)
+++ palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3802)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! unused subroutines commented out
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! unused variables removed
 ! 
@@ -1513,97 +1516,97 @@
       
       
-       SUBROUTINE interpolate_to_fine_flux
-      
-      
-           USE arrays_3d
-           USE control_parameters
-           USE grid_variables
-           USE indices
-           USE pegrid
-           
-      
-           IMPLICIT NONE
-
-           INTEGER(iwp)                       :: i
-           INTEGER(iwp)                       :: j
-           INTEGER(iwp)                       :: iif
-           INTEGER(iwp)                       :: jjf
-           INTEGER(iwp)                       :: bottomx
-           INTEGER(iwp)                       :: bottomy
-           INTEGER(iwp)                       :: topx
-           INTEGER(iwp)                       :: topy
-           REAL(wp)                           :: eps
-           REAL(wp)                           :: alpha
-           REAL(wp)                           :: eminus
-           REAL(wp)                           :: edot
-           REAL(wp)                           :: eplus
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uswspr
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: vswspr
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: tspr 
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uspr
-      
-           ALLOCATE( uswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
-           ALLOCATE( vswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
-           ALLOCATE( tspr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
-           ALLOCATE( uspr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
-      
-           !
-           !-- Initialisation of scalar variables (2D)
-      
-           !
-           !-- Interpolation in x-direction (quadratic, Clark and Farley)
-      
-           DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
-               DO i = bdims_rem(1,1), bdims_rem(1,2)
-              
-                   bottomx = (nxf+1)/(nxc+1) * i
-                   topx    = (nxf+1)/(nxc+1) * (i+1) - 1
-      
-                   DO iif = bottomx, topx
-      
-                       eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
-                       alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
-                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
-                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-      
-                       uswspr(j,iif) = eminus * work2dusws(j,i-1) &
-                           + edot  * work2dusws(j,i)   &
-                           + eplus  * work2dusws(j,i+1)
-      
-                       vswspr(j,iif) = eminus * work2dvsws(j,i-1) &
-                           + edot  * work2dvsws(j,i)   &
-                           + eplus  * work2dvsws(j,i+1)
-      
-                       tspr(j,iif)   = eminus * work2dts(j,i-1) &
-                           + edot  * work2dts(j,i)   &
-                           + eplus  * work2dts(j,i+1)
-      
-                       uspr(j,iif)   = eminus * work2dus(j,i-1) &
-                           + edot  * work2dus(j,i)   &
-                           + eplus  * work2dus(j,i+1)
-      
-                   END DO
-      
-               END DO
-           END DO
-      
-           !
-           !-- Interpolation in y-direction (quadratic, Clark and Farley)
-      
-           DO j = bdims_rem(2,1), bdims_rem(2,2)
-             
-               bottomy = (nyf+1)/(nyc+1) * j
-               topy    = (nyf+1)/(nyc+1) * (j+1) - 1
-      
-               DO iif = nxl, nxr
-                   DO jjf = bottomy, topy
-      
-                       eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
-                       alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
-                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
-                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
- 
+!       SUBROUTINE interpolate_to_fine_flux
+!      
+!      
+!           USE arrays_3d
+!           USE control_parameters
+!           USE grid_variables
+!           USE indices
+!           USE pegrid
+!           
+!      
+!           IMPLICIT NONE
+!
+!           INTEGER(iwp)                       :: i
+!           INTEGER(iwp)                       :: j
+!           INTEGER(iwp)                       :: iif
+!           INTEGER(iwp)                       :: jjf
+!           INTEGER(iwp)                       :: bottomx
+!           INTEGER(iwp)                       :: bottomy
+!           INTEGER(iwp)                       :: topx
+!           INTEGER(iwp)                       :: topy
+!           REAL(wp)                           :: eps
+!           REAL(wp)                           :: alpha
+!           REAL(wp)                           :: eminus
+!           REAL(wp)                           :: edot
+!           REAL(wp)                           :: eplus
+!           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uswspr
+!           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: vswspr
+!           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: tspr 
+!           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uspr
+!      
+!           ALLOCATE( uswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+!           ALLOCATE( vswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+!           ALLOCATE( tspr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+!           ALLOCATE( uspr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+!      
+!           !
+!           !-- Initialisation of scalar variables (2D)
+!      
+!           !
+!           !-- Interpolation in x-direction (quadratic, Clark and Farley)
+!      
+!           DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+!               DO i = bdims_rem(1,1), bdims_rem(1,2)
+!              
+!                   bottomx = (nxf+1)/(nxc+1) * i
+!                   topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+!      
+!                   DO iif = bottomx, topx
+!      
+!                       eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+!                       alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+!                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+!                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+!      
+!                       uswspr(j,iif) = eminus * work2dusws(j,i-1) &
+!                           + edot  * work2dusws(j,i)   &
+!                           + eplus  * work2dusws(j,i+1)
+!      
+!                       vswspr(j,iif) = eminus * work2dvsws(j,i-1) &
+!                           + edot  * work2dvsws(j,i)   &
+!                           + eplus  * work2dvsws(j,i+1)
+!      
+!                       tspr(j,iif)   = eminus * work2dts(j,i-1) &
+!                           + edot  * work2dts(j,i)   &
+!                           + eplus  * work2dts(j,i+1)
+!      
+!                       uspr(j,iif)   = eminus * work2dus(j,i-1) &
+!                           + edot  * work2dus(j,i)   &
+!                           + eplus  * work2dus(j,i+1)
+!      
+!                   END DO
+!      
+!               END DO
+!           END DO
+!      
+!           !
+!           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+!      
+!           DO j = bdims_rem(2,1), bdims_rem(2,2)
+!             
+!               bottomy = (nyf+1)/(nyc+1) * j
+!               topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+!      
+!               DO iif = nxl, nxr
+!                   DO jjf = bottomy, topy
+!      
+!                       eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+!                       alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+!                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+!                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+! 
 !!
 !!--   TODO
@@ -1625,16 +1628,16 @@
 !                        + edot  * uspr(j,iif)   &
 !                        + eplus  * uspr(j+1,iif)
-      
-                   END DO
-               END DO
-      
-           END DO
-      
-      
-           DEALLOCATE( uswspr, vswspr )
-           DEALLOCATE( tspr, uspr )
-      
-      
-       END SUBROUTINE interpolate_to_fine_flux
+!      
+!                   END DO
+!               END DO
+!      
+!           END DO
+!      
+!      
+!           DEALLOCATE( uswspr, vswspr )
+!           DEALLOCATE( tspr, uspr )
+!      
+!      
+!       END SUBROUTINE interpolate_to_fine_flux
    
    
@@ -2398,279 +2401,279 @@
     
     
-    CONTAINS
- 
-       SUBROUTINE vnest_set_topbc_kh
-       
-       
-           USE arrays_3d
-           USE control_parameters
-           USE grid_variables
-           USE indices
-           USE pegrid
-           
-       
-           IMPLICIT NONE
-
-           INTEGER(iwp)                       :: i
-           INTEGER(iwp)                       :: j
-           INTEGER(iwp)                       :: k
-           INTEGER(iwp)                       :: iif
-           INTEGER(iwp)                       :: jjf
-           INTEGER(iwp)                       :: bottomx
-           INTEGER(iwp)                       :: bottomy
-           INTEGER(iwp)                       :: topx
-           INTEGER(iwp)                       :: topy
-           REAL(wp)                           :: eps
-           REAL(wp)                           :: alpha
-           REAL(wp)                           :: eminus
-           REAL(wp)                           :: edot
-           REAL(wp)                           :: eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
-       
-           
-       
-           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
-           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
-       
-           !
-           !-- Determination of a boundary condition for the potential temperature pt:
-           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
-       
-           !
-           !-- Interpolation in x-direction
-       
-           DO k = bdims_rem(3,1), bdims_rem(3,2)
-       
-               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
-       
-                   DO i = bdims_rem(1,1), bdims_rem(1,2)
-       
-                       bottomx = (nxf+1)/(nxc+1) * i
-                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
-       
-                       DO iif = bottomx, topx
-       
-                           eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
-       
-                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
-       
-                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
-       
-                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-       
-                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
-       
-                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
-                               + edot  * work3d(k,j,i)   &
-                               + eplus  * work3d(k,j,i+1)
-                       END DO
-       
-                   END DO
-       
-               END DO
-       
-           END DO
-       
-           !
-           !-- Interpolation in y-direction
-       
-           DO k = bdims_rem(3,1), bdims_rem(3,2)
-       
-               DO j = bdims_rem(2,1), bdims_rem(2,2)
-       
-                   bottomy = (nyf+1)/(nyc+1) * j
-                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
-       
-                   DO iif = nxl, nxr
-       
-                       DO jjf = bottomy, topy
-       
-                           eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
-       
-                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
-                   
-                           eminus = eps * (eps - 1.0) / 2.0 + alpha
-       
-                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-       
-                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-       
-                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
-                               + edot  * ptprf(k,j,iif)   &
-                               + eplus  * ptprf(k,j+1,iif)
-                       END DO
-       
-                   END DO
-       
-               END DO
-       
-           END DO
-       
-           !
-           !-- Interpolation in z-direction
-       
-           DO jjf = nys, nyn
-               DO iif = nxl, nxr
-       
-                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
-       
-                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
-       
-                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
-       
-                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-       
-                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-       
-                   kh (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif)   &
-                       + edot  * ptprs(bdims_rem(3,1)+1,jjf,iif) &
-                       + eplus  * ptprs(bdims_rem(3,1)+2,jjf,iif)
-       
-               END DO
-           END DO          
-           
-           DEALLOCATE ( ptprf, ptprs )
-       
-       
-       
-       END SUBROUTINE vnest_set_topbc_kh
-       
-       SUBROUTINE vnest_set_topbc_km
-       
-       
-           USE arrays_3d
-           USE control_parameters
-           USE grid_variables
-           USE indices
-           USE pegrid
-           
-       
-           IMPLICIT NONE
-
-           INTEGER(iwp)                       :: i
-           INTEGER(iwp)                       :: j
-           INTEGER(iwp)                       :: k
-           INTEGER(iwp)                       :: iif
-           INTEGER(iwp)                       :: jjf
-           INTEGER(iwp)                       :: bottomx
-           INTEGER(iwp)                       :: bottomy
-           INTEGER(iwp)                       :: topx
-           INTEGER(iwp)                       :: topy
-           REAL(wp)                           :: eps
-           REAL(wp)                           :: alpha
-           REAL(wp)                           :: eminus
-           REAL(wp)                           :: edot
-           REAL(wp)                           :: eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
-       
-           
-       
-           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
-           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
-       
-           !
-           !-- Determination of a boundary condition for the potential temperature pt:
-           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
-       
-           !
-           !-- Interpolation in x-direction
-       
-           DO k = bdims_rem(3,1), bdims_rem(3,2)
-       
-               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
-       
-                   DO i = bdims_rem(1,1), bdims_rem(1,2)
-       
-                       bottomx = (nxf+1)/(nxc+1) * i
-                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
-       
-                       DO iif = bottomx, topx
-       
-                           eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
-       
-                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
-       
-                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
-       
-                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-       
-                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
-       
-                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
-                               + edot  * work3d(k,j,i)   &
-                               + eplus  * work3d(k,j,i+1)
-                       END DO
-       
-                   END DO
-       
-               END DO
-       
-           END DO
-       
-           !
-           !-- Interpolation in y-direction
-       
-           DO k = bdims_rem(3,1), bdims_rem(3,2)
-       
-               DO j = bdims_rem(2,1), bdims_rem(2,2)
-       
-                   bottomy = (nyf+1)/(nyc+1) * j
-                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
-       
-                   DO iif = nxl, nxr
-       
-                       DO jjf = bottomy, topy
-       
-                           eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
-       
-                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
-                   
-                           eminus = eps * (eps - 1.0) / 2.0 + alpha
-       
-                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-       
-                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-       
-                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
-                               + edot  * ptprf(k,j,iif)   &
-                               + eplus  * ptprf(k,j+1,iif)
-                       END DO
-       
-                   END DO
-       
-               END DO
-       
-           END DO
-       
-           !
-           !-- Interpolation in z-direction
-       
-           DO jjf = nys, nyn
-               DO iif = nxl, nxr
-       
-                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
-       
-                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
-       
-                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
-       
-                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
-       
-                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-       
-                   km (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif)   &
-                       + edot  * ptprs(bdims_rem(3,1)+1,jjf,iif) &
-                       + eplus  * ptprs(bdims_rem(3,1)+2,jjf,iif)
-       
-               END DO
-           END DO          
-           
-           DEALLOCATE ( ptprf, ptprs )
-       
-       
-       
-       END SUBROUTINE vnest_set_topbc_km
+!    CONTAINS
+! 
+!       SUBROUTINE vnest_set_topbc_kh
+!       
+!       
+!           USE arrays_3d
+!           USE control_parameters
+!           USE grid_variables
+!           USE indices
+!           USE pegrid
+!           
+!       
+!           IMPLICIT NONE
+!
+!           INTEGER(iwp)                       :: i
+!           INTEGER(iwp)                       :: j
+!           INTEGER(iwp)                       :: k
+!           INTEGER(iwp)                       :: iif
+!           INTEGER(iwp)                       :: jjf
+!           INTEGER(iwp)                       :: bottomx
+!           INTEGER(iwp)                       :: bottomy
+!           INTEGER(iwp)                       :: topx
+!           INTEGER(iwp)                       :: topy
+!           REAL(wp)                           :: eps
+!           REAL(wp)                           :: alpha
+!           REAL(wp)                           :: eminus
+!           REAL(wp)                           :: edot
+!           REAL(wp)                           :: eplus
+!           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
+!           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
+!       
+!           
+!       
+!           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+!           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
+!       
+!           !
+!           !-- Determination of a boundary condition for the potential temperature pt:
+!           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
+!       
+!           !
+!           !-- Interpolation in x-direction
+!       
+!           DO k = bdims_rem(3,1), bdims_rem(3,2)
+!       
+!               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+!       
+!                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+!       
+!                       bottomx = (nxf+1)/(nxc+1) * i
+!                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
+!       
+!                       DO iif = bottomx, topx
+!       
+!                           eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+!       
+!                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+!       
+!                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
+!       
+!                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!       
+!                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
+!       
+!                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
+!                               + edot  * work3d(k,j,i)   &
+!                               + eplus  * work3d(k,j,i+1)
+!                       END DO
+!       
+!                   END DO
+!       
+!               END DO
+!       
+!           END DO
+!       
+!           !
+!           !-- Interpolation in y-direction
+!       
+!           DO k = bdims_rem(3,1), bdims_rem(3,2)
+!       
+!               DO j = bdims_rem(2,1), bdims_rem(2,2)
+!       
+!                   bottomy = (nyf+1)/(nyc+1) * j
+!                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+!       
+!                   DO iif = nxl, nxr
+!       
+!                       DO jjf = bottomy, topy
+!       
+!                           eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+!       
+!                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+!                   
+!                           eminus = eps * (eps - 1.0) / 2.0 + alpha
+!       
+!                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!       
+!                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+!       
+!                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+!                               + edot  * ptprf(k,j,iif)   &
+!                               + eplus  * ptprf(k,j+1,iif)
+!                       END DO
+!       
+!                   END DO
+!       
+!               END DO
+!       
+!           END DO
+!       
+!           !
+!           !-- Interpolation in z-direction
+!       
+!           DO jjf = nys, nyn
+!               DO iif = nxl, nxr
+!       
+!                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+!       
+!                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+!       
+!                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+!       
+!                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!       
+!                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+!       
+!                   kh (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif)   &
+!                       + edot  * ptprs(bdims_rem(3,1)+1,jjf,iif) &
+!                       + eplus  * ptprs(bdims_rem(3,1)+2,jjf,iif)
+!       
+!               END DO
+!           END DO          
+!           
+!           DEALLOCATE ( ptprf, ptprs )
+!       
+!       
+!       
+!       END SUBROUTINE vnest_set_topbc_kh
+       
+!       SUBROUTINE vnest_set_topbc_km
+!       
+!       
+!           USE arrays_3d
+!           USE control_parameters
+!           USE grid_variables
+!           USE indices
+!           USE pegrid
+!           
+!       
+!           IMPLICIT NONE
+!
+!           INTEGER(iwp)                       :: i
+!           INTEGER(iwp)                       :: j
+!           INTEGER(iwp)                       :: k
+!           INTEGER(iwp)                       :: iif
+!           INTEGER(iwp)                       :: jjf
+!           INTEGER(iwp)                       :: bottomx
+!           INTEGER(iwp)                       :: bottomy
+!           INTEGER(iwp)                       :: topx
+!           INTEGER(iwp)                       :: topy
+!           REAL(wp)                           :: eps
+!           REAL(wp)                           :: alpha
+!           REAL(wp)                           :: eminus
+!           REAL(wp)                           :: edot
+!           REAL(wp)                           :: eplus
+!           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
+!           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
+!       
+!           
+!       
+!           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+!           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
+!       
+!           !
+!           !-- Determination of a boundary condition for the potential temperature pt:
+!           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
+!       
+!           !
+!           !-- Interpolation in x-direction
+!       
+!           DO k = bdims_rem(3,1), bdims_rem(3,2)
+!       
+!               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+!       
+!                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+!       
+!                       bottomx = (nxf+1)/(nxc+1) * i
+!                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
+!       
+!                       DO iif = bottomx, topx
+!       
+!                           eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+!       
+!                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+!       
+!                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
+!       
+!                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!       
+!                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
+!       
+!                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
+!                               + edot  * work3d(k,j,i)   &
+!                               + eplus  * work3d(k,j,i+1)
+!                       END DO
+!       
+!                   END DO
+!       
+!               END DO
+!       
+!           END DO
+!       
+!           !
+!           !-- Interpolation in y-direction
+!       
+!           DO k = bdims_rem(3,1), bdims_rem(3,2)
+!       
+!               DO j = bdims_rem(2,1), bdims_rem(2,2)
+!       
+!                   bottomy = (nyf+1)/(nyc+1) * j
+!                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+!       
+!                   DO iif = nxl, nxr
+!       
+!                       DO jjf = bottomy, topy
+!       
+!                           eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+!       
+!                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+!                   
+!                           eminus = eps * (eps - 1.0) / 2.0 + alpha
+!       
+!                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!       
+!                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+!       
+!                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+!                               + edot  * ptprf(k,j,iif)   &
+!                               + eplus  * ptprf(k,j+1,iif)
+!                       END DO
+!       
+!                   END DO
+!       
+!               END DO
+!       
+!           END DO
+!       
+!           !
+!           !-- Interpolation in z-direction
+!       
+!           DO jjf = nys, nyn
+!               DO iif = nxl, nxr
+!       
+!                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+!       
+!                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+!       
+!                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+!       
+!                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+!       
+!                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+!       
+!                   km (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif)   &
+!                       + edot  * ptprs(bdims_rem(3,1)+1,jjf,iif) &
+!                       + eplus  * ptprs(bdims_rem(3,1)+2,jjf,iif)
+!       
+!               END DO
+!           END DO          
+!           
+!           DEALLOCATE ( ptprf, ptprs )
+!       
+!       
+!       
+!       END SUBROUTINE vnest_set_topbc_km
     
  
