Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 708)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 709)
@@ -1,8 +1,8 @@
  MODULE advec_ws
-
 
 !-----------------------------------------------------------------------------!
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
@@ -44,5 +44,4 @@
 !-----------------------------------------------------------------------------!
 
-
     PRIVATE
     PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
@@ -52,23 +51,29 @@
        MODULE PROCEDURE ws_init
     END INTERFACE ws_init
+
     INTERFACE ws_statistics
        MODULE PROCEDURE ws_statistics
     END INTERFACE ws_statistics
+
     INTERFACE advec_s_ws
        MODULE PROCEDURE advec_s_ws
        MODULE PROCEDURE advec_s_ws_ij
     END INTERFACE advec_s_ws
+
     INTERFACE advec_u_ws
        MODULE PROCEDURE advec_u_ws
        MODULE PROCEDURE advec_u_ws_ij
     END INTERFACE advec_u_ws
+
     INTERFACE advec_v_ws
        MODULE PROCEDURE advec_v_ws
        MODULE PROCEDURE advec_v_ws_ij
     END INTERFACE advec_v_ws
+
     INTERFACE advec_w_ws
        MODULE PROCEDURE advec_w_ws
        MODULE PROCEDURE advec_w_ws_ij
     END INTERFACE advec_w_ws
+
     INTERFACE local_diss
        MODULE PROCEDURE local_diss
@@ -78,8 +83,10 @@
  CONTAINS
 
-!
-!-- Initialisation
-
+
+!------------------------------------------------------------------------------!
+! Initialization of WS-scheme
+!------------------------------------------------------------------------------!
     SUBROUTINE ws_init
+
        USE arrays_3d
        USE constants
@@ -91,15 +98,16 @@
 !--    Allocate arrays needed for dissipation control.
        IF ( dissipation_control )  THEN
-   !          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),  &
-    !                  var_y(nzb+1:nzt,nys:nyn,nxl:nxr),  &
-     !                 var_z(nzb+1:nzt,nys:nyn,nxl:nxr),  &
-      !                gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-       !               gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-        !              gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+!          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),        &
+!                   var_y(nzb+1:nzt,nys:nyn,nxl:nxr),        &
+!                   var_z(nzb+1:nzt,nys:nyn,nxl:nxr),        &
+!                   gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+!                   gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+!                   gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
        ENDIF
-              
+
+!
 !--    Set the appropriate factors for scalar and momentum advection.
        adv_sca_5 = 0.016666666666666
-       adv_sca_3 = 0.083333333333333 
+       adv_sca_3 = 0.083333333333333
        adv_mom_5 = 0.0083333333333333
        adv_mom_3 = 0.041666666666666
@@ -108,10 +116,124 @@
 !--    Arrays needed for statical evaluation of fluxes.
        IF ( ws_scheme_mom )  THEN
-          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),         &
-              sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),               &
-              sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),                &
-              sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),                &
-              sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
-
+
+          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),  &
+                    sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),  &
+                    sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),   &
+                    sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),   &
+                    sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
+
+          sums_wsus_ws_l = 0.0
+          sums_wsvs_ws_l = 0.0
+          sums_us2_ws_l  = 0.0
+          sums_vs2_ws_l  = 0.0
+          sums_ws2_ws_l  = 0.0
+
+       ENDIF
+
+       IF ( ws_scheme_sca )  THEN
+
+          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions) )
+          sums_wspts_ws_l = 0.0
+
+          IF ( humidity .OR. passive_scalar )  THEN
+             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions) )
+             sums_wsqs_ws_l = 0.0
+          ENDIF
+
+          IF ( ocean )  THEN
+             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions) )
+             sums_wssas_ws_l = 0.0
+          ENDIF
+
+       ENDIF
+
+!
+!--    Arrays needed for reasons of speed optimization for cache and noopt 
+!--    version. For the vector version the buffer arrays are not necessary, 
+!--    because the the fluxes can swapped directly inside the loops of the 
+!--    advection routines.
+       IF ( loop_optimization /= 'vector' )  THEN
+
+          IF ( ws_scheme_mom )  THEN
+
+             ALLOCATE( flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),  &
+                       flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),  &
+                       diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt) )
+             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn),               &
+                       flux_l_v(nzb+1:nzt,nys:nyn),               &
+                       flux_l_w(nzb+1:nzt,nys:nyn),               &
+                       diss_l_u(nzb+1:nzt,nys:nyn),               &
+                       diss_l_v(nzb+1:nzt,nys:nyn),               &
+                       diss_l_w(nzb+1:nzt,nys:nyn) )
+
+          ENDIF
+
+          IF ( ws_scheme_sca )  THEN
+
+             ALLOCATE( flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), &
+                       diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt) ) 
+             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn),              &
+                       flux_l_e(nzb+1:nzt,nys:nyn),               &
+                       diss_l_pt(nzb+1:nzt,nys:nyn),              &
+                       diss_l_e(nzb+1:nzt,nys:nyn) )
+
+             IF ( humidity .OR. passive_scalar )  THEN
+                ALLOCATE( flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt) )
+                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn),            &
+                          diss_l_q(nzb+1:nzt,nys:nyn) )
+             ENDIF
+
+             IF ( ocean )  THEN
+                ALLOCATE( flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt) )
+                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn),           &
+                          diss_l_sa(nzb+1:nzt,nys:nyn) )
+             ENDIF
+
+          ENDIF
+
+       ENDIF
+
+!
+!--    Determine the flags where the order of the scheme for horizontal 
+!--    advection has to be degraded.
+       ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+
+             boundary_flags(j,i) = 0
+             IF ( outflow_l )  THEN
+                IF ( i == nxlu  )  boundary_flags(j,i) = 5
+                IF ( i == nxl   )  boundary_flags(j,i) = 6
+             ELSEIF ( outflow_r )  THEN
+                IF ( i == nxr-1 )  boundary_flags(j,i) = 1
+                IF ( i == nxr   )  boundary_flags(j,i) = 2
+             ELSEIF ( outflow_n )  THEN
+                IF ( j == nyn-1 )  boundary_flags(j,i) = 3
+                IF ( j == nyn   )  boundary_flags(j,i) = 4
+             ELSEIF ( outflow_s )  THEN
+                IF ( j == nysv  )  boundary_flags(j,i) = 7
+                IF ( j == nys   )  boundary_flags(j,i) = 8
+             ENDIF
+
+          ENDDO
+       ENDDO
+       
+    END SUBROUTINE ws_init
+
+
+!------------------------------------------------------------------------------!
+! Initialize variables used for storing statistic qauntities (fluxes, variances)
+!------------------------------------------------------------------------------!
+    SUBROUTINE ws_statistics
+    
+       USE control_parameters
+       USE statistics
+
+       IMPLICIT NONE
+
+!       
+!--    The arrays needed for statistical evaluation are set to to 0 at the 
+!--    begin of prognostic_equations.
+       IF ( ws_scheme_mom )  THEN
           sums_wsus_ws_l = 0.0
           sums_wsvs_ws_l = 0.0
@@ -122,98 +244,11 @@
 
        IF ( ws_scheme_sca )  THEN
-          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions))
           sums_wspts_ws_l = 0.0
-          IF ( humidity  .OR.  passive_scalar )  THEN
-             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions))
-             sums_wsqs_ws_l = 0.0
-          ENDIF
-          IF ( ocean )  THEN
-             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions))
-             sums_wssas_ws_l = 0.0
-          ENDIF
+          IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
+          IF ( ocean )  sums_wssas_ws_l = 0.0
+
        ENDIF
-!
-!--    Arrays needed for reasons of speed optimization for cache and noopt 
-!--    version. For the vector version the buffer arrays are not necessary, 
-!--    because the the fluxes can swapped directly inside the loops of the 
-!--    advection routines.
-       IF ( loop_optimization /= 'vector' )  THEN
-          IF ( ws_scheme_mom )  THEN
-             ALLOCATE(flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),            &
-                flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),                  &
-                diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt))
-             ALLOCATE(flux_l_u(nzb+1:nzt,nys:nyn),                         &
-                flux_l_v(nzb+1:nzt,nys:nyn), flux_l_w(nzb+1:nzt,nys:nyn),  & 
-                diss_l_u(nzb+1:nzt,nys:nyn), diss_l_v(nzb+1:nzt,nys:nyn),  & 
-                diss_l_w(nzb+1:nzt,nys:nyn)) 
-          ENDIF
-          IF ( ws_scheme_sca )  THEN
-             ALLOCATE(flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt),           &
-                diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt))
-             ALLOCATE(flux_l_pt(nzb+1:nzt,nys:nyn),                        &
-                flux_l_e(nzb+1:nzt,nys:nyn), diss_l_pt(nzb+1:nzt,nys:nyn), & 
-                diss_l_e(nzb+1:nzt,nys:nyn))
-             IF ( humidity  .OR.  passive_scalar )  THEN
-                ALLOCATE(flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt))
-                ALLOCATE(flux_l_q(nzb+1:nzt,nys:nyn),                      &
-                   diss_l_q(nzb+1:nzt,nys:nyn))
-             ENDIF
-             IF ( ocean )  THEN
-                ALLOCATE(flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt))
-                ALLOCATE(flux_l_sa(nzb+1:nzt,nys:nyn),                     &
-                   diss_l_sa(nzb+1:nzt,nys:nyn))
-             ENDIF
-          ENDIF
-       ENDIF
-!
-!--    Determine the flags where the order of the scheme for horizontal 
-!--    advection should be degraded.
-       ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             boundary_flags(j,i) = 0
-             IF ( outflow_l )  THEN
-                IF ( i == nxlu )  boundary_flags(j,i) = 5
-                IF ( i == nxl )  boundary_flags(j,i) = 6
-             ELSEIF ( outflow_r )  THEN
-                IF ( i == nxr-1 )  boundary_flags(j,i) = 1
-                IF ( i == nxr )  boundary_flags(j,i) = 2
-             ELSEIF ( outflow_n )  THEN
-                IF ( j == nyn-1 )  boundary_flags(j,i) = 3
-                IF ( j == nyn )  boundary_flags(j,i) = 4
-             ELSEIF ( outflow_s )  THEN
-                IF ( j == nysv )  boundary_flags(j,i) = 7
-                IF ( j == nys )  boundary_flags(j,i) = 8
-             ENDIF
-          ENDDO
-       ENDDO
-       
-    END SUBROUTINE ws_init
-    
-
-    
-    SUBROUTINE ws_statistics
-    
-       USE control_parameters
-       USE statistics
-
-!       
-!--      The arrays needed for statistical evaluation are set to to 0 at the 
-!--      begin of prognostic_equations.
-          IF ( ws_scheme_mom )  THEN
-             sums_wsus_ws_l = 0.0
-             sums_wsvs_ws_l = 0.0
-             sums_us2_ws_l = 0.0
-             sums_vs2_ws_l = 0.0
-             sums_ws2_ws_l = 0.0
-          ENDIF
-          IF ( ws_scheme_sca )  THEN
-             sums_wspts_ws_l = 0.0
-             IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0
-             IF ( ocean )  sums_wssas_ws_l = 0.0
-          ENDIF
 
     END SUBROUTINE ws_statistics
