Index: palm/trunk/SOURCE/random_function_mod.f90
===================================================================
--- palm/trunk/SOURCE/random_function_mod.f90	(revision 3799)
+++ palm/trunk/SOURCE/random_function_mod.f90	(revision 3802)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! type conversion added to avoid compiler warning about constant integer
+! division truncation
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Corrected "Former revisions" section
 ! 
@@ -131,5 +135,5 @@
 
        PARAMETER ( ia=16807, im=2147483647, am=1.0_wp/im, iq=127773, ir=2836, &
-                   ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7_wp, rnmx=1.0_wp-eps )
+                   ntab=32, ndiv=1+INT(REAL(im-1)/ntab), eps=1.2e-7_wp, rnmx=1.0_wp-eps )
 
        IF ( idum .le. 0  .or.  random_iy .eq. 0 )  THEN
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3799)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3802)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! unused subroutine commented out
+! 
+! 3769 2019-02-28 10:16:49Z moh.hefny
 ! removed unused variables
 !
@@ -8629,99 +8632,99 @@
  
  
-     CONTAINS
- !------------------------------------------------------------------------------!
- ! Description:
- ! ------------
- !> Calculation of specific humidity of the skin layer (surface). It is assumend
- !> that the skin is always saturated.
- !------------------------------------------------------------------------------!
-        SUBROUTINE calc_q_surface_usm
- 
-           IMPLICIT NONE
- 
-           REAL(wp) :: resistance    !< aerodynamic and soil resistance term
- 
-           DO  m = 1, surf_usm_h%ns
- 
-              i   = surf_usm_h%i(m)            
-              j   = surf_usm_h%j(m)
-              k   = surf_usm_h%k(m)
- 
-!
-!--          Calculate water vapour pressure at saturation
-              e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                  &
-                                     ( t_surf_green_h_p(m) - 273.16_wp ) /  &
-                                     ( t_surf_green_h_p(m) - 35.86_wp  )    &
-                                          )
- 
-!
-!--          Calculate specific humidity at saturation
-              q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
- 
-!              surf_usm_h%r_a_green(m) = ( surf_usm_h%pt1(m) - t_surf_green_h(m) / exner(k) ) /  &
-!                    ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp )
-!                 
-! !--          make sure that the resistance does not drop to zero
-!              IF ( ABS(surf_usm_h%r_a_green(m)) < 1.0E-10_wp )  surf_usm_h%r_a_green(m) = 1.0E-10_wp
- 
-              resistance = surf_usm_h%r_a_green(m) / ( surf_usm_h%r_a_green(m) + surf_usm_h%r_s(m) + 1E-5_wp )
- 
-!
-!--          Calculate specific humidity at surface
-              IF ( bulk_cloud_model )  THEN
-                 q(k,j,i) = resistance * q_s +                   &
-                                            ( 1.0_wp - resistance ) *              &
-                                            ( q(k,j,i) - ql(k,j,i) )
-              ELSE
-                 q(k,j,i) = resistance * q_s +                   &
-                                            ( 1.0_wp - resistance ) *              &
-                                              q(k,j,i)
-              ENDIF
- 
-!
-!--          Update virtual potential temperature
-              vpt(k,j,i) = pt(k,j,i) *         &
-                         ( 1.0_wp + 0.61_wp * q(k,j,i) )
- 
-           ENDDO
-
-!
-!--       Now, treat vertical surface elements
-           DO  l = 0, 3
-              DO  m = 1, surf_usm_v(l)%ns
-!
-!--             Get indices of respective grid point
-                 i = surf_usm_v(l)%i(m)
-                 j = surf_usm_v(l)%j(m)
-                 k = surf_usm_v(l)%k(m)
- 
-!
-!--             Calculate water vapour pressure at saturation
-                 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                       &
-                                        ( t_surf_green_v_p(l)%t(m) - 273.16_wp ) /  &
-                                        ( t_surf_green_v_p(l)%t(m) - 35.86_wp  )    &
-                                             )
- 
-!
-!--             Calculate specific humidity at saturation
-                 q_s = 0.622_wp * e_s / ( surface_pressure -e_s )
- 
-!
-!--             Calculate specific humidity at surface
-                 IF ( bulk_cloud_model )  THEN
-                    q(k,j,i) = ( q(k,j,i) - ql(k,j,i) )
-                 ELSE
-                    q(k,j,i) = q(k,j,i)
-                 ENDIF
-!
-!--             Update virtual potential temperature
-                 vpt(k,j,i) = pt(k,j,i) *         &
-                            ( 1.0_wp + 0.61_wp * q(k,j,i) )
- 
-              ENDDO
- 
-           ENDDO
- 
-        END SUBROUTINE calc_q_surface_usm
+!     CONTAINS
+! !------------------------------------------------------------------------------!
+! ! Description:
+! ! ------------
+! !> Calculation of specific humidity of the skin layer (surface). It is assumend
+! !> that the skin is always saturated.
+! !------------------------------------------------------------------------------!
+!        SUBROUTINE calc_q_surface_usm
+! 
+!           IMPLICIT NONE
+! 
+!           REAL(wp) :: resistance    !< aerodynamic and soil resistance term
+! 
+!           DO  m = 1, surf_usm_h%ns
+! 
+!              i   = surf_usm_h%i(m)            
+!              j   = surf_usm_h%j(m)
+!              k   = surf_usm_h%k(m)
+! 
+!!
+!!--          Calculate water vapour pressure at saturation
+!              e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                  &
+!                                     ( t_surf_green_h_p(m) - 273.16_wp ) /  &
+!                                     ( t_surf_green_h_p(m) - 35.86_wp  )    &
+!                                          )
+! 
+!!
+!!--          Calculate specific humidity at saturation
+!              q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
+! 
+!!              surf_usm_h%r_a_green(m) = ( surf_usm_h%pt1(m) - t_surf_green_h(m) / exner(k) ) /  &
+!!                    ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp )
+!!                 
+!! !--          make sure that the resistance does not drop to zero
+!!              IF ( ABS(surf_usm_h%r_a_green(m)) < 1.0E-10_wp )  surf_usm_h%r_a_green(m) = 1.0E-10_wp
+! 
+!              resistance = surf_usm_h%r_a_green(m) / ( surf_usm_h%r_a_green(m) + surf_usm_h%r_s(m) + 1E-5_wp )
+! 
+!!
+!!--          Calculate specific humidity at surface
+!              IF ( bulk_cloud_model )  THEN
+!                 q(k,j,i) = resistance * q_s +                   &
+!                                            ( 1.0_wp - resistance ) *              &
+!                                            ( q(k,j,i) - ql(k,j,i) )
+!              ELSE
+!                 q(k,j,i) = resistance * q_s +                   &
+!                                            ( 1.0_wp - resistance ) *              &
+!                                              q(k,j,i)
+!              ENDIF
+! 
+!!
+!!--          Update virtual potential temperature
+!              vpt(k,j,i) = pt(k,j,i) *         &
+!                         ( 1.0_wp + 0.61_wp * q(k,j,i) )
+! 
+!           ENDDO
+!
+!!
+!!--       Now, treat vertical surface elements
+!           DO  l = 0, 3
+!              DO  m = 1, surf_usm_v(l)%ns
+!!
+!!--             Get indices of respective grid point
+!                 i = surf_usm_v(l)%i(m)
+!                 j = surf_usm_v(l)%j(m)
+!                 k = surf_usm_v(l)%k(m)
+! 
+!!
+!!--             Calculate water vapour pressure at saturation
+!                 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                       &
+!                                        ( t_surf_green_v_p(l)%t(m) - 273.16_wp ) /  &
+!                                        ( t_surf_green_v_p(l)%t(m) - 35.86_wp  )    &
+!                                             )
+! 
+!!
+!!--             Calculate specific humidity at saturation
+!                 q_s = 0.622_wp * e_s / ( surface_pressure -e_s )
+! 
+!!
+!!--             Calculate specific humidity at surface
+!                 IF ( bulk_cloud_model )  THEN
+!                    q(k,j,i) = ( q(k,j,i) - ql(k,j,i) )
+!                 ELSE
+!                    q(k,j,i) = q(k,j,i)
+!                 ENDIF
+!!
+!!--             Update virtual potential temperature
+!                 vpt(k,j,i) = pt(k,j,i) *         &
+!                            ( 1.0_wp + 0.61_wp * q(k,j,i) )
+! 
+!              ENDDO
+! 
+!           ENDDO
+! 
+!        END SUBROUTINE calc_q_surface_usm
         
      END SUBROUTINE usm_surface_energy_balance
Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3799)
+++ 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
     
  
