Index: palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_s.f90	(revision 1310)
+++ palm/trunk/SOURCE/diffusion_s.f90	(revision 1320)
@@ -20,5 +20,11 @@
 ! Current revisions:
 ! ------------------
-!
+! ONLY-attribute added to USE-statements,
+! kind-parameters added to all INTEGER and REAL declaration statements,
+! kinds are defined in new module kinds,
+! old module precision_kind is removed,
+! revision history before 2012 removed,
+! comment fields (!:) to be used for variable explanations added to
+! all variable declaration statements
 !
 ! Former revisions:
@@ -47,26 +53,4 @@
 ! 1001 2012-09-13 14:08:46Z raasch
 ! some arrays comunicated by module instead of parameter list
-!
-! 667 2010-12-23 12:06:00Z suehring/gryschka
-! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
-!
-! 183 2008-08-04 15:39:12Z letzel
-! bugfix: calculation of fluxes at vertical surfaces
-!
-! 129 2007-10-30 12:12:24Z letzel
-! replace wall_heatflux by wall_s_flux that is now included in the parameter 
-! list, bugfix for assignment of fluxes at walls
-!
-! 20 2007-02-26 00:12:32Z raasch
-! Bugfix: ddzw dimensioned 1:nzt"+1"
-! Calculation extended for gridpoint nzt, fluxes can be given at top,
-! +s_flux_t in parameter list, s_flux renamed s_flux_b
-!
-! RCS Log replace by Id keyword, revision history cleaned up
-!
-! Revision 1.8  2006/02/23 10:34:17  raasch
-! nzb_2d replaced by nzb_s_outer in horizontal diffusion and by nzb_s_inner
-! or nzb_diff_s_inner, respectively, in vertical diffusion, prescribed surface
-! fluxes at vertically oriented topography
 !
 ! Revision 1.1  2000/04/13 14:54:02  schroeter
@@ -99,18 +83,30 @@
     SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
 
-       USE arrays_3d
-       USE control_parameters
-       USE grid_variables
-       USE indices
+       USE arrays_3d,                                                          &
+           ONLY:  ddzu, ddzw, kh, tend
+       
+       USE control_parameters,                                                 & 
+           ONLY: use_surface_fluxes, use_top_fluxes
+       
+       USE grid_variables,                                                     &
+           ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+       
+       USE indices,                                                            &
+           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg,                  &
+                  nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, nzt, nzt_diff
+       
+       USE kinds
 
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
-       REAL    ::  wall_s_flux(0:4)
-       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
+       INTEGER(iwp) ::  i                 !: 
+       INTEGER(iwp) ::  j                 !: 
+       INTEGER(iwp) ::  k                 !: 
+       REAL(wp)     ::  wall_s_flux(0:4)  !: 
+       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t !:
 #if defined( __nopointer )
-       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !: 
 #else
-       REAL, DIMENSION(:,:,:), POINTER ::  s
+       REAL(wp), DIMENSION(:,:,:), POINTER ::  s  !: 
 #endif
 
@@ -121,12 +117,12 @@
              DO  k = nzb_s_outer(j,i)+1, nzt
 