-
 
 
@@ -222,5 +257,6 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
-               swap_diss_y_local, swap_flux_x_local, swap_diss_x_local  )
+                              swap_diss_y_local, swap_flux_x_local, &
+                              swap_diss_x_local  )
 
        USE arrays_3d
@@ -239,10 +275,11 @@
        REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
                                              flux_n, diss_n
-       REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &       
+       REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &
                                              swap_diss_y_local
        REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
                                              swap_diss_x_local
        CHARACTER (LEN = *), INTENT(IN)    :: sk_char
-       
+
+
        degraded_l = .FALSE.
        degraded_s = .FALSE.
@@ -250,168 +287,180 @@
        IF ( boundary_flags(j,i) /= 0 )  THEN
 !       
-!--      Degrade the order for Dirichlet bc. at the outflow boundary.
-         SELECT CASE ( boundary_flags(j,i) )
-         
-            CASE ( 1 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  u_comp    = u(k,j,i+1) - u_gtrans
-                  flux_r(k) = u_comp * (                                      &
-                              7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
-                              -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
-                  diss_r(k) = - abs(u_comp) * (                               &
-                              3. * ( sk(k,j,i+1)-sk(k,j,i)     )              & 
-                              -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+!--       Degrade the order for Dirichlet bc. at the outflow boundary
+          SELECT CASE ( boundary_flags(j,i) )
+
+             CASE ( 1 )
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   u_comp    = u(k,j,i+1) - u_gtrans
+                   flux_r(k) = u_comp * (                                      &
+                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
+                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
+                   diss_r(k) = -ABS( u_comp ) * (                              &
+                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
+                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+                ENDDO
+
+             CASE ( 2 )
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   u_comp    = u(k,j,i+1) - u_gtrans
+                   flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
+                   diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
+                                         sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
+                                         0.5, ddx )
+                ENDDO
+
+             CASE ( 3 )
+
+                DO  k = nzb_s_inner(j,i) + 1, nzt
+                   v_comp = v(k,j+1,i) - v_gtrans
+                   flux_n(k) = v_comp * (                                      &
+                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
+                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
+                   diss_n(k) = -ABS( v_comp ) * (                              &
+                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)     )           &
+                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
                ENDDO
-               
-            CASE ( 2 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  u_comp    = u(k,j,i+1) - u_gtrans
-                  flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
-                  diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
-                                        sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
-                                        0.5, ddx )
-              ENDDO
-              
-            CASE ( 3 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  v_comp = v(k,j+1,i) - v_gtrans
-                  flux_n(k) = v_comp * (                                      &
-                              7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
-                              -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( sk(k,j+1,i)-sk(k,j,i)    )               & 
-                              -    (sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
-              ENDDO
-            CASE ( 4 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  v_comp    = v(k,j+1,i) - v_gtrans
-                  flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
-                  diss_n(k) =  diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
+
+             CASE ( 4 )
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   v_comp    = v(k,j+1,i) - v_gtrans
+                   flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
+                   diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
                                          sk(k,j-1,i), sk(k,j-2,i), v_comp,    &
                                          0.5, ddy )     
-               ENDDO
+                ENDDO
+
+             CASE ( 5 )
+
+                DO  k = nzb_w_inner(j,i)+1, nzt
+                   u_comp    = u(k,j,i+1) - u_gtrans 
+                   flux_r(k) = u_comp  * (                                     &
+                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
+                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
+                   diss_r(k) = -ABS( u_comp ) * (                              &
+                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
+                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+                ENDDO  
+
+             CASE ( 6 )
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   u_comp    = u(k,j,i+1) - u_gtrans
+                   flux_r(k) = u_comp * (                                      &
+                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
+                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
+                   diss_r(k) = -ABS( u_comp ) * (                              &
+                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
+                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+!                              
+!--                Compute leftside fluxes for the left boundary of PE domain
+                   u_comp                 = u(k,j,i) - u_gtrans
+                   swap_flux_x_local(k,j) = u_comp * (                         &
+                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
+                   swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
+                                                      sk(k,j,i), sk(k,j,i-1),  &
+                                                      sk(k,j,i-1), u_comp,     &
+                                                      0.5, ddx ) 
+                ENDDO
+                degraded_l = .TRUE.
+
+             CASE ( 7 )
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   v_comp    = v(k,j+1,i)-v_gtrans
+                   flux_n(k) = v_comp * (                                      &
+                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
+                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
+                   diss_n(k) = -ABS( v_comp ) * (                              &
+                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
+                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
+                ENDDO
+
+             CASE ( 8 )
+
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   v_comp    = v(k,j+1,i) - v_gtrans 
+                   flux_n(k) = v_comp * (                                      &
+                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
+                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
+                   diss_n(k) = -ABS( v_comp ) * (                              &
+                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
+                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
+!                              
+!--                Compute southside fluxes for the south boundary of PE domain
+                   v_comp               = v(k,j,i) - v_gtrans
+                   swap_flux_y_local(k) = v_comp *                             &
+                                          ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 
+                   swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
+                                                    sk(k,j,i), sk(k,j-1,i),    &
+                                                    sk(k,j-1,i), v_comp,       &
+                                                    0.5, ddy )
+                ENDDO
+                degraded_s = .TRUE.
                
-            CASE ( 5 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
-                  u_comp    = u(k,j,i+1) - u_gtrans 
-                  flux_r(k) = u_comp  * (                                     &
-                              7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
-                              -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
-                  diss_r(k) = - abs(u_comp) * (                               &
-                              3. * ( sk(k,j,i+1)-sk(k,j,i) )                  & 
-                              -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
-               ENDDO  
-               
-            CASE ( 6 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  u_comp    = u(k,j,i+1) - u_gtrans
-                  flux_r(k) = u_comp * (                                      &
-                              7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
-                              -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
-                  diss_r(k) = - abs(u_comp) * (                               &
-                              3. * ( sk(k,j,i+1)-sk(k,j,i)     )              & 
-                              -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
-!                              
-!--               Compute leftside fluxes for the left boundary of PE domain.
-                  u_comp                 = u(k,j,i) - u_gtrans
-                  swap_flux_x_local(k,j) = u_comp * (                         &
-                                           sk(k,j,i) + sk(k,j,i-1) ) * 0.5
-                  swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
-                                                     sk(k,j,i), sk(k,j,i-1),  &
-                                                     sk(k,j,i-1), u_comp,     &
-                                                     0.5, ddx ) 
-               ENDDO
-               degraded_l = .TRUE.
-               
-            CASE ( 7 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  v_comp    = v(k,j+1,i)-v_gtrans
-                  flux_n(k) = v_comp * (                                      &
-                              7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
-                              -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( sk(k,j+1,i) - sk(k,j,i)   )              & 
-                              -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
-               ENDDO
-               
-            CASE ( 8 )
-               DO  k = nzb_s_inner(j,i) + 1, nzt
-                  v_comp    = v(k,j+1,i) - v_gtrans 
-                  flux_n(k) = v_comp * (                                      &
-                              7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
-                              -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( sk(k,j+1,i) - sk(k,j,i) )                & 
-                              -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
-!                              
-!--               Compute southside fluxes for the south boundary of PE domain           
-                  v_comp               = v(k,j,i) - v_gtrans
-                  swap_flux_y_local(k) = v_comp * (                           &
-                                         sk(k,j,i)+ sk(k,j-1,i) ) * 0.5 
-                  swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
-                                                   sk(k,j,i), sk(k,j-1,i),    &
-                                                   sk(k,j-1,i), v_comp,       &
-                                                   0.5, ddy )
-               ENDDO
-               degraded_s = .TRUE.
-               
-            CASE DEFAULT
+             CASE DEFAULT
             
-         END SELECT
+          END SELECT
+
 !         
-!--      Compute the crosswise 5th order fluxes at the outflow
-         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
-         .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
-            DO  k = nzb_s_inner(j,i) + 1, nzt
-               v_comp    = v(k,j+1,i) - v_gtrans
-               flux_n(k) = v_comp * (                                         &
-                           37.  * ( sk(k,j+1,i) + sk(k,j,i)   )               &
-                           - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
-                           +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
-               diss_n(k) = - abs(v_comp) * (                                  &
-                           10. *  ( sk(k,j+1,i) - sk(k,j,i)   )               &
-                           - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
-                           +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
-            ENDDO
+!--       Compute the crosswise 5th order fluxes at the outflow
+          IF ( boundary_flags(j,i) == 1  .OR.  boundary_flags(j,i) == 2  .OR. &
+               boundary_flags(j,i) == 5  .OR.  boundary_flags(j,i) == 6 )  THEN
+
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                v_comp    = v(k,j+1,i) - v_gtrans
+                flux_n(k) = v_comp * (                                         &
+                            37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )               &
+                          -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
+                          +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
+                diss_n(k) = -ABS( v_comp ) * (                                 &
+                            10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )               &
+                          -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
+                          +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
+             ENDDO
             
-         ELSE
+          ELSE
          
-            DO  k = nzb_s_inner(j,i) + 1, nzt
-               u_comp    = u(k,j,i+1) - u_gtrans  
-               flux_r(k) = u_comp * (                                         &
-                           37.  * ( sk(k,j,i+1) + sk(k,j,i)   )               &
-                           - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
-                           +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
-               diss_r(k) = - abs(u_comp) * (                                  &
-                           10.  * ( sk(k,j,i+1) - sk(k,j,i)   )               &
-                           - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
-                           +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
-            ENDDO
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                u_comp    = u(k,j,i+1) - u_gtrans  
+                flux_r(k) = u_comp * (                                         &
+                            37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )               &
+                          -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
+                          +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
+                diss_r(k) = -ABS( u_comp ) * (                                 &
+                            10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )               &
+                          -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
+                          +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
+             ENDDO
             
-         ENDIF
-!
-!--    Compute the fifth order fluxes for the interior of PE domain.
+          ENDIF
+
        ELSE
                
-          DO  k = nzb_u_inner(j,i) + 1, nzt
+!
+!--       Compute the fifth order fluxes for the interior of PE domain.
+          DO  k = nzb_u_inner(j,i)+1, nzt
              u_comp    = u(k,j,i+1) - u_gtrans
              flux_r(k) = u_comp * (                                           &
-                         37.  * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
-                         - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
-                         +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
-             diss_r(k) = - abs(u_comp) * (                                    &
-                         10.  * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
-                         - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
-                         +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
+                         37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
+                       -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
+                       +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
+             diss_r(k) = -ABS( u_comp ) * (                                   &
+                         10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
+                       -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
+                       +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
 
              v_comp    = v(k,j+1,i) - v_gtrans
              flux_n(k) = v_comp * (                                           &
-                         37.  * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
-                         - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
-                         +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
-             diss_n(k) = - abs(v_comp) * (                                    &
-                         10.  * ( sk(k,j+1,i) - sk(k,j,i)   )                 &                      
-                         - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
-                         +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5               
+                         37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
+                       -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
+                       +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
+             diss_n(k) = -ABS( v_comp ) * (                                   &
+                         10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
+                       -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
+                       +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
           ENDDO
           
@@ -424,43 +473,45 @@
               v_comp               = v(k,j,i) - v_gtrans
               swap_flux_y_local(k) = v_comp * (                               &
-                                     37. *    ( sk(k,j,i) + sk(k,j-1,i)   )   &
-                                     - 8. *   ( sk(k,j+1,i) + sk(k,j-2,i) )   &
-                                     +        ( sk(k,j+2,i) + sk(k,j-3,i) )   &
-                                     ) * adv_sca_5
-              swap_diss_y_local(k) = - abs(v_comp) * (                        &
-                                     10. *  ( sk(k,j,i) - sk(k,j-1,i)   )     &
-                                     - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
-                                     +        sk(k,j+2,i) - sk(k,j-3,i)       &
-                                     ) * adv_sca_5
+                                     37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
+                                   -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
+                                   +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
+                                              ) * adv_sca_5
+              swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
+                                     10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
+                                   -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
+                                   +          sk(k,j+2,i) - sk(k,j-3,i)       &
+                                                      ) * adv_sca_5
            ENDDO           
           
         ENDIF
         
-        IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
+        IF ( i == nxl  .AND.  .NOT. degraded_l )  THEN
         
-           DO  k = nzb_s_inner(j,i) + 1, nzt
+           DO  k = nzb_s_inner(j,i)+1, nzt
               u_comp                 = u(k,j,i) - u_gtrans
               swap_flux_x_local(k,j) = u_comp * (                             &
-                                       37.  * ( sk(k,j,i) + sk(k,j,i-1)   )   &
-                                       - 8. * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
-                                       +      ( sk(k,j,i+2) + sk(k,j,i-3) )   &
-                                       ) * adv_sca_5
-              swap_diss_x_local(k,j) = - abs(u_comp) * (                      &
-                                       10.  * ( sk(k,j,i) - sk(k,j,i-1)   )   &
-                                       - 5. * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
-                                       +      ( sk(k,j,i+2) - sk(k,j,i-3) )   &
-                                       ) * adv_sca_5
+                                       37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
+                                     -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
+                                     +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
+                                                ) * adv_sca_5
+              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
+                                       10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
+                                     -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
+                                     +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
+                                                        ) * adv_sca_5
            ENDDO
            
         ENDIF
+
 !        
-!--    Now compute the tendency terms for the horizontal parts.
-       DO  k = nzb_s_inner(j,i) + 1, nzt
+!--    Now compute the tendency terms for the horizontal parts
+       DO  k = nzb_s_inner(j,i)+1, nzt
           
-          tend(k,j,i) = tend(k,j,i) - (                                       &
-               ( flux_r(k) + diss_r(k)                                        &
-             - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx    &
-             + (flux_n(k)+diss_n(k)                                           &
-             - (swap_flux_y_local(k)    + swap_diss_y_local(k) ) )   * ddy )
+          tend(k,j,i) = tend(k,j,i) - (                                      &
+                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
+                            swap_diss_x_local(k,j) ) * ddx                   &
+                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   - &
+                            swap_diss_y_local(k)   ) * ddy                   &
+                                      )
 
           swap_flux_y_local(k)   = flux_n(k)
@@ -470,107 +521,123 @@
           
        ENDDO
-!
-!--    Vertical advection, degradation of order near surface and top. 
-!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
-!--    statistical evaluation the top flux at the surface should be 0
-
-       flux_t(nzb_s_inner(j,i)) = 0.
-       diss_t(nzb_s_inner(j,i)) = 0.
+
+!
+!--    Vertical advection, degradation of order near bottom and top. 
+!--    The fluxes flux_d and diss_d at the surface are 0. Due to later
+!--    calculation of statistics the top flux at the surface should be 0.
+       flux_t(nzb_s_inner(j,i)) = 0.0
+       diss_t(nzb_s_inner(j,i)) = 0.0
+
 !       
-!--    2nd order scheme       
+!--    2nd-order scheme (bottom)
        k         = nzb_s_inner(j,i) + 1
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
        flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
+
 !       
-!--    sk(k,j,i) is referenced three times to avoid a access below surface.
+!--    sk(k,j,i) is referenced three times to avoid an access below surface
        diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i),  &
                              sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
                    
-       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!--    WS3 as an intermediate step       
+       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
+                                   * ddzw(k)
+!
+!--    WS3 as an intermediate step (bottom)
        k         = nzb_s_inner(j,i) + 2
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
        flux_t(k) = w(k,j,i) * (                                               &
-                   7. * ( sk(k+1,j,i) + sk(k,j,i)   )                         &
-                   -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
-       diss_t(k) = - abs(w(k,j,i)) * (                                        &
-                   3. * ( sk(k+1,j,i) - sk(k,j,i)   )                         & 
-                   -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
-
-       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d ) ) * ddzw(k)
+                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )           &
+                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )           &
+                              ) * adv_sca_3
+       diss_t(k) = -ABS( w(k,j,i) ) * (                                       &
+                                3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )           & 
+                                    - ( sk(k+2,j,i) - sk(k-1,j,i) )           &
+                                      ) * adv_sca_3
+
+       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
+                                   * ddzw(k)
 !
 !--    WS5
-       DO  k = nzb_s_inner(j,i) + 3, nzt - 2
+       DO  k = nzb_s_inner(j,i)+3, nzt-2
        
           flux_d    = flux_t(k-1)
-          diss_d    = diss_t(k-1)
-          
-          flux_t(k) = w(k,j,i) * (                                            &  
-                      37.  * ( sk(k+1,j,i) + sk(k,j,i)   )                    &
-                      - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )                    &
-                      +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) *adv_sca_5
-          diss_t(k) = - abs(w(k,j,i)) * (                                     & 
-                      10.  * ( sk(k+1,j,i) - sk(k,j,i)   )                    &      
-                      - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )                    &
-                      +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
+          diss_d    = diss_t(k-1)          
+          flux_t(k) = w(k,j,i) * (                                            &
+                                   37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )       &
+                                 -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )       &
+                                 +        ( sk(k+3,j,i) + sk(k-2,j,i) )       &
+                                 ) * adv_sca_5
+          diss_t(k) = -ABS( w(k,j,i) ) * (                                     &
+                                           10.0 * ( sk(k+1,j,i) - sk(k,j,i)   )&
+                                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )&
+                                         +        ( sk(k+3,j,i) - sk(k-2,j,i) )&
+                                         ) * adv_sca_5
 
           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
                                     - ( flux_d + diss_d ) ) * ddzw(k)
+
        ENDDO
-!
-!--    WS3 as an intermediate step 
+
+!
+!--    WS3 as an intermediate step (top)
        k         = nzt - 1
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
-       flux_t(k) = w(k,j,i) * (                                               &
-                   7. * ( sk(k+1,j,i) + sk(k,j,i) )                           &
-                   -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
-       diss_t(k) = - abs(w(k,j,i)) * (                                        &
-                   3. * ( sk(k+1,j,i) - sk(k,j,i) )                           & 
-                   -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
-
-       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d ) ) * ddzw(k)
+       flux_t(k) = w(k,j,i) * (                                                &
+                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )            &
+                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )            &
+                              ) * adv_sca_3
+       diss_t(k) = -ABS( w(k,j,i) ) * (                                        &
+                                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )    &
+                                            - ( sk(k+2,j,i) - sk(k-1,j,i) )    &
+                                      ) * adv_sca_3
+
+       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
+                                   * ddzw(k)
 !                                 
-!--    2nd
+!--    2nd-order scheme (top)
        k         = nzt 
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
-       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) *0.5
+       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
+
 !       
-!--    sk(k+1) is referenced two times to avoid a segmentation fault at top.
-       diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i),&
+!--    sk(k+1) is referenced two times to avoid a segmentation fault at top
+       diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i), &
                              sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
 
-       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d ) ) * ddzw(k)
+       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
+                                   * ddzw(k)
 !       
-!--    evaluation of statistics
+!--    Evaluation of statistics
        SELECT CASE ( sk_char )
 
           CASE ( 'pt' )
-            DO  k = nzb_s_inner(j,i), nzt
-               sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                     &
-                 + (flux_t(k)+diss_t(k))                                       &
-                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+
+             DO  k = nzb_s_inner(j,i), nzt
+                sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) +                  &
+                                       ( flux_t(k) + diss_t(k) )               &
+                                 * weight_substep(intermediate_timestep_count) &
+                                       * rmask(j,i,:)
              ENDDO
              
           CASE ( 'sa' )
-            DO  k = nzb_s_inner(j,i), nzt
-               sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                     &
-                 +(flux_t(k)+diss_t(k))                                        &
-                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+
+             DO  k = nzb_s_inner(j,i), nzt
+                sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) +                  &
+                                       ( flux_t(k) + diss_t(k) )               &
+                                 * weight_substep(intermediate_timestep_count) &
+                                       * rmask(j,i,:)
              ENDDO
              
           CASE ( 'q' )
+
              DO  k = nzb_s_inner(j,i), nzt
-                sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                      &
-                  + (flux_t(k)+diss_t(k))                                      &
-                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+                sums_wsqs_ws_l(k,:)  = sums_wsqs_ws_l(k,:) +                   &
+                                      ( flux_t(k) + diss_t(k) )                &
+                                 * weight_substep(intermediate_timestep_count) &
+                                      * rmask(j,i,:)
              ENDDO
 
Index: palm/trunk/SOURCE/exchange_horiz.f90
===================================================================
--- palm/trunk/SOURCE/exchange_horiz.f90	(revision 708)
+++ palm/trunk/SOURCE/exchange_horiz.f90	(revision 709)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! formatting adjustments
 !
 ! Former revisions:
@@ -69,5 +69,5 @@
 
 !
-!-- Exchange of lateral boundary values for parallel computers
+!-- Exchange of lateral boundary values
     IF ( pdims(1) == 1  .OR.  mg_switch_to_pe0 )  THEN
 !
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 708)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 709)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
@@ -114,4 +115,5 @@
     USE arrays_3d
     USE cloud_parameters
+    USE control_parameters
     USE cpulog
     USE grid_variables
@@ -120,5 +122,4 @@
     USE pegrid
     USE statistics
-    USE control_parameters
 
     IMPLICIT NONE
@@ -738,9 +739,8 @@
           ENDDO
        ENDDO
-!-     for reasons of speed optimization the loop is splitted, to avoid if-else
-!-     statements inside the loop
-!-     Fluxes which have been computed in part directly inside the advection routines
-!-     treated seperatly.
-!-     First treat the momentum fluxes
+!
+!--    For speed optimization fluxes which have been computed in part directly
+!--    inside the WS advection routines are treated seperatly
+!--    Momentum fluxes first:
        IF ( .NOT. ws_scheme_mom )  THEN
          !$OMP DO
@@ -771,6 +771,7 @@
          DO  i = nxl, nxr
             DO  j = nys, nyn
-               DO  k = nzb_diff_s_inner(j,i) - 1, nzt_diff
-!-                vertical heat flux
+               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
+!
+!--               Vertical heat flux
                   sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
                            ( pt(k,j,i)   - hom(k,1,4,sr) + &
@@ -1221,5 +1222,2 @@
 
  END SUBROUTINE flow_statistics
-
-
-
Index: palm/trunk/SOURCE/inflow_turbulence.f90
===================================================================
--- palm/trunk/SOURCE/inflow_turbulence.f90	(revision 708)
+++ palm/trunk/SOURCE/inflow_turbulence.f90	(revision 709)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
@@ -54,5 +55,4 @@
 !
 !-- First, local averaging within the recycling domain
-
     i = recycling_plane
 
@@ -79,6 +79,6 @@
 !-- Now, averaging over all PEs
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr,  & 
-                    MPI_REAL, MPI_SUM, comm2d, ierr )
+    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr, MPI_REAL, &
+                        MPI_SUM, comm2d, ierr )
 
 #else
@@ -165,9 +165,10 @@
           DO  k = nzb, nzt + 1
 
-              u(k,j,-nbgp+1:0)   = mean_inflow_profiles(k,1) + &
+              u(k,j,-nbgp+1:0) = mean_inflow_profiles(k,1) + &
                            inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
               v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
                            inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
-              w(k,j,-nbgp:-1)  = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
+              w(k,j,-nbgp:-1)  =                             &
+                           inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
               pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
                            inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
@@ -181,40 +182,4 @@
     ENDIF
 
-!
-!-- Conserve the volume flow at the inflow in order to avoid generation of
-!-- waves in the stable layer
-!    IF ( conserve_volume_flow  .AND.  inflow_l )  THEN
-
-!       volume_flow(1)   = 0.0
-!       volume_flow_l(1) = 0.0
-
-!       i = 0
-
-!       DO  j = nys, nyn
-!
-!--       Sum up the volume flow through the south/north boundary
-!          DO  k = nzb_2d(j,i) + 1, nzt
-!             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
-!          ENDDO
-!       ENDDO
-
-#if defined( __parallel )   
-!       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-!       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
-!                           MPI_SUM, comm1dy, ierr )   
-#else
-!       volume_flow = volume_flow_l  
-#endif
-!       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
-!                               / volume_flow_area(1)
-
-!       DO  j = nys-1, nyn+1
-!          DO  k = nzb_v_inner(j,i) + 1, nzt
-!             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
-!          ENDDO
-!       ENDDO
-
-!    ENDIF
-
     CALL cpu_log( log_point(40), 'inflow_turbulence', 'stop' )
 
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 708)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 709)
@@ -7,5 +7,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! formatting adjustments
 !
 ! Former revisions:
