Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 2711)
+++ palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 2712)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Formatting and clean-up (SadiqHuq)
+! 
+! 2696 2017-12-14 17:12:51Z kanani
 ! renamed diffusivities to tcm_diffusivities (TG)
 ! 
@@ -48,4 +51,15 @@
 !>
 !> Definition of parameters and variables for vertical nesting
+!> The horizontal extent of the parent (Coarse Grid) and the child (Fine Grid)
+!> have to be identical. The vertical extent of the FG should be smaller than CG.
+!> Only integer nesting ratio supported. Odd nesting ratio preferred
+!> The code follows MPI-1 standards. The available processors are split into
+!> two groups using MPI_COMM_SPLIT. Exchange of data from CG to FG is called
+!> interpolation. FG initialization by interpolation is done once at the start.
+!> FG boundary conditions are set by interpolated at every timestep. 
+!> Exchange of data from CG to FG is called anterpolation, the two-way interaction
+!> occurs at every timestep.
+!> vnest_start_time set in PARIN allows delayed start of the coupling
+!> after spin-up of the CG
 !>
 !> @todo Ensure that code can be compiled for serial and parallel mode. Please
@@ -64,42 +78,70 @@
     IMPLICIT NONE
 
-    INTEGER(iwp),DIMENSION(3,2) ::  bdims = 0
-    INTEGER(iwp),DIMENSION(3,2) ::  bdims_rem = 0  !> Add description. It should not be longer than up to this point   |
-                                                   !> If really necessary, a second line can be added like this.
-    INTEGER(iwp)                             :: cg_nprocs, fg_nprocs
-    INTEGER(iwp),DIMENSION(:,:,:),ALLOCATABLE:: c2f_dims_cg, f2c_dims_cg
-    INTEGER(iwp),DIMENSION(:),ALLOCATABLE    :: c2f_dims_fg, f2c_dims_fg
-    INTEGER(iwp)                             :: TYPE_VNEST_BC, TYPE_VNEST_ANTER
-
-    INTEGER(iwp),DIMENSION(:,:),ALLOCATABLE  :: f_rnk_lst, c_rnk_lst
-    INTEGER(iwp),DIMENSION(3)                :: cfratio
-
-    INTEGER(iwp)                             :: nxc, nxf, nyc, nyf, nzc, nzf
-    INTEGER(iwp)                             :: ngp_c, ngp_f
-
-    INTEGER(iwp)                             :: target_idex, n_cell_c, n_cell_f
-    INTEGER(iwp),DIMENSION(2)                :: pdims_partner
-    INTEGER(iwp),DIMENSION(2)                :: offset,map_coord
-
-    REAL(wp)                                :: dxc, dyc, dxf, dyf, dzc, dzf,dtc,dtf
-
-    REAL(wp),DIMENSION(:,:,:), ALLOCATABLE  :: work3d 
-    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dshf
-    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dusws
-    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dvsws
-    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dts
-    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dus
-    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dz0
-
-
-    REAL(wp), DIMENSION(:),    ALLOCATABLE  :: zuc, zuf, zwc, zwf
-    REAL(wp), DIMENSION(:,:,:), POINTER      :: interpol3d,anterpol3d 
-!   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  interpol3d 
-
-    LOGICAL ::  vnest_init = .FALSE., vnested = .FALSE., &
-                vnest_twi = .FALSE., vnest_couple_rk3 = .FALSE.
-
-    ! PARIN
-    REAL(wp) ::  vnest_start_time = 9999999.9
+    LOGICAL                                   ::  vnested = .FALSE.            !> set to true when
+                                                                               !> mrun is called with -N option
+    LOGICAL                                   ::  vnest_init = .FALSE.         !> set to true when FG is initialized
+    REAL(wp)                                  ::  vnest_start_time = 9999999.9 !> simulated time when FG should be 
+                                                                               !> initialized. Should be
+                                                                               !> identical in PARIN & PARIN_N
+
+
+
+    INTEGER(iwp),DIMENSION(3,2)               :: bdims = 0        !> sub-domain grid topology of current PE
+    INTEGER(iwp),DIMENSION(3,2)               :: bdims_rem = 0    !> sub-domain grid topology of partner PE
+    INTEGER(iwp)                              :: cg_nprocs        !> no. of PE in CG. Set by mrun -N
+    INTEGER(iwp)                              :: fg_nprocs        !> no. of PE in FG. Set by mrun -N
+    INTEGER(iwp)                              :: TYPE_VNEST_BC    !> derived contiguous data type for interpolation
+    INTEGER(iwp)                              :: TYPE_VNEST_ANTER !> derived contiguous data type for anterpolation
+    INTEGER(iwp),DIMENSION(:,:,:),ALLOCATABLE :: c2f_dims_cg      !> One CG PE sends data to multiple FG PEs
+                                                                  !> list of grid-topology of partners
+    INTEGER(iwp),DIMENSION(:,:,:),ALLOCATABLE :: f2c_dims_cg      !> One CG PE receives data from multiple FG PEs
+                                                                  !> list of grid-topology of partners
+    INTEGER(iwp),DIMENSION(:),ALLOCATABLE     :: c2f_dims_fg      !> One FG PE sends data to multiple CG PE
+                                                                  !> list of grid-topology of partner
+    INTEGER(iwp),DIMENSION(:),ALLOCATABLE     :: f2c_dims_fg      !> One FG PE sends data to only one CG PE
+                                                                  !> list of grid-topology of partner
+
+    INTEGER(iwp),DIMENSION(:,:),ALLOCATABLE   :: f_rnk_lst        !> list storing rank of FG PE denoted by pdims
+    INTEGER(iwp),DIMENSION(:,:),ALLOCATABLE   :: c_rnk_lst        !> list storing rank of CG PE denoted by pdims
+    INTEGER(iwp),DIMENSION(3)                 :: cfratio          !> Nesting ratio in x,y and z-directions
+
+    INTEGER(iwp)                              :: nxc              !> no. of CG grid points in x-direction
+    INTEGER(iwp)                              :: nxf              !> no. of FG grid points in x-direction
+    INTEGER(iwp)                              :: nyc              !> no. of CG grid points in y-direction
+    INTEGER(iwp)                              :: nyf              !> no. of FG grid points in y-direction
+    INTEGER(iwp)                              :: nzc              !> no. of CG grid points in z-direction
+    INTEGER(iwp)                              :: nzf              !> no. of FG grid points in z-direction
+    INTEGER(iwp)                              :: ngp_c            !> no. of CG grid points in one vertical level
+    INTEGER(iwp)                              :: ngp_f            !> no. of FG grid points in one vertical level
+
+    INTEGER(iwp)                              :: n_cell_c         !> total no. of CG grid points in a PE
+    INTEGER(iwp)                              :: n_cell_f         !> total no. of FG grid points in a PE
+    INTEGER(iwp),DIMENSION(2)                 :: pdims_partner    !> processor topology of partner PE
+    INTEGER(iwp)                              :: target_idex      !> temporary variable
+    INTEGER(iwp),DIMENSION(2)                 :: offset           !> temporary variable
+    INTEGER(iwp),DIMENSION(2)                 :: map_coord        !> temporary variable
+
+    REAL(wp)                                  :: dxc              !> CG grid pacing in x-direction
+    REAL(wp)                                  :: dyc              !> FG grid pacing in x-direction
+    REAL(wp)                                  :: dxf              !> CG grid pacing in y-direction
+    REAL(wp)                                  :: dyf              !> FG grid pacing in y-direction
+    REAL(wp)                                  :: dzc              !> CG grid pacing in z-direction
+    REAL(wp)                                  :: dzf              !> FG grid pacing in z-direction
+    REAL(wp)                                  :: dtc              !> dt calculated for CG
+    REAL(wp)                                  :: dtf              !> dt calculated for FG
+
+    REAL(wp), DIMENSION(:),    ALLOCATABLE    :: zuc              !> CG vertical u-levels
+    REAL(wp), DIMENSION(:),    ALLOCATABLE    :: zuf              !> FG vertical u-levels
+    REAL(wp), DIMENSION(:),    ALLOCATABLE    :: zwc              !> CG vertical w-levels
+    REAL(wp), DIMENSION(:),    ALLOCATABLE    :: zwf              !> FG vertical w-levels
+    REAL(wp), DIMENSION(:,:,:), POINTER       :: interpol3d       !> pointers to simplify function calls
+    REAL(wp), DIMENSION(:,:,:), POINTER       :: anterpol3d       !> pointers to simplify function calls
+
+
+    REAL(wp),DIMENSION(:,:,:), ALLOCATABLE    :: work3d           !> temporary array for exchange of 3D data
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE    :: work2dusws       !> temporary array for exchange of 2D data
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE    :: work2dvsws       !> temporary array for exchange of 2D data
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE    :: work2dts         !> temporary array for exchange of 2D data
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE    :: work2dus         !> temporary array for exchange of 2D data
 
     SAVE
@@ -112,10 +154,8 @@
 
 !-- Public constants and variables
-    PUBLIC vnested, vnest_init, vnest_twi, vnest_couple_rk3,                  &
-           vnest_start_time
+    PUBLIC vnested, vnest_init, vnest_start_time
 
     PRIVATE bdims, bdims_rem,                                                 &