-                tend(k,j,i) = tend(k,j,i)                                     &
-                                          + 0.5 * (                           &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
-                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
-                                                  ) * ddx2                    &
-                                          + 0.5 * (                           &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
-                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
+                tend(k,j,i) = tend(k,j,i)                                      &
+                                          + 0.5 * (                            &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
+                                                  ) * ddx2                     &
+                                          + 0.5 * (                            &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
                                                   ) * ddy2
              ENDDO
@@ -138,18 +134,18 @@
                 DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
 
-                   tend(k,j,i) = tend(k,j,i)                                  &
-                                                + ( fwxp(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
-                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
-                                                   -fwxm(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
-                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
-                                                  ) * ddx2                    &
-                                                + ( fwyp(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
-                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
-                                                   -fwym(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
-                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                                + ( fwxp(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
+                                                   -fwxm(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
+                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
+                                                  ) * ddx2                     &
+                                                + ( fwyp(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
+                                                   -fwym(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
+                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
                                                   ) * ddy2
                 ENDDO
@@ -162,8 +158,8 @@
              DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
-                tend(k,j,i) = tend(k,j,i)                                     &
-                                       + 0.5 * (                              &
-            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
-          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
+                tend(k,j,i) = tend(k,j,i)                                      &
+                                       + 0.5 * (                               &
+            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
                                                ) * ddzw(k)
              ENDDO
@@ -176,9 +172,9 @@
                 k = nzb_s_inner(j,i)+1
 
-                tend(k,j,i) = tend(k,j,i)                                     &
-                                       + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )    &
-                                               * ( s(k+1,j,i)-s(k,j,i) )      &
-                                               * ddzu(k+1)                    &
-                                           + s_flux_b(j,i)                    &
+                tend(k,j,i) = tend(k,j,i)                                      &
+                                       + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )     &
+                                               * ( s(k+1,j,i)-s(k,j,i) )       &
+                                               * ddzu(k+1)                     &
+                                           + s_flux_b(j,i)                     &
                                          ) * ddzw(k)
 
@@ -192,9 +188,9 @@
                 k = nzt
 
-                tend(k,j,i) = tend(k,j,i)                                     &
-                                       + ( - s_flux_t(j,i)                    &
-                                           - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
-                                                 * ( s(k,j,i)-s(k-1,j,i) )    &
-                                                 * ddzu(k)                    &
+                tend(k,j,i) = tend(k,j,i)                                      &
+                                       + ( - s_flux_t(j,i)                     &
+                                           - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )   &
+                                                 * ( s(k,j,i)-s(k-1,j,i) )     &
+                                                 * ddzu(k)                     &
                                          ) * ddzw(k)
 
@@ -212,18 +208,31 @@
     SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
 
-       USE arrays_3d
-       USE control_parameters
-       USE grid_variables
-       USE indices
+       USE arrays_3d,                                                          &
+           ONLY:  ddzu, ddzw, kh, tend
+           
+       USE control_parameters,                                                 & 
+           ONLY: use_surface_fluxes, use_top_fluxes
+       
+       USE grid_variables,                                                     &
+           ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+       
+       USE indices, &
+           ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,    &
+                 nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, nzt, nzt_diff
+           
+       USE kinds
 
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
-       REAL    ::  wall_s_flux(0:4)
-       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
+       INTEGER(iwp) ::  i                 !: 
+       INTEGER(iwp) ::  j                 !: 
+       INTEGER(iwp) ::  k                 !: 
+       REAL(wp)     ::  wall_s_flux(0:4)  !: 
+       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b !: 
+       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_t !: 
 #if defined( __nopointer )
-       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !:
 #else
-       REAL, DIMENSION(:,:,:), POINTER ::  s
+       REAL(wp), DIMENSION(:,:,:), POINTER ::  s  !:
 #endif
 
@@ -239,12 +248,12 @@
                 IF ( k > nzb_s_outer(j,i) )  THEN
 
-                   tend(k,j,i) = tend(k,j,i)                                  &
-                                          + 0.5 * (                           &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
-                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
-                                                  ) * ddx2                    &
-                                          + 0.5 * (                           &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
-                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                          + 0.5 * (                            &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
+                                                  ) * ddx2                     &
+                                          + 0.5 * (                            &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
                                                   ) * ddy2
                 ENDIF
@@ -257,18 +266,18 @@
                      ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 ) )    &
                 THEN
-                   tend(k,j,i) = tend(k,j,i)                                  &
-                                                + ( fwxp(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
-                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
-                                                   -fwxm(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
-                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
-                                                  ) * ddx2                    &
-                                                + ( fwyp(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
-                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
-                                                   -fwym(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
-                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                                + ( fwxp(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
+                                                   -fwxm(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
+                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
+                                                  ) * ddx2                     &
+                                                + ( fwyp(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
+                                                   -fwym(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
+                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
                                                   ) * ddy2
                 ENDIF
@@ -281,8 +290,8 @@
              DO  k = 1, nzt_diff
                 IF ( k >= nzb_diff_s_inner(j,i) )  THEN
-                   tend(k,j,i) = tend(k,j,i)                                  &
-                                       + 0.5 * (                              &
-            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
-          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                       + 0.5 * (                               &
+            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
                                                ) * ddzw(k)
                 ENDIF
@@ -294,9 +303,9 @@
              DO  k = 1, nzt
                 IF ( use_surface_fluxes  .AND.  k == nzb_s_inner(j,i)+1 )  THEN
-                   tend(k,j,i) = tend(k,j,i)                                  &
-                                          + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
-                                                  * ( s(k+1,j,i)-s(k,j,i) )   &
-                                                  * ddzu(k+1)                 &
-                                              + s_flux_b(j,i)                 &
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                          + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
+                                                  * ( s(k+1,j,i)-s(k,j,i) )    &
+                                                  * ddzu(k+1)                  &
+                                              + s_flux_b(j,i)                  &
                                             ) * ddzw(k)
                 ENDIF
@@ -327,18 +336,31 @@
     SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
 
-       USE arrays_3d
-       USE control_parameters
-       USE grid_variables
-       USE indices
+       USE arrays_3d,                                                          &
+           ONLY:  ddzu, ddzw, kh, tend
+           
+       USE control_parameters,                                                 & 
+           ONLY: use_surface_fluxes, use_top_fluxes
+       
+       USE grid_variables,                                                     &
+           ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
+       
+       USE indices,                                                            &
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb_diff_s_inner, nzb_s_inner,       &
+                  nzb_s_outer, nzt, nzt_diff
+       
+       USE kinds
 
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
-       REAL    ::  wall_s_flux(0:4)
-       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
+       INTEGER(iwp) ::  i                 !: 
+       INTEGER(iwp) ::  j                 !: 
+       INTEGER(iwp) ::  k                 !: 
+       REAL(wp)     ::  wall_s_flux(0:4)  !: 
+       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b  !: 
+       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_t  !: 
 #if defined( __nopointer )
-       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s !: 
 #else
-       REAL, DIMENSION(:,:,:), POINTER ::  s
+       REAL(wp), DIMENSION(:,:,:), POINTER ::  s  !: 
 #endif
 
@@ -347,12 +369,12 @@
        DO  k = nzb_s_outer(j,i)+1, nzt
 
-          tend(k,j,i) = tend(k,j,i)                                           &
-                                          + 0.5 * (                           &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
-                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
-                                                  ) * ddx2                    &
-                                          + 0.5 * (                           &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
-                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
+          tend(k,j,i) = tend(k,j,i)                                            &
+                                          + 0.5 * (                            &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
+                                                  ) * ddx2                     &
+                                          + 0.5 * (                            &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
                                                   ) * ddy2
        ENDDO
@@ -360,22 +382,22 @@
 !
 !--    Apply prescribed horizontal wall heatflux where necessary
-       IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
+       IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) )       &
        THEN
           DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
 
-             tend(k,j,i) = tend(k,j,i)                                        &
-                                                + ( fwxp(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
-                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
-                                                   -fwxm(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
-                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
-                                                  ) * ddx2                    &
-                                                + ( fwyp(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
-                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
-                                                   -fwym(j,i) * 0.5 *         &
-                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
-                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
+             tend(k,j,i) = tend(k,j,i)                                         &
+                                                + ( fwxp(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
+                                                   -fwxm(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
+                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
+                                                  ) * ddx2                     &
+                                                + ( fwyp(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
+                                                   -fwym(j,i) * 0.5 *          &
+                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
+                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
                                                   ) * ddy2
           ENDDO
@@ -388,8 +410,8 @@
        DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
-          tend(k,j,i) = tend(k,j,i)                                           &
-                                       + 0.5 * (                              &
-            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
-          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
+          tend(k,j,i) = tend(k,j,i)                                            &
+                                       + 0.5 * (                               &
+            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
+          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
                                                ) * ddzw(k)
        ENDDO
@@ -401,8 +423,8 @@
           k = nzb_s_inner(j,i)+1
 
-          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
-                                            * ( s(k+1,j,i)-s(k,j,i) )    &
-                                            * ddzu(k+1)                  &
-                                        + s_flux_b(j,i)                  &
+          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )        &
+                                            * ( s(k+1,j,i)-s(k,j,i) )          &
+                                            * ddzu(k+1)                        &
+                                        + s_flux_b(j,i)                        &
                                       ) * ddzw(k)
 
@@ -415,8 +437,8 @@
           k = nzt
 
-          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                  &
-                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
-                                            * ( s(k,j,i)-s(k-1,j,i) )    &
-                                            * ddzu(k)                    &
+          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                        &
+                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )        &
+                                            * ( s(k,j,i)-s(k-1,j,i) )          &
+                                            * ddzu(k)                          &
                                       ) * ddzw(k)
 