@@ -29,5 +29,5 @@
 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and 
 ! allocation of arrays. Calls of exchange_horiz are modified.
-! Call ws_init to initialize arrays needed for statistical evaluation and 
+! Call ws_init to initialize arrays needed for calculating statisticas and for
 ! optimization when ws-scheme is used.
 ! Initial volume flow is now calculated by using the variable hom_sum.
@@ -35,6 +35,6 @@
 ! removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc) 
 ! Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
-! mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is 
-! representative for the height z0 
+! mirror to Dirichlet boundary conditions (u=v=0), so that k=nzb is 
+! representative for the height z0.
 ! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner) 
 !
@@ -453,13 +453,14 @@
     
 !
-!-- Allocate arrays containing the RK coefficient for right evaluation of 
-!-- perturbation pressure and turbulent fluxes. At this point it is needed
-!-- for right pressure correction during initialization. Further below 
-!-- the real values will be set.
-    ALLOCATE (weight_substep(1:intermediate_timestep_count_max),           &
+!-- Allocate arrays containing the RK coefficient for calculation of 
+!-- perturbation pressure and turbulent fluxes. At this point values are
+!-- set for pressure calculation during initialization (where no timestep
+!-- is done). Further below the values needed within the timestep scheme
+!-- will be set.
+    ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
               weight_pres(1:intermediate_timestep_count_max) )