-            work3d, work2dshf, work2dusws, work2dvsws,                        &
-            work2dts, work2dus, work2dz0,                                     &
+            work3d, work2dusws, work2dvsws, work2dts, work2dus,               &
             dxc, dyc, dxf, dyf, dzc, dzf, dtc, dtf,                           &
             zuc, zuf, zwc, zwf, interpol3d, anterpol3d,                       &
@@ -176,4 +216,5 @@
     
     SUBROUTINE vnest_init_fine
+#if defined( __parallel )
    
         !--------------------------------------------------------------------------------!
@@ -191,19 +232,25 @@
         USE interfaces
         USE pegrid
-        USE surface_mod,                                                        &
+        USE surface_mod,                                                       &
             ONLY :  surf_def_h, surf_def_v
-        USE turbulence_closure_mod,                                             &
+        USE turbulence_closure_mod,                                            &
             ONLY :  tcm_diffusivities
+        
    
         IMPLICIT NONE
    
-        REAL(wp)    ::  time_since_reference_point_rem
-   
-        INTEGER(iwp)   :: i, j, k,im,jn,ko
-        INTEGER(iwp)                             ::  if, jf, kf
-   
-#if defined( __parallel )
-   
-        if (myid ==0 )print *, ' TIME TO INIT FINE from COARSE', simulated_time
+        REAL(wp)                              :: time_since_reference_point_rem
+        INTEGER(iwp)                          :: i
+        INTEGER(iwp)                          :: j
+        INTEGER(iwp)                          :: k
+        INTEGER(iwp)                          :: im
+        INTEGER(iwp)                          :: jn
+        INTEGER(iwp)                          :: ko
+        INTEGER(iwp)                          :: iif
+        INTEGER(iwp)                          :: jjf
+        INTEGER(iwp)                          :: kkf
+   
+   
+        if (myid ==0 ) print *, ' TIME TO INIT FINE from COARSE', simulated_time
    
         !
@@ -344,37 +391,30 @@
                         (bdims(2,2)-bdims(2,1)+3)
    
-!-- WARNING
-!-- shf,z0 not interpolated
-!-- line commented in interpolate_to_fine_flux
-!-- FG needs to read it's own data file   
-!MERGE-WIP                  CALL MPI_SEND(surf_def_h(0)%shf ( bdims(2,1)-1:bdims(2,2)+1, &
-!MERGE-WIP                      bdims(1,1)-1:bdims(1,2)+1),&
-!MERGE-WIP                      n_cell_c, MPI_REAL, target_idex,  &
-!MERGE-WIP                      109, comm_inter, ierr   )
-!MERGE-WIP 
-!MERGE-WIP                  CALL MPI_SEND(surf_def_h(0)%usws( bdims(2,1)-1:bdims(2,2)+1, &
-!MERGE-WIP                      bdims(1,1)-1:bdims(1,2)+1),&
-!MERGE-WIP                      n_cell_c, MPI_REAL, target_idex,  &
-!MERGE-WIP                      110, comm_inter, ierr   )
-!MERGE-WIP 
-!MERGE-WIP                  CALL MPI_SEND(surf_def_h(0)%vsws( bdims(2,1)-1:bdims(2,2)+1, &
-!MERGE-WIP                      bdims(1,1)-1:bdims(1,2)+1),&
-!MERGE-WIP                      n_cell_c, MPI_REAL, target_idex,  &
-!MERGE-WIP                      111, comm_inter, ierr   )
-!MERGE-WIP 
-!MERGE                   CALL MPI_SEND(ts  ( bdims(2,1)-1:bdims(2,2)+1, &
-!MERGE                       bdims(1,1)-1:bdims(1,2)+1),&
-!MERGE                       n_cell_c, MPI_REAL, target_idex,  &
-!MERGE                       112, comm_inter, ierr   )
-!MERGE  
-!MERGE                   CALL MPI_SEND(us  ( bdims(2,1)-1:bdims(2,2)+1, &
-!MERGE                       bdims(1,1)-1:bdims(1,2)+1),&
-!MERGE                       n_cell_c, MPI_REAL, target_idex,  &
-!MERGE                       113, comm_inter, ierr   )
-!MERGE  
-!MERGE                   CALL MPI_SEND(z0  ( bdims(2,1)-1:bdims(2,2)+1, &
-!MERGE                       bdims(1,1)-1:bdims(1,2)+1),&
-!MERGE                       n_cell_c, MPI_REAL, target_idex,  &
-!MERGE                       114, comm_inter, ierr   )
+!
+!--     shf and z0 for CG / FG need to initialized in input file or user_code
+!--     TODO
+!--     initialization of usws, vsws, ts and us not vital to vnest FG
+!--     variables are not compatible with the new surface layer module
+!
+!                 CALL MPI_SEND(surf_def_h(0)%usws( bdims(2,1)-1:bdims(2,2)+1, &
+!                     bdims(1,1)-1:bdims(1,2)+1),&
+!                     n_cell_c, MPI_REAL, target_idex,  &
+!                     110, comm_inter, ierr   )
+!
+!                 CALL MPI_SEND(surf_def_h(0)%vsws( bdims(2,1)-1:bdims(2,2)+1, &
+!                     bdims(1,1)-1:bdims(1,2)+1),&
+!                     n_cell_c, MPI_REAL, target_idex,  &
+!                     111, comm_inter, ierr   )
+!
+!                 CALL MPI_SEND(ts  ( bdims(2,1)-1:bdims(2,2)+1, &
+!                     bdims(1,1)-1:bdims(1,2)+1),&
+!                     n_cell_c, MPI_REAL, target_idex,  &
+!                     112, comm_inter, ierr   )
+!   
+!                 CALL MPI_SEND(us  ( bdims(2,1)-1:bdims(2,2)+1, &
+!                     bdims(1,1)-1:bdims(1,2)+1),&
+!                     n_cell_c, MPI_REAL, target_idex,  &
+!                     113, comm_inter, ierr   )
+!   
             ENDIF
    
@@ -465,6 +505,4 @@
                 (bdims_rem(2,2)-bdims_rem(2,1)+3)
    