-    weight_substep = 1.
-    weight_pres = 1.
-    intermediate_timestep_count = 1 ! needed for simulated_time=0
+    weight_substep = 1.0
+    weight_pres    = 1.0
+    intermediate_timestep_count = 1  ! needed when simulated_time = 0.0
        
 !
@@ -1011,5 +1012,4 @@
 !
 !--    Read binary data from restart file
-
        CALL read_3d_binary
 
@@ -1119,5 +1119,5 @@
 !
 !-- Calculate the initial volume flow at the right and north boundary
-    IF ( conserve_volume_flow ) THEN
+    IF ( conserve_volume_flow )  THEN
 
        volume_flow_initial_l = 0.0
@@ -1128,5 +1128,5 @@
           IF ( nxr == nx )  THEN
              DO  j = nys, nyn
-                DO  k = nzb_2d(j,nx) + 1, nzt
+                DO  k = nzb_2d(j,nx)+1, nzt
                    volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
                                               hom_sum(k,1,0) * dzw(k)
@@ -1138,7 +1138,7 @@
           IF ( nyn == ny )  THEN
              DO  i = nxl, nxr
-                DO  k = nzb_2d(ny,i) + 1, nzt  
+                DO  k = nzb_2d(ny,i)+1, nzt  
                    volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
-                                               hom_sum(k,2,0) * dzw(k)
+                                              hom_sum(k,2,0) * dzw(k)
                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
                 ENDDO
@@ -1150,7 +1150,7 @@
           IF ( nxr == nx )  THEN
              DO  j = nys, nyn
-                DO  k = nzb_2d(j,nx) + 1, nzt
+                DO  k = nzb_2d(j,nx)+1, nzt
                    volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
-                                               u(k,j,nx) * dzw(k)
+                                              u(k,j,nx) * dzw(k)
                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
                 ENDDO
@@ -1160,5 +1160,5 @@
           IF ( nyn == ny )  THEN
              DO  i = nxl, nxr
-                DO  k = nzb_2d(ny,i) + 1, nzt  
+                DO  k = nzb_2d(ny,i)+1, nzt  
                    volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
                                               v(k,ny,i) * dzw(k)
@@ -1171,7 +1171,7 @@
 
 #if defined( __parallel )
-       CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
+       CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), &
                            2, MPI_REAL, MPI_SUM, comm2d, ierr )
-       CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
+       CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),       &
                            2, MPI_REAL, MPI_SUM, comm2d, ierr )
 
@@ -1182,6 +1182,6 @@
 
 !
-!--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
-!--       and calculated from u|v_bulk instead.
+!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
+!--    from u|v_bulk instead
        IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
           volume_flow_initial(1) = u_bulk * volume_flow_area(1)
@@ -1208,5 +1208,5 @@
 !
 !-- Initialization of the leaf area density
-     IF ( plant_canopy ) THEN
+    IF ( plant_canopy )  THEN
  
        SELECT CASE ( TRIM( canopy_mode ) )
@@ -1218,5 +1218,5 @@
                    lad_s(:,j,i) = lad(:)
                    cdc(:,j,i)   = drag_coefficient
-                   IF ( passive_scalar ) THEN
+                   IF ( passive_scalar )  THEN
                       sls(:,j,i) = leaf_surface_concentration
                       sec(:,j,i) = scalar_exchange_coefficient
@@ -1240,5 +1240,5 @@
        CALL exchange_horiz( cdc, nbgp )
 
-       IF ( passive_scalar ) THEN
+       IF ( passive_scalar )  THEN
           CALL exchange_horiz( sls, nbgp )
           CALL exchange_horiz( sec, nbgp )
@@ -1255,5 +1255,5 @@
           DO  j = nys, nyn
              DO  k = nzb, nzt+1 
-                IF ( lad_s(k,j,i) > 0.0 ) THEN
+                IF ( lad_s(k,j,i) > 0.0 )  THEN
                    lad_u(k,j,i)   = lad_s(k,j,i) 
                    lad_u(k,j,i+1) = lad_s(k,j,i)
@@ -1277,5 +1277,5 @@
 !
 !--    Initialisation of the canopy heat source distribution
-       IF ( cthf /= 0.0 ) THEN
+       IF ( cthf /= 0.0 )  THEN
 !
 !--       Piecewise evaluation of the leaf area index by 
@@ -1313,5 +1313,5 @@
           shf(:,:) = canopy_heat_flux(0,:,:)
 
-          IF ( ASSOCIATED( shf_m ) ) shf_m = shf
+          IF ( ASSOCIATED( shf_m ) )  shf_m = shf
 
        ENDIF
@@ -1346,23 +1346,29 @@
 
 !
-!-- Setting weighting factors for right evaluation of perturbation pressure 
-!-- and turbulent quantities during the RK substeps.                
-    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN    ! RK3       
+!-- Setting weighting factors for calculation of perturbation pressure 
+!-- and turbulent quantities from the RK substeps                
+    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
+
        weight_substep(1) = 0.166666666666666
        weight_substep(2) = 0.3
        weight_substep(3) = 0.533333333333333
-          
-       weight_pres(1) = 0.333333333333333
-       weight_pres(2) = 0.416666666666666
-       weight_pres(3) = 0.25          
-    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! RK2       
+
+       weight_pres(1)    = 0.333333333333333
+       weight_pres(2)    = 0.416666666666666
+       weight_pres(3)    = 0.25
+
+    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
+
        weight_substep(1) = 0.5
        weight_substep(2) = 0.5
           
-       weight_pres(1) = 0.5
-       weight_pres(2) = 0.5          
-    ELSE                                                  ! Euler and Leapfrog       
+       weight_pres(1)    = 0.5
+       weight_pres(2)    = 0.5          
+
+    ELSE                                     ! for Euler- and leapfrog-method
+
        weight_substep(1) = 1.0      
-       weight_pres(1) = 1.0                   
+       weight_pres(1)    = 1.0                   
+
     ENDIF
 
@@ -1476,8 +1482,8 @@
 
 !
-!-- Initialize local summation arrays for UP flow_statistics. This is necessary
-!-- because they may not yet have been initialized when they are called from 
-!-- flow_statistics (or - depending on the chosen model run - are never 
-!-- initialized)
+!-- Initialize local summation arrays for routine flow_statistics.
+!-- This is necessary because they may not yet have been initialized when they
+!-- are called from flow_statistics (or - depending on the chosen model run -
+!-- are never initialized)
     sums_divnew_l      = 0.0
     sums_divold_l      = 0.0
@@ -1492,5 +1498,5 @@
 !
 !-- User-defined initializing actions. Check afterwards, if maximum number
-!-- of allowed timeseries is not exceeded
+!-- of allowed timeseries is exceeded
     CALL user_init
 
Index: palm/trunk/SOURCE/init_coupling.f90
===================================================================
--- palm/trunk/SOURCE/init_coupling.f90	(revision 708)
+++ palm/trunk/SOURCE/init_coupling.f90	(revision 709)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 708)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 709)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! formatting adjustments
 !
 ! Former revisions:
@@ -89,9 +89,10 @@
     
 !
-!   Computation of the array bounds.
+!-- Calculation of horizontal array bounds including ghost layers
     nxlg = nxl - nbgp
     nxrg = nxr + nbgp
     nysg = nys - nbgp
     nyng = nyn + nbgp
+
 !
 !-- Allocate grid arrays
@@ -196,13 +197,13 @@
     
 !    
-!-- In case of FFT method or SOR swap inverse grid lenght ddzu to ddzu_fft and 
-!-- modify the lowest entry. This is necessary for atmosphere runs, because
-!-- zu(0) and so ddzu(1) changed. Accompanied with this modified arrays 
-!-- pressure solver uses wrong values and this causes kinks in the profiles
-!-- of turbulent quantities.  
-    IF ( psolver /= 'multigrid' ) THEN
+!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
+!-- everywhere. For the actual grid, the grid spacing at the lowest level
+!-- is only dz/2, but should be dz. Therefore, an additional array
+!-- containing with appropriate grid information is created for these
+!-- solvers.
+    IF ( psolver /= 'multigrid' )  THEN
        ALLOCATE( ddzu_pres(1:nzt+1) )
        ddzu_pres = ddzu
-       IF( .NOT. ocean ) ddzu_pres(1) = ddzu_pres(2)
+       IF( .NOT. ocean )  ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
     ENDIF   
 
@@ -221,8 +222,10 @@
        dzu_mg(:,maximum_grid_level) = dzu
 !       
-!--    To ensure a equally spaced grid. For ocean runs this is not necessary, 
+!--    Next line to ensure an equally spaced grid. For ocean runs this is not
+!--    necessary, 
 !--    because zu(0) does not changed so far. Also this would cause errors 
 !--    if a vertical stretching for ocean runs is used.
        IF ( .NOT. ocean ) dzu_mg(1,maximum_grid_level) = dzu(2)
+
        dzw_mg(:,maximum_grid_level) = dzw
        nzt_l = nzt
@@ -1132,6 +1135,5 @@
 
 !
-!-- Need to set lateral boundary conditions for l_wall
-
+!-- Set lateral boundary conditions for l_wall
     CALL exchange_horiz( l_wall, nbgp )
 
Index: palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- palm/trunk/SOURCE/init_pegrid.f90	(revision 708)
+++ palm/trunk/SOURCE/init_pegrid.f90	(revision 709)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! formatting adjustments
 ! 
 ! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
@@ -577,17 +577,7 @@
 
           CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
-!
-!--       TEST OUTPUT (TO BE REMOVED)
-          WRITE(9,*)  TRIM( coupling_mode ),  &
-               ', ierr after MPI_OPEN_PORT: ', ierr
-          CALL LOCAL_FLUSH( 9 )
 
           CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
                                  ierr )
-!
-!--       TEST OUTPUT (TO BE REMOVED)
-          WRITE(9,*)  TRIM( coupling_mode ),  &
-               ', ierr after MPI_PUBLISH_NAME: ', ierr
-          CALL LOCAL_FLUSH( 9 )
 
 !
@@ -614,10 +604,4 @@
 
           CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
-!
-!--       TEST OUTPUT (TO BE REMOVED)
-          WRITE(9,*)  TRIM( coupling_mode ),  &
-               ', ierr after MPI_LOOKUP_NAME: ', ierr
-          CALL LOCAL_FLUSH( 9 )
-
 
        ENDIF
@@ -631,21 +615,12 @@
     IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
 
-       PRINT*, '... before COMM_ACCEPT'
        CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
                              comm_inter, ierr )
-       PRINT*, '--- ierr = ', ierr
-       PRINT*, '--- comm_inter atmosphere = ', comm_inter
-
        coupling_mode_remote = 'ocean_to_atmosphere'
 
     ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
 
-       IF ( myid == 0 )  PRINT*, '*** read: ', port_name, '  ierr = ', ierr
-       PRINT*, '... before COMM_CONNECT'
        CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
                               comm_inter, ierr )
-       PRINT*, '--- ierr = ', ierr
-       PRINT*, '--- comm_inter ocean      = ', comm_inter
-
        coupling_mode_remote = 'atmosphere_to_ocean'
 
@@ -654,26 +629,20 @@
 
 ! 
-!-- Determine the number of ghost points
-    IF (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme') THEN
+!-- Determine the number of ghost point layers
+    IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' )  THEN
        nbgp = 3
     ELSE
        nbgp = 1
-    END IF 
-
-!
-!-- In case of coupled runs, create a new MPI derived datatype for the
-!-- exchange of surface (xy) data .
-!-- Gridpoint number for the exchange of ghost points (xy-plane)
-
+    ENDIF 
+
+!
+!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
+!-- which is needed for coupled atmosphere-ocean runs.
+!-- First, calculate number of grid points of an xy-plane.
     ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
-
-!
-!-- Define a new MPI derived datatype for the exchange of ghost points in
-!-- y-direction for 2D-arrays (line)
     CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
     CALL MPI_TYPE_COMMIT( type_xy, ierr )
 
-
-    IF ( TRIM( coupling_mode ) .NE. 'uncoupled' ) THEN
+    IF ( TRIM( coupling_mode ) /= 'uncoupled' )  THEN
    
 !
@@ -685,22 +654,23 @@
           ny_a = ny
 
-          IF ( myid == 0 ) THEN
-             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, &
-                            comm_inter, ierr )
-             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, &
-                            comm_inter, ierr )
-             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, &
-                            comm_inter, ierr )
-             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, &
-                            comm_inter, status, ierr )
-             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, &
-                            comm_inter, status, ierr )
-             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
+          IF ( myid == 0 )  THEN
+
+             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter,  &
+                            ierr )
+             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter,  &
+                            ierr )
+             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
+                            ierr )
+             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter,  &
+                            status, ierr )
+             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter,  &
+                            status, ierr )
+             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6,      &
                             comm_inter, status, ierr )
           ENDIF
 
-          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr)
-          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr) 
-          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
+          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
+          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr ) 
+          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
        
        ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
@@ -710,16 +680,14 @@
 
           IF ( myid == 0 ) THEN
-             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, &
-                            comm_inter, status, ierr )
-             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, &
-                            comm_inter, status, ierr )
-             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, &
-                            comm_inter, status, ierr )
-             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, &
-                            comm_inter, ierr )
-             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, &
-                            comm_inter, ierr )
-             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, &
-                            comm_inter, ierr )
+
+             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
+                            ierr )
+             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
+                            ierr )
+             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
+                            status, ierr )
+             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
+             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
+             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
           ENDIF
 