-            ALLOCATE( work2dshf  ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
-                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
             ALLOCATE( work2dusws ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
                 bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
@@ -475,42 +513,38 @@
             ALLOCATE( work2dus   ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
                 bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
-            ALLOCATE( work2dz0   ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
-                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
-   
-!MERGE-WIP          CALL MPI_RECV( work2dshf ,n_cell_c, MPI_REAL, target_idex, 109, &
-!MERGE-WIP              comm_inter,status, ierr )
-!MERGE-WIP 
-!MERGE-WIP          CALL MPI_RECV( work2dusws,n_cell_c, MPI_REAL, target_idex, 110, &
-!MERGE-WIP              comm_inter,status, ierr )
-!MERGE-WIP 
-!MERGE-WIP          CALL MPI_RECV( work2dvsws,n_cell_c, MPI_REAL, target_idex, 111, &
-!MERGE-WIP              comm_inter,status, ierr )
-!MERGE-WIP 
-!MERGE           CALL MPI_RECV( work2dts  ,n_cell_c, MPI_REAL, target_idex, 112, &
-!MERGE               comm_inter,status, ierr )
-!MERGE  
-!MERGE           CALL MPI_RECV( work2dus  ,n_cell_c, MPI_REAL, target_idex, 113, &
-!MERGE               comm_inter,status, ierr )
-!MERGE  
-!MERGE           CALL MPI_RECV( work2dz0  ,n_cell_c, MPI_REAL, target_idex, 114, &
-!MERGE               comm_inter,status, ierr )
-!MERGE  
-!MERGE-WIP           CALL interpolate_to_fine_flux ( 108 )
-   
-            DEALLOCATE(   work2dshf    )
+   
+!
+!--     shf and z0 for CG / FG need to initialized in input file or user_code
+!--     TODO
+!--     initialization of usws, vsws, ts and us not vital to vnest FG
+!--     variables are not compatible with the new surface layer module
+! 
+!          CALL MPI_RECV( work2dusws,n_cell_c, MPI_REAL, target_idex, 110, &
+!              comm_inter,status, ierr )
+! 
+!          CALL MPI_RECV( work2dvsws,n_cell_c, MPI_REAL, target_idex, 111, &
+!              comm_inter,status, ierr )
+! 
+!          CALL MPI_RECV( work2dts  ,n_cell_c, MPI_REAL, target_idex, 112, &
+!              comm_inter,status, ierr )
+!   
+!          CALL MPI_RECV( work2dus  ,n_cell_c, MPI_REAL, target_idex, 113, &
+!              comm_inter,status, ierr )
+!   
+!           CALL interpolate_to_fine_flux ( 108 )
+   
             DEALLOCATE(   work2dusws   )
             DEALLOCATE(   work2dvsws   )
             DEALLOCATE(   work2dts     )
             DEALLOCATE(   work2dus     )
-            DEALLOCATE(   work2dz0     )
           ENDIF
    
           IF ( .NOT. constant_diffusion ) THEN
-             DO kf =  bdims(3,1)+1,bdims(3,2)+1
-                 DO jf =  bdims(2,1),bdims(2,2)
-                     DO if =  bdims(1,1),bdims(1,2)
-   
-                         IF ( e(kf,jf,if) < 0.0 ) THEN
-                              e(kf,jf,if) = 1E-15_wp
+             DO kkf =  bdims(3,1)+1,bdims(3,2)+1
+                 DO jjf =  bdims(2,1),bdims(2,2)
+                     DO iif =  bdims(1,1),bdims(1,2)
+   
+                         IF ( e(kkf,jjf,iif) < 0.0 ) THEN
+                              e(kkf,jjf,iif) = 1E-15_wp
                          END IF
    
@@ -586,10 +620,8 @@
 
         if (myid==0) print *, '** Fine Initalized ** simulated_time:', simulated_time
-#endif
+
     CONTAINS
    
        SUBROUTINE interpolate_to_fine_w( tag )
-      
-#if defined( __parallel )
       
            USE arrays_3d
@@ -602,14 +634,26 @@
            IMPLICIT NONE
       
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  wprs, wprf
-      
-           INTEGER(iwp) :: nzbottom, nztop
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: eps
+           REAL(wp)                           :: alpha
+           REAL(wp)                           :: eminus
+           REAL(wp)                           :: edot
+           REAL(wp)                           :: eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprf
       
       
@@ -632,7 +676,7 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                       DO if = bottomx, topx
-      
-                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                       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
@@ -640,5 +684,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
                       
-                           wprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           wprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -657,8 +701,8 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
       
-                   DO if = nxl, nxr
-                       DO jf = bottomy, topy
+                   DO iif = nxl, nxr
+                       DO jjf = bottomy, topy
                       
-                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                           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
@@ -666,7 +710,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           wprs(k,jf,if) = eminus * wprf(k,j-1,if) &
-                               + edot  * wprf(k,j,if)   &
-                               + eplus  * wprf(k,j+1,if)
+                           wprs(k,jjf,iif) = eminus * wprf(k,j-1,iif) &
+                               + edot  * wprf(k,j,iif)   &
+                               + eplus  * wprf(k,j+1,iif)
       
                        END DO
@@ -684,10 +728,10 @@
                topz    = (dzc/dzf) * (k+1) - 1
       
-               DO jf = nys, nyn
-                   DO if = nxl, nxr
-                       DO kf = bottomz, topz
-      
-                           w(kf,jf,if) = wprs(k,jf,if) + ( zwf(kf) - zwc(k) ) &
-                               * ( wprs(k+1,jf,if) - wprs(k,jf,if) ) / dzc
+               DO jjf = nys, nyn
+                   DO iif = nxl, nxr
+                       DO kkf = bottomz, topz
+      
+                           w(kkf,jjf,iif) = wprs(k,jjf,iif) + ( zwf(kkf) - zwc(k) ) &
+                               * ( wprs(k+1,jjf,iif) - wprs(k,jjf,iif) ) / dzc
       
                        END DO
@@ -697,8 +741,8 @@
            END DO
       
-           DO jf = nys, nyn
-               DO if = nxl, nxr
-      
-                   w(nzt,jf,if) = wprs(nztop,jf,if)
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
+      
+                   w(nzt,jjf,iif) = wprs(nztop,jjf,iif)
        
                END DO
@@ -709,10 +753,8 @@
            DEALLOCATE( wprf, wprs )
       
-#endif
        END SUBROUTINE interpolate_to_fine_w
       
        SUBROUTINE interpolate_to_fine_u( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -725,14 +767,27 @@
            IMPLICIT NONE
       
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprs, uprf
-      
-           INTEGER(iwp) :: nzbottom, nztop
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: eps
+           REAL(wp)                           :: alpha
+           REAL(wp)                           :: eminus
+           REAL(wp)                           :: edot
+           REAL(wp)                           :: eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprs
+      
       
       
@@ -756,7 +811,7 @@
       
                    DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
-                       DO jf = bottomy, topy
-      
-                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                       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
@@ -764,5 +819,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           uprf(k,jf,i) = eminus * work3d(k,j-1,i) &
+                           uprf(k,jjf,i) = eminus * work3d(k,j-1,i) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j+1,i)
@@ -782,9 +837,9 @@
                topz    = (dzc/dzf) * k
       
-               DO jf = nys, nyn
+               DO jjf = nys, nyn
                    DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
-                       DO kf = bottomz, topz
+                       DO kkf = bottomz, topz
                       
-                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
                            alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
                            eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
@@ -792,7 +847,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           uprs(kf,jf,i) = eminus * uprf(k-1,jf,i) &
-                               + edot  * uprf(k,jf,i)   &
-                               + eplus  * uprf(k+1,jf,i)
+                           uprs(kkf,jjf,i) = eminus * uprf(k-1,jjf,i) &
+                               + edot  * uprf(k,jjf,i)   &
+                               + eplus  * uprf(k+1,jjf,i)
       
                        END DO
@@ -802,5 +857,5 @@
            END DO
       
-           DO jf = nys, nyn
+           DO jjf = nys, nyn
                DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
       
@@ -811,7 +866,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                   uprs(nzt+1,jf,i)  = eminus * uprf(nztop,jf,i)   &
-                       + edot  * uprf(nztop+1,jf,i) &
-                       + eplus  * uprf(nztop+2,jf,i)
+                   uprs(nzt+1,jjf,i)  = eminus * uprf(nztop,jjf,i)   &
+                       + edot  * uprf(nztop+1,jjf,i) &
+                       + eplus  * uprf(nztop+2,jjf,i)
       
                END DO
@@ -821,6 +876,6 @@
            !-- Interpolation in x-direction (linear)
       
-           DO kf = nzb+1, nzt+1
-               DO jf = nys, nyn
+           DO kkf = nzb+1, nzt+1
+               DO jjf = nys, nyn
                    DO i = bdims_rem(1,1), bdims_rem(1,2)
       
@@ -828,7 +883,7 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                       DO if = bottomx, topx
-                           u(kf,jf,if)  = uprs(kf,jf,i) + ( if * dxf - i * dxc ) &
-                               * ( uprs(kf,jf,i+1) - uprs(kf,jf,i) ) / dxc
+                       DO iif = bottomx, topx
+                           u(kkf,jjf,iif)  = uprs(kkf,jjf,i) + ( iif * dxf - i * dxc ) &
+                               * ( uprs(kkf,jjf,i+1) - uprs(kkf,jjf,i) ) / dxc
                        END DO
       
@@ -843,5 +898,4 @@
            DEALLOCATE( uprf, uprs )
       
-#endif
        END SUBROUTINE interpolate_to_fine_u
       
@@ -849,5 +903,4 @@
        SUBROUTINE interpolate_to_fine_v( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -859,16 +912,28 @@
       
            IMPLICIT NONE
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: eps
+           REAL(wp)                           :: alpha
+           REAL(wp)                           :: eminus
+           REAL(wp)                           :: edot
+           REAL(wp)                           :: eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprs  
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprf
           
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprs, vprf
-      
-           INTEGER(iwp) :: nzbottom, nztop
-      
       
            nzbottom = bdims_rem (3,1)
@@ -890,7 +955,7 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                       DO if = bottomx, topx
-      
-                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                       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
@@ -898,5 +963,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           vprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           vprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -917,8 +982,8 @@
       
                DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
-                   DO if = nxl, nxr
-                       DO kf = bottomz, topz
-      
-                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                   DO iif = nxl, nxr
+                       DO kkf = bottomz, topz
+      
+                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
                            alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
                            eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
@@ -926,7 +991,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           vprs(kf,j,if) = eminus * vprf(k-1,j,if) &
-                               + edot  * vprf(k,j,if)   &
-                               + eplus  * vprf(k+1,j,if)
+                           vprs(kkf,j,iif) = eminus * vprf(k-1,j,iif) &
+                               + edot  * vprf(k,j,iif)   &
+                               + eplus  * vprf(k+1,j,iif)
       
                        END DO
@@ -937,5 +1002,5 @@
              
            DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
-               DO if = nxl, nxr
+               DO iif = nxl, nxr
       
                    eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
@@ -945,7 +1010,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                   vprs(nzt+1,j,if)  = eminus * vprf(nztop,j,if)   &
-                       + edot  * vprf(nztop+1,j,if) &
-                       + eplus  * vprf(nztop+2,j,if)
+                   vprs(nzt+1,j,iif)  = eminus * vprf(nztop,j,iif)   &
+                       + edot  * vprf(nztop+1,j,iif) &
+                       + eplus  * vprf(nztop+2,j,iif)
       
                END DO
@@ -955,5 +1020,5 @@
            !-- Interpolation in y-direction (linear)
       
-           DO kf = nzb+1, nzt+1
+           DO kkf = nzb+1, nzt+1
                DO j = bdims_rem(2,1), bdims_rem(2,2)
       
@@ -961,8 +1026,8 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
       
-                   DO if = nxl, nxr
-                       DO jf = bottomy, topy
-                           v (kf,jf,if) = vprs(kf,j,if) + ( jf * dyf - j * dyc ) &
-                               * ( vprs(kf,j+1,if) - vprs(kf,j,if) ) / dyc
+                   DO iif = nxl, nxr
+                       DO jjf = bottomy, topy
+                           v (kkf,jjf,iif) = vprs(kkf,j,iif) + ( jjf * dyf - j * dyc ) &
+                               * ( vprs(kkf,j+1,iif) - vprs(kkf,j,iif) ) / dyc
                        END DO
                    END DO
@@ -977,5 +1042,4 @@
            DEALLOCATE( vprf, vprs )
       
-#endif
        END SUBROUTINE interpolate_to_fine_v
       
@@ -983,5 +1047,4 @@
        SUBROUTINE interpolate_to_fine_s( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -993,16 +1056,27 @@
       
            IMPLICIT NONE
-      
-          
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprs, ptprf
-      
-           INTEGER(iwp) :: nzbottom, nztop
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: eps
+           REAL(wp)                           :: alpha
+           REAL(wp)                           :: eminus
+           REAL(wp)                           :: edot
+           REAL(wp)                           :: eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
       
       
@@ -1026,7 +1100,7 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                       DO if = bottomx, topx
-      
-                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                       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
@@ -1034,5 +1108,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -1052,8 +1126,8 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
       
-                   DO if = nxl, nxr
-                       DO jf = bottomy, topy
-      
-                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                   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
@@ -1061,7 +1135,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
-                               + edot  * ptprf(k,j,if)   &
-                               + eplus  * ptprf(k,j+1,if)
+                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+                               + edot  * ptprf(k,j,iif)   &
+                               + eplus  * ptprf(k,j+1,iif)
       
                        END DO
@@ -1079,9 +1153,9 @@
                topz    = (dzc/dzf) * k
       
-               DO jf = nys, nyn
-                   DO if = nxl, nxr
-                       DO kf = bottomz, topz
+               DO jjf = nys, nyn
+                   DO iif = nxl, nxr
+                       DO kkf = bottomz, topz
                        
-                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
                            alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
                            eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
@@ -1089,7 +1163,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           interpol3d(kf,jf,if) = eminus * ptprs(k-1,jf,if) &
-                               + edot  * ptprs(k,jf,if)   &
-                               + eplus  * ptprs(k+1,jf,if)
+                           interpol3d(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
+                               + edot  * ptprs(k,jjf,iif)   &
+                               + eplus  * ptprs(k+1,jjf,iif)
       
                        END DO
@@ -1099,6 +1173,6 @@
            END DO
                
-           DO jf = nys, nyn
-               DO if = nxl, nxr
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
       
                    eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
@@ -1108,7 +1182,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                   interpol3d(nzt+1,jf,if) = eminus * ptprs(nztop,jf,if)   &
-                       + edot  * ptprs(nztop+1,jf,if) &
-                       + eplus  * ptprs(nztop+2,jf,if)
+                   interpol3d(nzt+1,jjf,iif) = eminus * ptprs(nztop,jjf,iif)   &
+                       + edot  * ptprs(nztop+1,jjf,iif) &
+                       + eplus  * ptprs(nztop+2,jjf,iif)
       
                END DO
@@ -1118,5 +1192,4 @@
            DEALLOCATE( ptprf, ptprs ) 
       
-#endif
        END SUBROUTINE interpolate_to_fine_s
       
@@ -1124,5 +1197,4 @@
        SUBROUTINE interpolate_to_fine_kh( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -1134,24 +1206,43 @@
       
            IMPLICIT NONE
-      
-          
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprs, uprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprs, vprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  wprs, wprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprs, ptprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  eprs, eprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kmprs, kmprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  khprs, khprf
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  shfpr, uswspr, vswspr
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tspr, uspr, z0pr
-      
-           INTEGER(iwp) :: nzbottom, nztop
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: eps
+           REAL(wp)                           :: alpha
+           REAL(wp)                           :: eminus
+           REAL(wp)                           :: edot
+           REAL(wp)                           :: eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: eprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: kmprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: khprs
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: tspr 
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: eprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: kmprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: khprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: vswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uspr
       
       
@@ -1178,7 +1269,7 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                       DO if = bottomx, topx
-      
-                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                       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
@@ -1186,5 +1277,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -1204,8 +1295,8 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
       
-                   DO if = nxl, nxr
-                       DO jf = bottomy, topy
-      
-                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                   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
@@ -1213,7 +1304,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
-                               + edot  * ptprf(k,j,if)   &
-                               + eplus  * ptprf(k,j+1,if)
+                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+                               + edot  * ptprf(k,j,iif)   &
+                               + eplus  * ptprf(k,j+1,iif)
       
                        END DO
@@ -1231,9 +1322,9 @@
                topz    = (dzc/dzf) * k
       
-               DO jf = nys, nyn
-                   DO if = nxl, nxr
-                       DO kf = bottomz, topz
+               DO jjf = nys, nyn
+                   DO iif = nxl, nxr
+                       DO kkf = bottomz, topz
                        
-                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
                            alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
                            eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
@@ -1241,7 +1332,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           kh(kf,jf,if) = eminus * ptprs(k-1,jf,if) &
-                               + edot  * ptprs(k,jf,if)   &
-                               + eplus  * ptprs(k+1,jf,if)
+                           kh(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
+                               + edot  * ptprs(k,jjf,iif)   &
+                               + eplus  * ptprs(k+1,jjf,iif)
       
                        END DO
@@ -1251,6 +1342,6 @@
            END DO
                
-           DO jf = nys, nyn
-               DO if = nxl, nxr
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
       
                    eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
@@ -1260,7 +1351,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                   kh(nzt+1,jf,if) = eminus * ptprs(nztop,jf,if)   &
-                       + edot  * ptprs(nztop+1,jf,if) &
-                       + eplus  * ptprs(nztop+2,jf,if)
+                   kh(nzt+1,jjf,iif) = eminus * ptprs(nztop,jjf,iif)   &
+                       + edot  * ptprs(nztop+1,jjf,iif) &
+                       + eplus  * ptprs(nztop+2,jjf,iif)
       
                END DO
@@ -1270,10 +1361,8 @@
            DEALLOCATE( ptprf, ptprs ) 
       
-#endif
        END SUBROUTINE interpolate_to_fine_kh
       
        SUBROUTINE interpolate_to_fine_km( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -1285,24 +1374,42 @@
       
            IMPLICIT NONE
-      
-          
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprs, uprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprs, vprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  wprs, wprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprs, ptprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  eprs, eprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kmprs, kmprf
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  khprs, khprf
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  shfpr, uswspr, vswspr
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tspr, uspr, z0pr
-      
-           INTEGER(iwp) :: nzbottom, nztop
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: eps
+           REAL(wp)                           :: alpha
+           REAL(wp)                           :: eminus
+           REAL(wp)                           :: edot
+           REAL(wp)                           :: eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprs 
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprs 
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprs 
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: eprs 
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: kmprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: khprs
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: eprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: kmprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: khprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: vswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: tspr 
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uspr
       
       
@@ -1329,7 +1436,7 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                       DO if = bottomx, topx
-      
-                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                       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
@@ -1337,5 +1444,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -1355,8 +1462,8 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
       
-                   DO if = nxl, nxr
-                       DO jf = bottomy, topy
-      
-                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                   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
@@ -1364,7 +1471,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
-                               + edot  * ptprf(k,j,if)   &
-                               + eplus  * ptprf(k,j+1,if)
+                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+                               + edot  * ptprf(k,j,iif)   &
+                               + eplus  * ptprf(k,j+1,iif)
       
                        END DO
@@ -1382,9 +1489,9 @@
                topz    = (dzc/dzf) * k
       
-               DO jf = nys, nyn
-                   DO if = nxl, nxr
-                       DO kf = bottomz, topz
+               DO jjf = nys, nyn
+                   DO iif = nxl, nxr
+                       DO kkf = bottomz, topz
                        
-                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
                            alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
                            eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
@@ -1392,7 +1499,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           km(kf,jf,if) = eminus * ptprs(k-1,jf,if) &
-                               + edot  * ptprs(k,jf,if)   &
-                               + eplus  * ptprs(k+1,jf,if)
+                           km(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
+                               + edot  * ptprs(k,jjf,iif)   &
+                               + eplus  * ptprs(k+1,jjf,iif)
       
                        END DO
@@ -1402,6 +1509,6 @@
            END DO
                
-           DO jf = nys, nyn
-               DO if = nxl, nxr
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
       
                    eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
@@ -1411,7 +1518,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                   km(nzt+1,jf,if)  = eminus * ptprs(nztop,jf,if)   &
-                       + edot  * ptprs(nztop+1,jf,if) &
-                       + eplus  * ptprs(nztop+2,jf,if)
+                   km(nzt+1,jjf,iif)  = eminus * ptprs(nztop,jjf,iif)   &
+                       + edot  * ptprs(nztop+1,jjf,iif) &
+                       + eplus  * ptprs(nztop+2,jjf,iif)
       
                END DO
@@ -1421,5 +1528,4 @@
            DEALLOCATE( ptprf, ptprs ) 
       
-#endif
        END SUBROUTINE interpolate_to_fine_km
       
@@ -1429,5 +1535,4 @@
        SUBROUTINE interpolate_to_fine_flux( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -1439,25 +1544,34 @@
       
            IMPLICIT NONE
-      
-          
-           INTEGER(iwp), intent(in) ::  tag
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf, kf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           INTEGER(iwp)                             ::  bottomz, topz
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  shfpr, uswspr, vswspr
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tspr, uspr, z0pr
-      
-      
-           INTEGER(iwp) :: nzbottom, nztop
-      
-           ALLOCATE( shfpr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: nzbottom
+           INTEGER(iwp)                       :: nztop
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           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) )
-           ALLOCATE( z0pr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
       
            !
@@ -1473,7 +1587,7 @@
                    topx    = (nxf+1)/(nxc+1) * (i+1) - 1
       
-                   DO if = bottomx, topx
-      
-                       eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                   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
@@ -1481,27 +1595,19 @@
                        eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
       
-                       shfpr(j,if)  = eminus * work2dshf(j,i-1) &
-                           + edot  * work2dshf(j,i)   &
-                           + eplus  * work2dshf(j,i+1)
-      
-                       uswspr(j,if) = eminus * work2dusws(j,i-1) &
+                       uswspr(j,iif) = eminus * work2dusws(j,i-1) &
                            + edot  * work2dusws(j,i)   &
                            + eplus  * work2dusws(j,i+1)
       
-                       vswspr(j,if) = eminus * work2dvsws(j,i-1) &
+                       vswspr(j,iif) = eminus * work2dvsws(j,i-1) &
                            + edot  * work2dvsws(j,i)   &
                            + eplus  * work2dvsws(j,i+1)
       
-                       tspr(j,if)   = eminus * work2dts(j,i-1) &
+                       tspr(j,iif)   = eminus * work2dts(j,i-1) &
                            + edot  * work2dts(j,i)   &
                            + eplus  * work2dts(j,i+1)
       
-                       uspr(j,if)   = eminus * work2dus(j,i-1) &
+                       uspr(j,iif)   = eminus * work2dus(j,i-1) &
                            + edot  * work2dus(j,i)   &
                            + eplus  * work2dus(j,i+1)
-      
-                       z0pr(j,if)   = eminus * work2dz0(j,i-1) &
-                           + edot  * work2dz0(j,i)   &
-                           + eplus  * work2dz0(j,i+1)
       
                    END DO
@@ -1518,40 +1624,32 @@
                topy    = (nyf+1)/(nyc+1) * (j+1) - 1
       
-               DO if = nxl, nxr
-                   DO jf = bottomy, topy
-      
-                       eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+               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
-      
-!-- WARNING
-!-- shf,z0 not interpolated
-!-- line commented in interpolate_to_fine_flux
-!-- FG needs to read it's own data file      
-!MERGE-WIP!                    surf_def_h(0)%shf(jf,if)  = eminus * shfpr(j-1,if) &
-!MERGE-WIP!                        + edot  * shfpr(j,if)   &
-!MERGE-WIP!                        + eplus  * shfpr(j+1,if)
-!MERGE-WIP    
-!MERGE-WIP                     surf_def_h(0)%usws(jf,if) = eminus * uswspr(j-1,if) &
-!MERGE-WIP                         + edot  * uswspr(j,if)   &
-!MERGE-WIP                         + eplus  * uswspr(j+1,if)
-!MERGE-WIP    
-!MERGE-WIP                     surf_def_h(0)%vsws(jf,if) = eminus * vswspr(j-1,if) &
-!MERGE-WIP                         + edot  * vswspr(j,if)   &
-!MERGE-WIP                         + eplus  * vswspr(j+1,if)
-!MERGE-WIP    
-!MERGE                      ts(jf,if)   = eminus * tspr(j-1,if) &
-!MERGE                          + edot  * tspr(j,if)   &
-!MERGE                          + eplus  * tspr(j+1,if)
-!MERGE     
-!MERGE                      us(jf,if)   = eminus * uspr(j-1,if) &
-!MERGE                          + edot  * uspr(j,if)   &
-!MERGE                          + eplus  * uspr(j+1,if)
-!MERGE     
-!MERGE!                     z0(jf,if)   = eminus * z0pr(j-1,if) &
-!MERGE!                         + edot  * z0pr(j,if)   &
-!MERGE!                         + eplus  * z0pr(j+1,if)
+ 
+!!
+!!--   TODO
+!--    variables are not compatible with the new surface layer module
+!   
+!                    surf_def_h(0)%usws(jjf,iif) = eminus * uswspr(j-1,if) &
+!                        + edot  * uswspr(j,iif)   &
+!                        + eplus  * uswspr(j+1,iif)
+!   
+!                    surf_def_h(0)%vsws(jjf,iif) = eminus * vswspr(j-1,if) &
+!                        + edot  * vswspr(j,iif)   &
+!                        + eplus  * vswspr(j+1,iif)
+!   
+!                    ts(jjf,iif)   = eminus * tspr(j-1,if) &
+!                        + edot  * tspr(j,iif)   &
+!                        + eplus  * tspr(j+1,iif)
+!
+!                    us(jjf,iif)   = eminus * uspr(j-1,if) &
+!                        + edot  * uspr(j,iif)   &
+!                        + eplus  * uspr(j+1,iif)
       
                    END DO
@@ -1561,15 +1659,16 @@
       
       
-           DEALLOCATE( shfpr, uswspr, vswspr )
-           DEALLOCATE( tspr, uspr, z0pr )
-      
-      
-#endif
+           DEALLOCATE( uswspr, vswspr )
+           DEALLOCATE( tspr, uspr )
+      
+      
        END SUBROUTINE interpolate_to_fine_flux
    
    
+#endif       
     END SUBROUTINE vnest_init_fine
    
     SUBROUTINE vnest_boundary_conds
+#if defined( __parallel )
         !------------------------------------------------------------------------------!
         ! Description:
@@ -1591,10 +1690,11 @@
         IMPLICIT NONE
     
-        INTEGER(iwp) ::  i, j, k
-        INTEGER(iwp) ::  if, jf
-    
-        REAL(wp)    ::  c_max, denom
-    
-#if defined( __parallel )
+        INTEGER(iwp)                          :: i
+        INTEGER(iwp)                          :: j
+        INTEGER(iwp)                          :: iif
+        INTEGER(iwp)                          :: jjf
+        REAL(wp)                              :: c_max
+        REAL(wp)                              :: denom
+    
     
         !
@@ -1707,7 +1807,7 @@
     
 !-- TKE Neumann BC for FG top
-            DO jf = nys, nyn
-                DO if = nxl, nxr
-                   e(nzt+1,jf,if) = e(nzt,jf,if)
+            DO jjf = nys, nyn
+                DO iif = nxl, nxr
+                   e(nzt+1,jjf,iif) = e(nzt,jjf,iif)
                 END DO
             END DO
@@ -1731,10 +1831,8 @@
     
     
-#endif
     CONTAINS
     
        SUBROUTINE vnest_set_topbc_w
        
-#if defined( __parallel )
        
            USE arrays_3d
@@ -1746,12 +1844,22 @@
        
            IMPLICIT NONE
-       
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  wprf
-         
+
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           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 :: wprf
+       
            
            ALLOCATE( wprf(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
@@ -1773,12 +1881,12 @@
                    topx    = (nxf+1)/(nxc+1) * (i+1) - 1
        
-                   DO if = bottomx, topx
-       
-                       eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                   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
-                       wprf(j,if) = eminus * work3d(bdims_rem(3,1),j,i-1) &
+                       wprf(j,iif) = eminus * work3d(bdims_rem(3,1),j,i-1) &
                            + edot  * work3d(bdims_rem(3,1),j,i)   &
                            + eplus  * work3d(bdims_rem(3,1),j,i+1)
@@ -1798,9 +1906,9 @@
                topy    = (nyf+1)/(nyc+1) * (j+1) - 1
        
-               DO if = nxl, nxr
-       
-                   DO jf = bottomy, topy
-       
-                       eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+               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
@@ -1812,7 +1920,7 @@
                        eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                       w(nzt,jf,if) = eminus * wprf(j-1,if) &
-                           + edot  * wprf(j,if)   &
-                           + eplus  * wprf(j+1,if)
+                       w(nzt,jjf,iif) = eminus * wprf(j-1,iif) &
+                           + edot  * wprf(j,iif)   &
+                           + eplus  * wprf(j+1,iif)
        
                    END DO
@@ -1823,5 +1931,4 @@
 
            DEALLOCATE( wprf )
-#endif 
        
        END SUBROUTINE vnest_set_topbc_w
@@ -1830,5 +1937,4 @@
        SUBROUTINE vnest_set_topbc_u
        
-#if defined( __parallel )
        
            USE arrays_3d
@@ -1840,14 +1946,22 @@
        
            IMPLICIT NONE
-       
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprf
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  uprs
-         
-           
+
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           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 :: uprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uprs
        
            ALLOCATE( uprf(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
@@ -1865,7 +1979,7 @@
        
                    DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
-                       DO jf = bottomy, topy
-       
-                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+                       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
@@ -1873,5 +1987,5 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           uprf(k,jf,i) = eminus * work3d(k,j-1,i) &
+                           uprf(k,jjf,i) = eminus * work3d(k,j-1,i) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j+1,i)
@@ -1885,5 +1999,5 @@
            !-- Interpolation in z-direction
        
-           DO jf = nys, nyn
+           DO jjf = nys, nyn
                DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
@@ -1892,7 +2006,7 @@
                    edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-                   uprs(jf,i) = eminus * uprf(bdims_rem(3,1),jf,i)   &
-                       + edot  * uprf(bdims_rem(3,1)+1,jf,i) &
-                       + eplus  * uprf(bdims_rem(3,1)+2,jf,i)
+                   uprs(jjf,i) = eminus * uprf(bdims_rem(3,1),jjf,i)   &
+                       + edot  * uprf(bdims_rem(3,1)+1,jjf,i) &
+                       + eplus  * uprf(bdims_rem(3,1)+2,jjf,i)
                END DO
            END DO
@@ -1901,5 +2015,5 @@
            !-- Interpolation in x-direction
        
-           DO jf = nys, nyn
+           DO jjf = nys, nyn
                DO i = bdims_rem(1,1), bdims_rem(1,2)
        
@@ -1907,6 +2021,6 @@
                    topx    = (nxf+1)/(nxc+1) * (i+1) - 1
        
-                   DO if = bottomx, topx
-                       u(nzt+1,jf,if) = uprs(jf,i) + ( if * dxf - i * dxc ) * ( uprs(jf,i+1) - uprs(jf,i) ) / dxc
+                   DO iif = bottomx, topx
+                       u(nzt+1,jjf,iif) = uprs(jjf,i) + ( iif * dxf - i * dxc ) * ( uprs(jjf,i+1) - uprs(jjf,i) ) / dxc
                    END DO
        
@@ -1917,5 +2031,4 @@
        
            DEALLOCATE ( uprf, uprs )
-#endif 
        
        END SUBROUTINE vnest_set_topbc_u
@@ -1924,5 +2037,4 @@
        SUBROUTINE vnest_set_topbc_v
        
-#if defined( __parallel )
        
            USE arrays_3d
@@ -1934,12 +2046,22 @@
        
            IMPLICIT NONE
-       
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprf
-           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  vprs
+
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           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 :: vprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: vprs
          
            
@@ -1963,12 +2085,12 @@
                        topx    = (nxf+1)/(nxc+1) * (i+1) - 1
        
-                       DO if = bottomx, topx
-       
-                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                       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
-                           vprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           vprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -1983,5 +2105,5 @@
        
            DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
-               DO if = nxl, nxr
+               DO iif = nxl, nxr
        
                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
@@ -1990,7 +2112,7 @@
                    edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
-                   vprs(j,if) = eminus * vprf(bdims_rem(3,1),j,if)   &
-                       + edot  * vprf(bdims_rem(3,1)+1,j,if) &
-                       + eplus  * vprf(bdims_rem(3,1)+2,j,if)
+                   vprs(j,iif) = eminus * vprf(bdims_rem(3,1),j,iif)   &
+                       + edot  * vprf(bdims_rem(3,1)+1,j,iif) &
+                       + eplus  * vprf(bdims_rem(3,1)+2,j,iif)
        
                END DO
@@ -2001,12 +2123,12 @@
        
            DO j = bdims_rem(2,1), bdims_rem(2,2)
-               DO if = nxl, nxr
+               DO iif = nxl, nxr
        
                    bottomy = (nyf+1)/(nyc+1) * j
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
        
-                   DO jf = bottomy, topy
-       
-                       v(nzt+1,jf,if) = vprs(j,if) + ( jf * dyf - j * dyc ) * ( vprs(j+1,if) - vprs(j,if) ) / dyc
+                   DO jjf = bottomy, topy
+       
+                       v(nzt+1,jjf,iif) = vprs(j,iif) + ( jjf * dyf - j * dyc ) * ( vprs(j+1,iif) - vprs(j,iif) ) / dyc
        
                    END DO
@@ -2018,5 +2140,4 @@
        
        
-#endif 
        
        END SUBROUTINE vnest_set_topbc_v
@@ -2025,5 +2146,4 @@
        SUBROUTINE vnest_set_topbc_s
        
-#if defined( __parallel )
        
            USE arrays_3d
@@ -2035,12 +2155,23 @@
        
            IMPLICIT NONE
-       
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprf, ptprs
-         
+
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           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
+       
            
        
@@ -2064,7 +2195,7 @@
                        topx    = (nxf+1)/(nxc+1) *(i+1) - 1
        
-                       DO if = bottomx, topx
-       
-                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                       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
@@ -2076,5 +2207,5 @@
                            eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -2097,9 +2228,9 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
        
-                   DO if = nxl, nxr
-       
-                       DO jf = bottomy, topy
-       
-                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+                   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
@@ -2111,7 +2242,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
-                               + edot  * ptprf(k,j,if)   &
-                               + eplus  * ptprf(k,j+1,if)
+                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+                               + edot  * ptprf(k,j,iif)   &
+                               + eplus  * ptprf(k,j+1,iif)
                        END DO
        
@@ -2125,6 +2256,6 @@
            !-- Interpolation in z-direction
        
-           DO jf = nys, nyn
-               DO if = nxl, nxr
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
        
                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
@@ -2138,7 +2269,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                   interpol3d (nzt+1,jf,if) = eminus * ptprs(bdims_rem(3,1),jf,if)   &
-                       + edot  * ptprs(bdims_rem(3,1)+1,jf,if) &
-                       + eplus  * ptprs(bdims_rem(3,1)+2,jf,if)
+                   interpol3d (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
@@ -2148,11 +2279,12 @@
        
        
-#endif 
        
        END SUBROUTINE vnest_set_topbc_s
+#endif
     END SUBROUTINE vnest_boundary_conds
     
  
     SUBROUTINE vnest_boundary_conds_khkm
+#if defined( __parallel )
     
         !--------------------------------------------------------------------------------!
@@ -2175,10 +2307,11 @@
         IMPLICIT NONE
     
-        INTEGER(iwp) ::  i, j, k
-        INTEGER(iwp) ::  if, jf
-    
-        REAL(wp)    ::  c_max, denom
-    
-#if defined( __parallel )
+        INTEGER(iwp)                          :: i
+        INTEGER(iwp)                          :: j
+        INTEGER(iwp)                          :: iif
+        INTEGER(iwp)                          :: jjf
+        REAL(wp)                              :: c_max
+        REAL(wp)                              :: denom
+    
     
         IF ( coupling_mode == 'vnested_crse' )  THEN
@@ -2265,7 +2398,7 @@
     
         ! Neumann BC for FG kh
-        DO jf = nys, nyn
-            DO if = nxl, nxr
-               kh(nzt+1,jf,if) = kh(nzt,jf,if)
+        DO jjf = nys, nyn
+            DO iif = nxl, nxr
+               kh(nzt+1,jjf,iif) = kh(nzt,jjf,iif)
             END DO
         END DO
@@ -2275,7 +2408,7 @@
     
         ! Neumann BC for FG kh
-        DO jf = nys, nyn
-            DO if = nxl, nxr
-               km(nzt+1,jf,if) = km(nzt,jf,if)
+        DO jjf = nys, nyn
+            DO iif = nxl, nxr
+               km(nzt+1,jjf,iif) = km(nzt,jjf,iif)
             END DO
         END DO
@@ -2284,9 +2417,9 @@
         !
         !-- The following evaluation can only be performed, if the fine grid is situated below the inversion
-        !!    DO jf = nys-1, nyn+1
-        !!       DO if = nxl-1, nxr+1
+        !!    DO jjf = nys-1, nyn+1
+        !!       DO iif = nxl-1, nxr+1
         !!
-        !!          km(nzt+1,jf,if) = 0.1 * l_grid(nzt+1) * SQRT( e(nzt+1,jf,if) )
-        !!          kh(nzt+1,jf,if) = 3.0 * km(nzt+1,jf,if)
+        !!          km(nzt+1,jjf,iif) = 0.1 * l_grid(nzt+1) * SQRT( e(nzt+1,jjf,iif) )
+        !!          kh(nzt+1,jjf,iif) = 3.0 * km(nzt+1,jjf,iif)
         !!
         !!       END DO
@@ -2300,5 +2433,4 @@
         ENDIF
     
-#endif 
     
     CONTAINS
@@ -2306,5 +2438,4 @@
        SUBROUTINE vnest_set_topbc_kh
        
-#if defined( __parallel )
        
            USE arrays_3d
@@ -2316,12 +2447,23 @@
        
            IMPLICIT NONE
-       
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprf, ptprs
-         
+
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           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
+       
            
        
@@ -2345,7 +2487,7 @@
                        topx    = (nxf+1)/(nxc+1) *(i+1) - 1
        
-                       DO if = bottomx, topx
-       
-                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                       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
@@ -2357,5 +2499,5 @@
                            eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -2378,9 +2520,9 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
        
-                   DO if = nxl, nxr
-       
-                       DO jf = bottomy, topy
-       
-                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+                   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
@@ -2392,7 +2534,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
-                               + edot  * ptprf(k,j,if)   &
-                               + eplus  * ptprf(k,j+1,if)
+                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+                               + edot  * ptprf(k,j,iif)   &
+                               + eplus  * ptprf(k,j+1,iif)
                        END DO
        
@@ -2406,6 +2548,6 @@
            !-- Interpolation in z-direction
        
-           DO jf = nys, nyn
-               DO if = nxl, nxr
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
        
                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
@@ -2419,7 +2561,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                   kh (nzt+1,jf,if) = eminus * ptprs(bdims_rem(3,1),jf,if)   &
-                       + edot  * ptprs(bdims_rem(3,1)+1,jf,if) &
-                       + eplus  * ptprs(bdims_rem(3,1)+2,jf,if)
+                   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
@@ -2429,5 +2571,4 @@
        
        
-#endif 
        
        END SUBROUTINE vnest_set_topbc_kh
@@ -2435,5 +2576,4 @@
        SUBROUTINE vnest_set_topbc_km
        
-#if defined( __parallel )
        
            USE arrays_3d
@@ -2445,12 +2585,24 @@
        
            IMPLICIT NONE
-       
-           INTEGER(iwp)                             ::  i, j, k
-           INTEGER(iwp)                             ::  if, jf
-           INTEGER(iwp)                             ::  bottomx, topx
-           INTEGER(iwp)                             ::  bottomy, topy
-           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
-           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprf, ptprs
-         
+
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           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
+       
            
        
@@ -2474,7 +2626,7 @@
                        topx    = (nxf+1)/(nxc+1) *(i+1) - 1
        
-                       DO if = bottomx, topx
-       
-                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                       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
@@ -2486,5 +2638,5 @@
                            eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
                                + edot  * work3d(k,j,i)   &
                                + eplus  * work3d(k,j,i+1)
@@ -2507,9 +2659,9 @@
                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
        
-                   DO if = nxl, nxr
-       
-                       DO jf = bottomy, topy
-       
-                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+                   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
@@ -2521,7 +2673,7 @@
                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
-                               + edot  * ptprf(k,j,if)   &
-                               + eplus  * ptprf(k,j+1,if)
+                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
+                               + edot  * ptprf(k,j,iif)   &
+                               + eplus  * ptprf(k,j+1,iif)
                        END DO
        
@@ -2535,6 +2687,6 @@
            !-- Interpolation in z-direction
        
-           DO jf = nys, nyn
-               DO if = nxl, nxr
+           DO jjf = nys, nyn
+               DO iif = nxl, nxr
        
                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
@@ -2548,7 +2700,7 @@
                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
        
-                   km (nzt+1,jf,if) = eminus * ptprs(bdims_rem(3,1),jf,if)   &
-                       + edot  * ptprs(bdims_rem(3,1)+1,jf,if) &
-                       + eplus  * ptprs(bdims_rem(3,1)+2,jf,if)
+                   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
@@ -2558,9 +2710,9 @@
        
        
-#endif 
        
        END SUBROUTINE vnest_set_topbc_km
     
  
+#endif
     END SUBROUTINE vnest_boundary_conds_khkm
    
@@ -2568,4 +2720,6 @@
    
     SUBROUTINE vnest_anterpolate
+
+#if defined( __parallel )
    
         !--------------------------------------------------------------------------------!
@@ -2587,14 +2741,15 @@
         IMPLICIT NONE
    
-        REAL(wp)    ::  time_since_reference_point_rem
-        INTEGER(iwp)   :: i, j, k, im, jn, ko
-
-!---    INTEGER(iwp) ::  j  !< grid index y direction
-!--     INTEGER(iwp) ::  k  !< grid index z direction
-        INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing. 
-        INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
-        INTEGER(iwp) ::  m  !< running index surface elements
-   
-#if defined( __parallel )
+        REAL(wp)                              ::  time_since_reference_point_rem
+        INTEGER(iwp)                          ::  i
+        INTEGER(iwp)                          ::  j
+        INTEGER(iwp)                          ::  k
+        INTEGER(iwp)                          ::  im
+        INTEGER(iwp)                          ::  jn
+        INTEGER(iwp)                          ::  ko
+        INTEGER(iwp)                          ::  kb !< variable to set respective boundary value, depends on facing. 
+        INTEGER(iwp)                          ::  l  !< running index boundary type, for up- and downward-facing walls
+        INTEGER(iwp)                          ::  m  !< running index surface elements
+   
    
    
@@ -2835,10 +2990,8 @@
    
    
-#endif
 
     CONTAINS
        SUBROUTINE anterpolate_to_crse_u( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -2850,12 +3003,19 @@
       
            IMPLICIT NONE
-      
-           INTEGER(iwp)               ::  i, j, k
-           INTEGER(iwp)               ::  if, jf, kf
-           INTEGER(iwp)               ::  bottomx, topx
-           INTEGER(iwp)               ::  bottomy, topy
-           INTEGER(iwp)               ::  bottomz, topz
-           REAL(wp)                  ::  aweight
-           INTEGER(iwp), intent(in)   ::  tag
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: aweight
       
            !
@@ -2876,12 +3036,12 @@
                    DO i = bdims_rem(1,1),bdims_rem(1,2)
       
-                       if = (nxf+1) / (nxc+1) * i
+                       iif = (nxf+1) / (nxc+1) * i
       
                        aweight   = 0.0
       
-                       DO kf = bottomz, topz
-                           DO jf = bottomy, topy
-      
-                               aweight   = aweight + anterpol3d(kf,jf,if) *        &
+                       DO kkf = bottomz, topz
+                           DO jjf = bottomy, topy
+      
+                               aweight   = aweight + anterpol3d(kkf,jjf,iif) *        &
                                    (dzf/dzc) * (dyf/dyc)
       
@@ -2898,5 +3058,4 @@
       
       
-#endif
       
        END SUBROUTINE anterpolate_to_crse_u
@@ -2905,5 +3064,4 @@
        SUBROUTINE anterpolate_to_crse_v( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -2915,12 +3073,20 @@
       
            IMPLICIT NONE
-      
-           INTEGER(iwp)               ::  i, j, k
-           INTEGER(iwp)               ::  if, jf, kf
-           INTEGER(iwp)               ::  bottomx, topx
-           INTEGER(iwp)               ::  bottomy, topy
-           INTEGER(iwp)               ::  bottomz, topz
-           REAL(wp)                  ::  aweight
-           INTEGER(iwp), intent(in) ::  tag
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: aweight
+
            !
            !-- Anterpolation of the velocity components v
@@ -2935,5 +3101,5 @@
                DO j = bdims_rem(2,1), bdims_rem(2,2)
       
-                   jf = (nyf+1) / (nyc+1) * j
+                   jjf = (nyf+1) / (nyc+1) * j
       
                    DO i = bdims_rem(1,1), bdims_rem(1,2)
@@ -2944,8 +3110,8 @@
                        aweight   = 0.0
       
-                       DO kf = bottomz, topz
-                           DO if = bottomx, topx
-      
-                               aweight   = aweight + anterpol3d(kf,jf,if) *        &
+                       DO kkf = bottomz, topz
+                           DO iif = bottomx, topx
+      
+                               aweight   = aweight + anterpol3d(kkf,jjf,iif) *        &
                                    (dzf/dzc) * (dxf/dxc)
       
@@ -2961,5 +3127,4 @@
       
       
-#endif
       
        END SUBROUTINE anterpolate_to_crse_v
@@ -2968,5 +3133,4 @@
        SUBROUTINE anterpolate_to_crse_w( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -2978,12 +3142,20 @@
       
            IMPLICIT NONE
-      
-           INTEGER(iwp)               ::  i, j, k
-           INTEGER(iwp)               ::  if, jf, kf
-           INTEGER(iwp)               ::  bottomx, topx
-           INTEGER(iwp)               ::  bottomy, topy
-           INTEGER(iwp)               ::  bottomz, topz
-           REAL(wp)                  ::  aweight
-           INTEGER(iwp), intent(in) ::  tag
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: aweight
+
            !
            !-- Anterpolation of the velocity components w
@@ -2993,5 +3165,5 @@
            DO k = bdims_rem(3,1), bdims_rem(3,2)-1
       
-               kf = cfratio(3) * k
+               kkf = cfratio(3) * k
       
                DO j = bdims_rem(2,1), bdims_rem(2,2)
@@ -3007,8 +3179,8 @@
                        aweight   = 0.0
       
-                       DO jf = bottomy, topy
-                           DO if = bottomx, topx
-      
-                               aweight   = aweight + anterpol3d (kf,jf,if) *        &
+                       DO jjf = bottomy, topy
+                           DO iif = bottomx, topx
+      
+                               aweight   = aweight + anterpol3d (kkf,jjf,iif) *        &
                                    (dxf/dxc) * (dyf/dyc)
       
@@ -3024,5 +3196,4 @@
            END DO
       
-#endif
       
        END SUBROUTINE anterpolate_to_crse_w
@@ -3031,5 +3202,4 @@
        SUBROUTINE anterpolate_to_crse_s( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -3041,12 +3211,19 @@
       
            IMPLICIT NONE
-      
-           INTEGER(iwp)               ::  i, j, k
-           INTEGER(iwp)               ::  if, jf, kf
-           INTEGER(iwp)               ::  bottomx, topx
-           INTEGER(iwp)               ::  bottomy, topy
-           INTEGER(iwp)               ::  bottomz, topz
-           REAL(wp)                  ::  aweight
-           INTEGER(iwp), intent(in)   ::  tag
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: aweight
       
            !
@@ -3071,9 +3248,9 @@
                        aweight   = 0.0
           
-                       DO kf = bottomz, topz
-                           DO jf = bottomy, topy
-                               DO if = bottomx, topx
-      
-                                   aweight   = aweight  + anterpol3d(kf,jf,if) *                  &
+                       DO kkf = bottomz, topz
+                           DO jjf = bottomy, topy
+                               DO iif = bottomx, topx
+      
+                                   aweight   = aweight  + anterpol3d(kkf,jjf,iif) *                  &
                                        (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
       
@@ -3090,7 +3267,7 @@
            END DO
       
-#endif
       
        END SUBROUTINE anterpolate_to_crse_s
+#endif       
     END SUBROUTINE vnest_anterpolate
    
@@ -3098,4 +3275,5 @@
    
     SUBROUTINE vnest_anterpolate_e
+#if defined( __parallel )
    
         !--------------------------------------------------------------------------------!
@@ -3115,8 +3293,12 @@
         IMPLICIT NONE
    
-        REAL(wp)    ::  time_since_reference_point_rem
-        INTEGER(iwp)   :: i, j, k, im, jn, ko
-   
-#if defined( __parallel )
+        REAL(wp)                              :: time_since_reference_point_rem
+        INTEGER(iwp)                          :: i
+        INTEGER(iwp)                          :: j
+        INTEGER(iwp)                          :: k
+        INTEGER(iwp)                          :: im
+        INTEGER(iwp)                          :: jn
+        INTEGER(iwp)                          :: ko
+   
    
         !
@@ -3243,5 +3425,4 @@
         ENDIF
    
-#endif
    
     CONTAINS 
@@ -3253,5 +3434,4 @@
        SUBROUTINE anterpolate_to_crse_e( tag )
       
-#if defined( __parallel )
       
            USE arrays_3d
@@ -3263,14 +3443,24 @@
       
            IMPLICIT NONE
-      
-           INTEGER(iwp)               ::  i, j, k
-           INTEGER(iwp)               ::  if, jf, kf
-           INTEGER(iwp)               ::  bottomx, topx
-           INTEGER(iwp)               ::  bottomy, topy
-           INTEGER(iwp)               ::  bottomz, topz
-           REAL(wp)                  ::  aweight_a, aweight_b, aweight_c, aweight_d, aweight_e
-           REAL(wp)                  ::  energ
-           INTEGER(iwp), intent(in)   ::  tag
-      
+
+           INTEGER(iwp), intent(in)           :: tag
+           INTEGER(iwp)                       :: i
+           INTEGER(iwp)                       :: j
+           INTEGER(iwp)                       :: k
+           INTEGER(iwp)                       :: iif
+           INTEGER(iwp)                       :: jjf
+           INTEGER(iwp)                       :: kkf
+           INTEGER(iwp)                       :: bottomx
+           INTEGER(iwp)                       :: bottomy
+           INTEGER(iwp)                       :: bottomz
+           INTEGER(iwp)                       :: topx
+           INTEGER(iwp)                       :: topy
+           INTEGER(iwp)                       :: topz
+           REAL(wp)                           :: aweight_a
+           REAL(wp)                           :: aweight_b
+           REAL(wp)                           :: aweight_c
+           REAL(wp)                           :: aweight_d
+           REAL(wp)                           :: aweight_e
+           REAL(wp)                           :: energ 
           
            DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
@@ -3295,26 +3485,26 @@
                        aweight_e   = 0.0
       
-                       DO kf = bottomz, topz
-                           DO jf = bottomy, topy
-                               DO if = bottomx, topx
-      
-                                   aweight_a = aweight_a + anterpol3d(kf,jf,if)  *                     &
+                       DO kkf = bottomz, topz
+                           DO jjf = bottomy, topy
+                               DO iif = bottomx, topx
+      
+                                   aweight_a = aweight_a + anterpol3d(kkf,jjf,iif)  *                     &
                                        (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
                               
       
-                                   energ = ( 0.5 * ( u(kf,jf,if)   + u(kf,jf,if+1) ) )**2.0 +          &
-                                       ( 0.5 * ( v(kf,jf,if)   + v(kf,jf+1,if) ) )**2.0 +              &
-                                       ( 0.5 * ( w(kf-1,jf,if) + w(kf,jf,if) ) )**2.0
+                                   energ = ( 0.5 * ( u(kkf,jjf,iif)   + u(kkf,jjf,iif+1) ) )**2.0 +          &
+                                       ( 0.5 * ( v(kkf,jjf,iif)   + v(kkf,jjf+1,iif) ) )**2.0 +              &
+                                       ( 0.5 * ( w(kkf-1,jjf,iif) + w(kkf,jjf,iif) ) )**2.0
       
                                    aweight_b   =   aweight_b + energ         *                         &
                                        (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
        
-                                   aweight_c   =   aweight_c + 0.5 * ( u(kf,jf,if) + u(kf,jf,if+1) ) * &
+                                   aweight_c   =   aweight_c + 0.5 * ( u(kkf,jjf,iif) + u(kkf,jjf,iif+1) ) * &
                                        (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
       
-                                   aweight_d   =   aweight_d + 0.5 * ( v(kf,jf,if) + v(kf,jf+1,if) ) * &
+                                   aweight_d   =   aweight_d + 0.5 * ( v(kkf,jjf,iif) + v(kkf,jjf+1,iif) ) * &
                                        (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
        
-                                   aweight_e   =   aweight_e + 0.5 * ( w(kf-1,jf,if) + w(kf,jf,if) ) * &
+                                   aweight_e   =   aweight_e + 0.5 * ( w(kkf-1,jjf,iif) + w(kkf,jjf,iif) ) * &
                                        (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
       
@@ -3336,10 +3526,11 @@
       
        
-#endif
       
        END SUBROUTINE anterpolate_to_crse_e
+#endif       
     END SUBROUTINE vnest_anterpolate_e
 
     SUBROUTINE vnest_init_pegrid_rank
+#if defined( __parallel )
 ! Domain decomposition and exchange of grid variables between coarse and fine
 ! Given processor coordinates as index f_rnk_lst(pcoord(1), pcoord(2))
@@ -3369,6 +3560,6 @@
        IMPLICIT NONE
 
-       INTEGER(iwp)                    :: dest_rnk
-       INTEGER(iwp)               ::  i                        !<
+       INTEGER(iwp)                           :: dest_rnk
+       INTEGER(iwp)                           ::  i                        !<
 
        IF (myid == 0) THEN
@@ -3429,4 +3620,5 @@
 
 
+#endif
  
     END SUBROUTINE vnest_init_pegrid_rank
@@ -3434,4 +3626,5 @@
 
     SUBROUTINE vnest_init_pegrid_domain
+#if defined( __parallel )
 
        USE control_parameters,                                                    &
@@ -3451,9 +3644,17 @@
        IMPLICIT NONE
 
-       INTEGER(iwp)              :: dest_rnk
-       INTEGER(iwp) ::  i, j                               !<
-       INTEGER(iwp) ::  tempx, tempy                               !<
-       INTEGER(iwp) :: TYPE_INT_YZ,  SIZEOFREAL
-       INTEGER(iwp) :: MTV_X,MTV_Y,MTV_Z,MTV_RX,MTV_RY,MTV_RZ
+       INTEGER(iwp)                           :: dest_rnk
+       INTEGER(iwp)                           :: i              !<
+       INTEGER(iwp)                           :: j              !<
+       INTEGER(iwp)                           :: tempx
+       INTEGER(iwp)                           :: tempy
+       INTEGER(iwp)                           :: TYPE_INT_YZ
+       INTEGER(iwp)                           :: SIZEOFREAL
+       INTEGER(iwp)                           :: MTV_X
+       INTEGER(iwp)                           :: MTV_Y
+       INTEGER(iwp)                           :: MTV_Z
+       INTEGER(iwp)                           :: MTV_RX
+       INTEGER(iwp)                           :: MTV_RY
+       INTEGER(iwp)                           :: MTV_RZ
    
        !
@@ -3692,5 +3893,5 @@
 
         ENDIF
-   
+#endif   
        END SUBROUTINE vnest_init_pegrid_domain
 
@@ -3698,4 +3899,5 @@
        SUBROUTINE vnest_init_grid
  
+#if defined( __parallel )
           USE arrays_3d,                                                             &
               ONLY:  zu, zw
@@ -3761,28 +3963,19 @@
           ENDIF
 
+#endif
        END SUBROUTINE vnest_init_grid
 
 
        SUBROUTINE vnest_check_parameters
-
-          USE arrays_3d,                                                             &
-              ONLY:  zu, zw
-
-          USE control_parameters,                                                    &
-              ONLY:  coupling_mode
-             
-          USE indices,                                                               &
-              ONLY:  nzt
-         
-          USE kinds
-          
-          USE pegrid
+#if defined( __parallel )
+
+          USE pegrid,                                                                &
+              ONLY:  myid
 
           IMPLICIT NONE
 
-
           IF (myid==0) PRINT*, '*** vnest: check parameters not implemented yet ***'
       
-
+#endif
        END SUBROUTINE vnest_check_parameters
 
@@ -3790,6 +3983,7 @@
        SUBROUTINE vnest_timestep_sync
 
+#if defined( __parallel )
          USE control_parameters,                                                    &
-             ONLY:  coupling_mode, dt_3d, dt_coupling
+             ONLY:  coupling_mode, dt_3d, old_dt
       
          USE interfaces
@@ -3825,10 +4019,10 @@
 !-- Identical timestep for coarse and fine grids
           dt_3d = MIN( dtc, dtf )
-!-- Nest coupling at every timestep
-          dt_coupling = dt_3d
-
+          old_dt = dt_3d
+#endif
        END SUBROUTINE vnest_timestep_sync
        
        SUBROUTINE vnest_deallocate
+#if defined( __parallel )
           USE control_parameters,                                                    &
               ONLY:  coupling_mode
@@ -3846,5 +4040,5 @@
              IF ( ALLOCATED (f2c_dims_fg) ) DEALLOCATE (f2c_dims_fg)
           ENDIF
-       
+#endif 
        END SUBROUTINE vnest_deallocate
 