@@ -730,12 +698,11 @@
        ENDIF
   
-       ngp_a = (nx_a+1+2*nbgp)*(ny_a+1+2*nbgp)
-       ngp_o = (nx_o+1+2*nbgp)*(ny_o+1+2*nbgp)
-
-!
-!--    determine if the horizontal grid and the number of PEs 
-!--    in ocean and atmosphere is same or not
-!--    (different number of PEs still not implemented)
-       IF ( nx_o == nx_a .AND. ny_o == ny_a .AND.  &
+       ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
+       ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
+
+!
+!--    Determine if the horizontal grid and the number of PEs in ocean and
+!--    atmosphere is same or not
+       IF ( nx_o == nx_a  .AND.  ny_o == ny_a  .AND.  &
             pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
        THEN
@@ -748,6 +715,9 @@
 !--    Determine the target PEs for the exchange between ocean and
 !--    atmosphere (comm2d)
-       IF ( coupling_topology == 0) THEN
-          IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
+       IF ( coupling_topology == 0 )  THEN
+!
+!--       In case of identical topologies, every atmosphere PE has exactly one
+!--       ocean PE counterpart and vice versa
+          IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
              target_id = myid + numprocs
           ELSE
@@ -756,18 +726,19 @@
 
        ELSE
-
 !
 !--       In case of nonequivalent topology in ocean and atmosphere only for
 !--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
-!--       data echxchange between ocean and atmosphere will be done only by
-!--       those PEs.    
-          IF ( myid == 0 ) THEN
-             IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
+!--       data echxchange between ocean and atmosphere will be done only
+!--       between these PEs.    
+          IF ( myid == 0 )  THEN
+
+             IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' )  THEN
                 target_id = numprocs 
              ELSE
                 target_id = 0
              ENDIF
- print*, coupling_mode, myid, " -> ", target_id, "numprocs: ", numprocs 
+
           ENDIF
+
        ENDIF
 
@@ -861,7 +832,9 @@
 !
 !--    Find out, if the total domain allows more levels. These additional
-!--    levels are processed on PE0 only.
+!--    levels are identically processed on all PEs.
        IF ( numprocs > 1  .AND.  mg_switch_to_pe0_level /= -1 )  THEN
+
           IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
+
              mg_switch_to_pe0_level_l = maximum_grid_level
 
@@ -889,7 +862,10 @@
                 mg_switch_to_pe0_level_l = 0
              ENDIF
+
           ELSE
+
              mg_switch_to_pe0_level_l = 0
              maximum_grid_level_l = maximum_grid_level
+
           ENDIF
 
@@ -920,4 +896,5 @@
 
              ENDIF
+
           ENDIF
 
@@ -939,5 +916,4 @@
 !--          Save the grid size of the subdomain at the switch level, because
 !--          it is needed in poismg.
-!--          Array bounds of the local subdomain grids are gathered on PE0
              ind(1) = nxl_l; ind(2) = nxr_l
              ind(3) = nys_l; ind(4) = nyn_l
@@ -953,5 +929,5 @@
              DEALLOCATE( ind_all )
 !
-!--          Calculate the grid size of the total domain gathered on PE0
+!--          Calculate the grid size of the total domain
              nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
              nxl_l = 0
@@ -1006,9 +982,11 @@
 
 !
-!-- Define a new MPI derived datatype for the exchange of ghost points in
-!-- y-direction for 2D-arrays (line)
-    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, ierr )
+!-- Define new MPI derived datatypes for the exchange of ghost points in
+!-- x- and y-direction for 2D-arrays (line)
+    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
+                          ierr )
     CALL MPI_TYPE_COMMIT( type_x, ierr )
-    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, type_x_int, ierr )
+    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
+                          type_x_int, ierr )
     CALL MPI_TYPE_COMMIT( type_x_int, ierr )
 
@@ -1029,37 +1007,43 @@
 
     nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
+
 !
 !-- Discern between the model grid, which needs nbgp ghost points and
 !-- grid levels for the multigrid scheme. In the latter case only one
 !-- ghost point is necessary.
-!-- First definition of mpi-vectors for exchange of ghost layers on normal 
+!-- First definition of MPI-datatypes for exchange of ghost layers on normal 
 !-- grid. The following loop is needed for data exchange in poismg.f90.
 !
 !-- Determine number of grid points of yz-layer for exchange
     ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
-!
-!-- Define a new mpi datatype for the exchange of left - right boundaries.
-!-- Indeed the data are connected in the physical memory and no mpi-vector
-!-- is necessary, but the data exchange between left and right PE's using
-!-- mpi-vectors is 10% faster than without.
+
+!
+!-- Define an MPI-datatype for the exchange of left/right boundaries.
+!-- Although data are contiguous in physical memory (which does not
+!-- necessarily require an MPI-derived datatype), the data exchange between
+!-- left and right PE's using the MPI-derived type is 10% faster than without.
     CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz(0), &
-                             MPI_REAL, type_xz(0), ierr )
+                          MPI_REAL, type_xz(0), ierr )
     CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
 
-    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), ierr) 
+    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
+                          ierr ) 
     CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
-!
-!-- Definition of mpi-vectors for multigrid
+
+!
+!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
     IF ( psolver == 'multigrid' )  THEN
 !    
-!--   The definition of mpi-vectors as aforementioned, but only 1 ghost point is used.
-       DO i = maximum_grid_level, 1 , -1
+!--    Definition of MPI-datatyoe as above, but only 1 ghost level is used
+       DO  i = maximum_grid_level, 1 , -1
+
           ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
 
           CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
-                             MPI_REAL, type_xz(i), ierr )
+                                MPI_REAL, type_xz(i), ierr )
           CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
 
-          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), ierr)
+          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
+                                ierr )
           CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
 
@@ -1069,6 +1053,8 @@
           nyn_l = nyn_l / 2
           nzt_l = nzt_l / 2
+
        ENDDO
-    END IF
+
+    ENDIF
 #endif
 
Index: palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 708)
+++ palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 709)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
@@ -10,6 +11,6 @@
 !
 ! 667 2010-12-23 12:06:00Z suehring/gryschka
-! Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
-! therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
+! Changed surface boundary conditions for u and v from mirror to Dirichlet.
+! Therefore u(uzb,:,:) and v(nzb,:,:) are now representative for height z0.
 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 708)
+++ palm/trunk/SOURCE/pres.f90	(revision 709)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! formatting adjustments
 !
 ! Former revisions:
@@ -104,5 +104,5 @@
 
     ddt_3d = 1.0 / dt_3d
-    d_weight_pres = 1. / weight_pres(intermediate_timestep_count)
+    d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)
 
 !
@@ -146,6 +146,5 @@
 !
 !-- Left/right
-
-    IF ( conserve_volume_flow  .AND.  ( outflow_l  .OR. outflow_r ) )  THEN
+    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
 
        volume_flow(1)   = 0.0
@@ -161,5 +160,5 @@
 !
 !--       Sum up the volume flow through the south/north boundary
-          DO  k = nzb_2d(j,i) + 1, nzt
+          DO  k = nzb_2d(j,i)+1, nzt
              volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
           ENDDO
@@ -173,9 +172,9 @@
        volume_flow = volume_flow_l  
 #endif
-       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
+       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
                                / volume_flow_area(1)
 
        DO  j = nysg, nyng
-          DO  k = nzb_2d(j,i) + 1, nzt
+          DO  k = nzb_2d(j,i)+1, nzt
              u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
           ENDDO
@@ -186,5 +185,5 @@
 !
 !-- South/north
-    IF ( conserve_volume_flow  .AND.  ( outflow_n  .OR. outflow_s ) )  THEN
+    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
 
        volume_flow(2)   = 0.0
@@ -200,5 +199,5 @@
 !
 !--       Sum up the volume flow through the south/north boundary
-          DO  k = nzb_2d(j,i) + 1, nzt
+          DO  k = nzb_2d(j,i)+1, nzt
              volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
           ENDDO
@@ -216,5 +215,5 @@
 
        DO  i = nxlg, nxrg
-          DO  k = nzb_v_inner(j,i) + 1, nzt
+          DO  k = nzb_v_inner(j,i)+1, nzt
              v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
           ENDDO
@@ -226,5 +225,5 @@
 !-- Remove mean vertical velocity
     IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
-       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner is not yet known
+       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner not yet known
           w_l = 0.0;  w_l_l = 0.0
           DO  i = nxl, nxr
@@ -237,6 +236,6 @@
 #if defined( __parallel )   
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
-                              ierr )
+          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
+                              comm2d, ierr )
 #else
           w_l = w_l_l  
@@ -574,5 +573,5 @@
 !-- Correction of the provisional velocities with the current perturbation 
 !-- pressure just computed
-    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .OR.  bc_ns_cyc ) )  THEN
+    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
        volume_flow_l(1) = 0.0
        volume_flow_l(2) = 0.0
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 708)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 709)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
@@ -152,6 +153,6 @@
     IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
     IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
-    IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
-       intermediate_timestep_count == 1 )  CALL ws_statistics
+    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
+         intermediate_timestep_count == 1 )  CALL ws_statistics
 
 !
@@ -927,6 +928,6 @@
     IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
     IF ( .NOT. constant_diffusion )  CALL production_e_init
-    IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
-       intermediate_timestep_count == 1 )  CALL ws_statistics
+    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
+         intermediate_timestep_count == 1 )  CALL ws_statistics
 
 
@@ -1459,6 +1460,6 @@
     IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
     IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
-    IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
-       intermediate_timestep_count == 1 )  CALL ws_statistics
+    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
+         intermediate_timestep_count == 1 )  CALL ws_statistics
 
 !
Index: palm/trunk/SOURCE/surface_coupler.f90
===================================================================
--- palm/trunk/SOURCE/surface_coupler.f90	(revision 708)
+++ palm/trunk/SOURCE/surface_coupler.f90	(revision 709)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! formatting adjustments
 !
 ! Former revisions:
@@ -10,6 +11,6 @@
 !
 ! 667 2010-12-23 12:06:00Z suehring/gryschka
-! additional case for nonequivalent processor and grid topopolgy in ocean and 
-! atmosphere added (coupling_topology = 1)
+! Additional case for nonequivalent processor and grid topopolgy in ocean and 
+! atmosphere added (coupling_topology = 1).
 ! Added exchange of u and v from Ocean to Atmosphere
 !
@@ -60,7 +61,7 @@
 
     IF ( coupling_topology == 0 ) THEN
-       CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, target_id,  &
-                          0,                                                    &
-                          terminate_coupled_remote, 1, MPI_INTEGER, target_id,  &
+       CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, target_id, &
+                          0,                                                   &
+                          terminate_coupled_remote, 1, MPI_INTEGER, target_id, &
                           0, comm_inter, status, ierr )
     ELSE
@@ -72,5 +73,6 @@
                              comm_inter, status, ierr )
        ENDIF
-       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
+       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
+                       ierr )
 
        ALLOCATE( total_2d_a(-nbgp:ny_a+nbgp,-nbgp:nx_a+nbgp),       &
@@ -97,25 +99,25 @@
 !-- Exchange the current simulated time between the models,
 !-- currently just for total_2ding
-    IF ( coupling_topology == 0 ) THEN   
-       CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, &
-                      target_id, 11, comm_inter, ierr )
-       CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
-                      target_id, 11, comm_inter, status, ierr )
+    IF ( coupling_topology == 0 ) THEN
+   
+       CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, 11, &
+                      comm_inter, ierr )
+       CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, target_id, &
+                      11, comm_inter, status, ierr )
     ELSE
+
        IF ( myid == 0 ) THEN
-          CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, &
-                         target_id, 11, comm_inter, ierr )
-          CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
+
+          CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
+                         11, comm_inter, ierr )
+          CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL,        &
                          target_id, 11, comm_inter, status, ierr )
+
        ENDIF
-       CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, &
-                       0, comm2d, ierr )
-    ENDIF
-    WRITE ( 9, * ) 'simulated time: ', simulated_time
-    WRITE ( 9, * ) 'time since start of coupling: ', &
-                  time_since_reference_point, ' remote: ', &
-                  time_since_reference_point_rem
-   CALL local_flush( 9 )
-  
+
+       CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
+                       ierr )
+
+    ENDIF
 
 !
@@ -124,47 +126,38 @@
     
 !
-!--    Horizontal grid size and number of processors is equal
-!--    in ocean and atmosphere
-       IF ( coupling_topology == 0 ) THEN
-
-!
-!--       Send heat flux at bottom surface to the ocean model
-          CALL MPI_SEND( shf(nysg,nxlg), ngp_xy, MPI_REAL, &
-                         target_id, 12, comm_inter, ierr )
-
-!
-!--       Send humidity flux at bottom surface to the ocean model
+!--    Horizontal grid size and number of processors is equal in ocean and
+!--    atmosphere
+       IF ( coupling_topology == 0 )  THEN
+
+!
+!--       Send heat flux at bottom surface to the ocean
+          CALL MPI_SEND( shf(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
+                         comm_inter, ierr )
+!
+!--       Send humidity flux at bottom surface to the ocean
           IF ( humidity )  THEN
-             CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, &
-                            target_id, 13, comm_inter, ierr )
-          ENDIF
-
-!
-!--       Receive temperature at the bottom surface from the ocean model
-          WRITE ( 9, * )  '*** receive pt from ocean'
-          CALL local_flush( 9 )
-          CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, &
-                         target_id, 14, comm_inter, status, ierr )
-
-!
-!--       Send the momentum flux (u) at bottom surface to the ocean model
-          CALL MPI_SEND( usws(nysg,nxlg), ngp_xy, MPI_REAL, &
-                         target_id, 15, comm_inter, ierr )
-
-!
-!--       Send the momentum flux (v) at bottom surface to the ocean model
-          CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, &
-                         target_id, 16, comm_inter, ierr )
-
-!
-!--       Receive u at the bottom surface from the ocean model
-          CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, &
-                         target_id, 17, comm_inter, status, ierr )
-
-!
-!--       Receive v at the bottom surface from the ocean model
-          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, &
-                         target_id, 18,  comm_inter, status, ierr )
-
+             CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 13, &
+                            comm_inter, ierr )
+          ENDIF
+!
+!--       Receive temperature at the bottom surface from the ocean
+          CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, target_id, 14, &
+                         comm_inter, status, ierr )
+!
+!--       Send the momentum flux (u) at bottom surface to the ocean
+          CALL MPI_SEND( usws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
+                         comm_inter, ierr )
+!
+!--       Send the momentum flux (v) at bottom surface to the ocean
+          CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
+                         comm_inter, ierr )
+!
+!--       Receive u at the bottom surface from the ocean
+          CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, target_id, 17, &
+                         comm_inter, status, ierr )
+!
+!--       Receive v at the bottom surface from the ocean
+          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18, &
+                         comm_inter, status, ierr )
 !
 !--    Horizontal grid size or number of processors differs between 
@@ -173,72 +166,68 @@
       
 !
-!--       Send heat flux at bottom surface to the ocean model
+!--       Send heat flux at bottom surface to the ocean
           total_2d_a = 0.0
-          total_2d = 0.0
+          total_2d   = 0.0
           total_2d(nys:nyn,nxl:nxr) = shf(nys:nyn,nxl:nxr)
-          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
-                           MPI_SUM, 0, comm2d, ierr )
-          CALL interpolate_to_ocean(12)
-   
-!
-!--       Send humidity flux at bottom surface to the ocean model
-          IF ( humidity ) THEN
+
+          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
+                           comm2d, ierr )
+          CALL interpolate_to_ocean( 12 )   
+!
+!--       Send humidity flux at bottom surface to the ocean
+          IF ( humidity )  THEN
              total_2d_a = 0.0
-             total_2d = 0.0
+             total_2d   = 0.0
              total_2d(nys:nyn,nxl:nxr) = qsws(nys:nyn,nxl:nxr)
-             CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
-                              MPI_SUM, 0, comm2d, ierr )
-             CALL interpolate_to_ocean(13)
-          ENDIF
-
-!
-!--       Receive temperature at the bottom surface from the ocean model
-          IF ( myid == 0 ) THEN
+
+             CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, &
+                              0, comm2d, ierr )
+             CALL interpolate_to_ocean( 13 )
+          ENDIF
+!
+!--       Receive temperature at the bottom surface from the ocean
+          IF ( myid == 0 )  THEN
              CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
                             target_id, 14, comm_inter, status, ierr )   
           ENDIF
           CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
-                          0, comm2d, ierr )
+          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
+                          ierr )
           pt(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
-
-!
-!--       Send momentum flux (u) at bottom surface to the ocean model
+!
+!--       Send momentum flux (u) at bottom surface to the ocean
           total_2d_a = 0.0 
-          total_2d = 0.0
+          total_2d   = 0.0
           total_2d(nys:nyn,nxl:nxr) = usws(nys:nyn,nxl:nxr)
-          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
-                           MPI_SUM, 0, comm2d, ierr )
-          CALL interpolate_to_ocean(15)
-
-!
-!--       Send momentum flux (v) at bottom surface to the ocean model
+          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
+                           comm2d, ierr )
+          CALL interpolate_to_ocean( 15 )
+!
+!--       Send momentum flux (v) at bottom surface to the ocean
           total_2d_a = 0.0
-          total_2d = 0.0
+          total_2d   = 0.0
           total_2d(nys:nyn,nxl:nxr) = vsws(nys:nyn,nxl:nxr)
-          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
-                           MPI_SUM, 0, comm2d, ierr )
-          CALL interpolate_to_ocean(16)
-
-!
-!--       Receive u at the bottom surface from the ocean model
-          IF ( myid == 0 ) THEN
+          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
+                           comm2d, ierr )
+          CALL interpolate_to_ocean( 16 )
+!
+!--       Receive u at the bottom surface from the ocean
+          IF ( myid == 0 )  THEN
              CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
-                            target_id, 17, comm_inter, status, ierr )            
+                            target_id, 17, comm_inter, status, ierr )
           ENDIF
           CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
-                          0, comm2d, ierr )
+          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
+                          ierr )
           u(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
-   
-!
-!--       Receive v at the bottom surface from the ocean model
-          IF ( myid == 0 ) THEN
+!
+!--       Receive v at the bottom surface from the ocean
+          IF ( myid == 0 )  THEN
              CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
-                            target_id, 18, comm_inter, status, ierr )            
+                            target_id, 18, comm_inter, status, ierr )
           ENDIF
           CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
-                          0, comm2d, ierr )
+          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
+                          ierr )
           v(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
 
@@ -252,128 +241,105 @@
        IF ( coupling_topology == 0 ) THEN
 !
-!--       Receive heat flux at the sea surface (top) from the atmosphere model
-          CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, &
-                         target_id, 12, comm_inter, status, ierr )
-
-
-!
-!--       Receive humidity flux from the atmosphere model (bottom)
+!--       Receive heat flux at the sea surface (top) from the atmosphere
+          CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
+                         comm_inter, status, ierr )
+!
+!--       Receive humidity flux from the atmosphere (bottom)
 !--       and add it to the heat flux at the sea surface (top)...
           IF ( humidity_remote )  THEN
              CALL MPI_RECV( qswst_remote(nysg,nxlg), ngp_xy, MPI_REAL, &
                             target_id, 13, comm_inter, status, ierr )
-
-          ENDIF
-
+          ENDIF
 !
 !--       Send sea surface temperature to the atmosphere model
-          CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, &
-                         target_id, 14, comm_inter, ierr )
-
+          CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, target_id, 14, &
+                         comm_inter, ierr )
 !
 !--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
-!--       model
-          WRITE ( 9, * )  '*** receive uswst from atmosphere'
-          CALL local_flush( 9 )
-          CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, &
-                         target_id, 15, comm_inter, status, ierr )
-
+          CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
+                         comm_inter, status, ierr )
 !
 !--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
-!--       model
-          CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, &
-                         target_id, 16, comm_inter, status, ierr )
-
-!--       Send u to the atmosphere model
-          CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, &
-                         target_id, 17, comm_inter, ierr )
-
-!
-!--       Send v to the atmosphere model
-          CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, &
-                         target_id, 18, comm_inter, ierr )
-
+          CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
+                         comm_inter, status, ierr )
+!
+!--       Send u to the atmosphere
+          CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, target_id, 17, &
+                         comm_inter, ierr )
+!
+!--       Send v to the atmosphere
+          CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, target_id, 18, &
+                         comm_inter, ierr )
 !
 !--    Horizontal gridsize or number of processors differs between 
 !--    ocean and atmosphere
        ELSE
-
-!
-!--       Receive heat flux at the sea surface (top) from the atmosphere model
-          IF ( myid == 0 ) THEN
+!
+!--       Receive heat flux at the sea surface (top) from the atmosphere
+          IF ( myid == 0 )  THEN
              CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                            target_id, 12, comm_inter, status, ierr )            
+                            target_id, 12, comm_inter, status, ierr )
+          ENDIF
+          CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
+                          ierr )
+          tswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
+!
+!--       Receive humidity flux at the sea surface (top) from the atmosphere
+          IF ( humidity_remote )  THEN
+             IF ( myid == 0 )  THEN
+                CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
+                               target_id, 13, comm_inter, status, ierr )
+             ENDIF
+             CALL MPI_BARRIER( comm2d, ierr )
+             CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, &
+                             comm2d, ierr)
+             qswst_remote(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
+          ENDIF
+!
+!--       Send surface temperature to atmosphere
+          total_2d_o = 0.0
+          total_2d   = 0.0
+          total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
+
+          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
+                           comm2d, ierr) 
+          CALL interpolate_to_atmos( 14 )
+!
+!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
+          IF ( myid == 0 )  THEN
+             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
+                            target_id, 15, comm_inter, status, ierr )
           ENDIF
           CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                          0, comm2d, ierr)
-          tswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
-
-!
-!--       Receive humidity flux at the sea surface (top) from the 
-!--       atmosphere model
-          IF ( humidity_remote ) THEN
-             IF ( myid == 0 ) THEN
-                CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                               target_id, 13, comm_inter, status, ierr )            
-             ENDIF
-             CALL MPI_BARRIER( comm2d, ierr )
-             CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                             0, comm2d, ierr)
-             qswst_remote(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
-          ENDIF
-
-!
-!--       Send surface temperature to atmosphere
-          total_2d_o = 0.0
-          total_2d = 0.0
-          total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
-
-          CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, &
-                          MPI_REAL, MPI_SUM, 0, comm2d, ierr) 
-
-          CALL interpolate_to_atmos(14)
-
-!
-!--       Receive momentum flux (u) at the sea surface (top) from the 
-!--       atmosphere model
-          IF ( myid == 0 ) THEN
+                          0, comm2d, ierr )
+          uswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
+!
+!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
+          IF ( myid == 0 )  THEN
              CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                            target_id, 15, comm_inter, status, ierr )            
+                            target_id, 16, comm_inter, status, ierr )
           ENDIF
           CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                          0, comm2d, ierr)
-          uswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
-
-!
-!--       Receive momentum flux (v) at the sea surface (top) from the 
-!--       atmosphere model
-          IF ( myid == 0 ) THEN
-             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                            target_id, 16, comm_inter, status, ierr )            
-          ENDIF
-          CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
-                          0, comm2d, ierr)
+          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
+                          ierr )
           vswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
-
 !
 !--       Send u to atmosphere
           total_2d_o = 0.0  
-          total_2d = 0.0
+          total_2d   = 0.0
           total_2d(nys:nyn,nxl:nxr) = u(nzt,nys:nyn,nxl:nxr)
-          CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, MPI_REAL, &
-                          MPI_SUM, 0, comm2d, ierr)  
-          CALL interpolate_to_atmos(17)
-
+          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
+                           comm2d, ierr )
+          CALL interpolate_to_atmos( 17 )
 !
 !--       Send v to atmosphere
           total_2d_o = 0.0
-          total_2d = 0.0
+          total_2d   = 0.0
           total_2d(nys:nyn,nxl:nxr) = v(nzt,nys:nyn,nxl:nxr)
-          CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, MPI_REAL, &
-                          MPI_SUM, 0, comm2d, ierr)  
-          CALL interpolate_to_atmos(18)
+          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
+                           comm2d, ierr )
+          CALL interpolate_to_atmos( 18 )
         
        ENDIF
@@ -382,8 +348,9 @@
 !--    Conversions of fluxes received from atmosphere
        IF ( humidity_remote )  THEN
-          !here tswst is still the sum of atmospheric bottom heat fluxes
-          tswst = tswst + qswst_remote * 2.2626108e6 / 1005.0
-          !*latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
-          !/(rho_atm(=1.0)*c_p)
+!
+!--       Here tswst is still the sum of atmospheric bottom heat fluxes,
+!--       * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
+!--       /(rho_atm(=1.0)*c_p)
+          tswst = tswst + qswst_remote * 2.2626108E6 / 1005.0
 !
 !--        ...and convert it to a salinity flux at the sea surface (top)
@@ -405,8 +372,7 @@
        vswst = vswst / rho(nzt,:,:)
 
-
-    ENDIF
-
-    IF ( coupling_topology == 1 ) THEN
+    ENDIF
+
+    IF ( coupling_topology == 1 )  THEN
        DEALLOCATE( total_2d_o, total_2d_a )
     ENDIF
@@ -420,5 +386,5 @@
 
 
-  SUBROUTINE interpolate_to_atmos(tag)
+  SUBROUTINE interpolate_to_atmos( tag )
 
     USE arrays_3d
@@ -430,5 +396,4 @@
     IMPLICIT NONE
 
-  
     INTEGER             ::  dnx, dnx2, dny, dny2, i, ii, j, jj
     INTEGER, intent(in) ::  tag
@@ -436,8 +401,7 @@
     CALL MPI_BARRIER( comm2d, ierr )
 
-    IF ( myid == 0 ) THEN
-
-! 
-!--    cyclic boundary conditions for the total 2D-grid
+    IF ( myid == 0 )  THEN
+!
+!--    Cyclic boundary conditions for the total 2D-grid
        total_2d_o(-nbgp:-1,:) = total_2d_o(ny+1-nbgp:ny,:)
        total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx+1-nbgp:nx)
@@ -452,6 +416,6 @@
 
 !
-!--    Distance for interpolation around coarse grid points within the fine grid 
-!--    (note: 2*dnx2 must not be equal with dnx)  
+!--    Distance for interpolation around coarse grid points within the fine
+!--    grid (note: 2*dnx2 must not be equal with dnx)
        dnx2 = 2 * ( dnx / 2 )
        dny2 = 2 * ( dny / 2 )
@@ -472,5 +436,5 @@
        ENDDO
 ! 
-!--    cyclic boundary conditions for atmosphere grid
+!--    Cyclic boundary conditions for atmosphere grid
        total_2d_a(-nbgp:-1,:) = total_2d_a(ny_a+1-nbgp:ny_a,:)
        total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx_a+1-nbgp:nx_a)
@@ -480,6 +444,6 @@
 !
 !--    Transfer of the atmosphere-grid-layer to the atmosphere
-       CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
-                      target_id, tag, comm_inter, ierr )
+       CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, target_id, &
+                      tag, comm_inter, ierr )
 
     ENDIF
@@ -490,5 +454,5 @@
 
 
-  SUBROUTINE interpolate_to_ocean(tag)
+  SUBROUTINE interpolate_to_ocean( tag )
 
     USE arrays_3d
@@ -500,19 +464,20 @@
     IMPLICIT NONE
 
-    REAL                ::  fl, fr, myl, myr
     INTEGER             ::  dnx, dny, i, ii, j, jj
     INTEGER, intent(in) ::  tag
+    REAL                ::  fl, fr, myl, myr
+
 
     CALL MPI_BARRIER( comm2d, ierr )
 
-    IF ( myid == 0 ) THEN   
-
-!
-!      Number of gridpoints of the fine grid within one mesh of the coarse grid
+    IF ( myid == 0 )  THEN   
+
+!
+!--    Number of gridpoints of the fine grid within one mesh of the coarse grid
        dnx = ( nx_o + 1 ) / ( nx_a + 1 ) 
        dny = ( ny_o + 1 ) / ( ny_a + 1 ) 
 
 ! 
-!--    cyclic boundary conditions for atmosphere grid
+!--    Cyclic boundary conditions for atmosphere grid
        total_2d_a(-nbgp:-1,:) = total_2d_a(ny+1-nbgp:ny,:)
        total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx+1-nbgp:nx)
@@ -521,5 +486,5 @@
        total_2d_a(:,nx+1:nx+nbgp) = total_2d_a(:,0:nbgp-1)
 !
-!--    Bilinear Interpolation from atmosphere-grid-layer to ocean-grid-layer
+!--    Bilinear Interpolation from atmosphere grid-layer to ocean grid-layer
        DO  j = 0, ny
           DO  i = 0, nx
@@ -527,6 +492,6 @@
              myr = ( total_2d_a(j+1,i+1) - total_2d_a(j,i+1) ) / dny
              DO  jj = 0, dny-1
-                fl = myl*jj  + total_2d_a(j,i)  
-                fr = myr*jj  + total_2d_a(j,i+1)  
+                fl = myl*jj + total_2d_a(j,i)  
+                fr = myr*jj + total_2d_a(j,i+1)  
                 DO  ii = 0, dnx-1
                    total_2d_o(j*dny+jj,i*dnx+ii) = ( fr - fl ) / dnx * ii + fl
@@ -536,5 +501,5 @@
        ENDDO
 ! 
-!--    cyclic boundary conditions for ocean grid
+!--    Cyclic boundary conditions for ocean grid
        total_2d_o(-nbgp:-1,:) = total_2d_o(ny_o+1-nbgp:ny_o,:)
        total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx_o+1-nbgp:nx_o)
@@ -542,5 +507,4 @@
        total_2d_o(ny_o+1:ny_o+nbgp,:) = total_2d_o(0:nbgp-1,:)
        total_2d_o(:,nx_o+1:nx_o+nbgp) = total_2d_o(:,0:nbgp-1)
-       
 
        CALL MPI_SEND( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
