Index: /palm/trunk/DOC/app/chapter_5.1.html
===================================================================
--- /palm/trunk/DOC/app/chapter_5.1.html	(revision 93)
+++ /palm/trunk/DOC/app/chapter_5.1.html	(revision 94)
@@ -36,5 +36,5 @@
 ~/palm/release-3.1c</span><br style="font-family: Courier New,Courier,monospace;"><span style="font-family: Courier New,Courier,monospace;">svn
 checkout --username &lt;your username&gt;
-svn://130.75.105.45:palm/tags/release-3.1c trunk</span></p>
+svn://130.75.105.45/palm/tags/release-3.1c trunk</span></p>
 <p>However, this would require to carry out again the complete
 installation process described in chapter 5.0. So far, differet
@@ -53,5 +53,8 @@
 -u<br>mbuild</blockquote>
 <p>on all remote hosts listed in the configuration file<span style="font-family: Courier New,Courier,monospace;">
-.mrun.config</span>.</p><p>You can use "subversion" for code comparison between the different versions. Also,&nbsp;modified code can be committed to the repository, but this is restricted to&nbsp;PALM developers. </p>
+.mrun.config</span>.</p><p>You
+can use "subversion" for code comparison between the different
+versions. Also,&nbsp;modified code can be committed to the repository,
+but this is restricted to&nbsp;PALM developers. </p>
 <p>As a last step, a suitable test run should be carried out. It
 should
Index: /palm/trunk/DOC/app/index.html
===================================================================
--- /palm/trunk/DOC/app/index.html	(revision 93)
+++ /palm/trunk/DOC/app/index.html	(revision 94)
@@ -148,5 +148,5 @@
 <p><span style="font-family: Thorndale;" lang="EN-GB">Current
 model
-version:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 3.2a<br>For date of last change see bottom line of each page. <o:p></o:p></span></p>
+version:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 3.2b<br>For date of last change see bottom line of each page. <o:p></o:p></span></p>
 <div style="text-align: center;" class="MsoNormal" align="center"><span style="font-family: Thorndale;">
 <hr align="center" size="2" width="100%"></span></div>
Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 93)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 94)
@@ -1,4 +1,8 @@
 New:
 ---
+ocean version including prognostic equation for salinity and equation of state for seawater
++ inipar-parameters ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinity_flux
+
+check_parameters, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, parin, prognostic_equations, read_var_list, write_var_list, write_3d_binary
 
 
@@ -7,4 +11,27 @@
 
 
+
+
+
 Errors:
 ------
+
+Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).
+
+diffusion_e, prognostic_equations
+
+
+
+
+To be completed (adjusted to ocean version):
+----------------
+
+boundary conditions for salinity (check_parameters)
+output of salinity profile (header)
+
+initialization of ug,vg profiles
+
+calculation of w*
+output of theta* and z_i (e.g. run_control file)
+
+default disturbance level
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 93)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 94)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Initial salinity profile is calculated.
+! z_max_do1d is checked only in case of ocean = .f.,
+! +initial temperature profile for the ocean version,
 !
 ! Former revisions:
@@ -153,4 +155,21 @@
        ENDIF
     ENDIF
+
+!
+!-- Check ocean setting
+    IF ( ocean )  THEN
+       action = ' '
+       IF ( timestep_scheme(1:8) == 'leapfrog' )  THEN
+          WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
+       ENDIF
+       IF ( action /= ' ' )  THEN
+          IF ( myid == 0 )  THEN
+             PRINT*, '+++ check_parameters:'
+             PRINT*, '    ocean = .T. does not allow ', TRIM( action )
+          ENDIF
+          CALL local_stop
+       ENDIF
+    ENDIF
+
 !
 !-- Check whether there are any illegal values
@@ -418,6 +437,7 @@
        v_init  = vg_surface
        pt_init = pt_surface
-       IF ( humidity )        q_init = q_surface
-       IF ( passive_scalar )  q_init = s_surface
+       IF ( humidity )        q_init  = q_surface
+       IF ( ocean )           sa_init = sa_surface
+       IF ( passive_scalar )  q_init  = s_surface
 
 !
@@ -504,43 +524,76 @@
 
 !
-!--    If required, compute initial temperature profile using the given
-!--    temperature gradient
+!--    Compute initial temperature profile using the given temperature gradients
        i = 1
        gradient = 0.0
-       pt_vertical_gradient_level_ind(1) = 0
-       DO  k = 1, nzt+1
-          IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
-               pt_vertical_gradient_level(i) >= 0.0 )  THEN
-             gradient = pt_vertical_gradient(i) / 100.0
-             pt_vertical_gradient_level_ind(i) = k - 1
-             i = i + 1
-             IF ( i > 10 )  THEN
-                IF ( myid == 0 )  THEN
-                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
-                           ' "pt_vertical_gradient_level_ind" exceeded'
-                ENDIF
-                CALL local_stop
-             ENDIF
-          ENDIF
-          IF ( gradient /= 0.0 )  THEN
-             IF ( k /= 1 )  THEN
-                pt_init(k) = pt_init(k-1) + dzu(k) * gradient
+
+       IF ( .NOT. ocean )  THEN
+
+          pt_vertical_gradient_level_ind(1) = 0
+          DO  k = 1, nzt+1
+             IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
+                  pt_vertical_gradient_level(i) >= 0.0 )  THEN
+                gradient = pt_vertical_gradient(i) / 100.0
+                pt_vertical_gradient_level_ind(i) = k - 1
+                i = i + 1
+                IF ( i > 10 )  THEN
+                   IF ( myid == 0 )  THEN
+                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
+                              ' "pt_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= 1 )  THEN
+                   pt_init(k) = pt_init(k-1) + dzu(k) * gradient
+                ELSE
+                   pt_init(k) = pt_surface   + 0.5 * dzu(k) * gradient
+                ENDIF
              ELSE
-                pt_init(k) = pt_init(k-1) + 0.5 * dzu(k) * gradient
-             ENDIF
-          ELSE
-             pt_init(k) = pt_init(k-1)
-          ENDIF
-       ENDDO
+                pt_init(k) = pt_init(k-1)
+             ENDIF
+          ENDDO
+
+       ELSE
+
+          pt_vertical_gradient_level_ind(1) = nzt+1
+          DO  k = nzt, 0, -1
+             IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
+                  pt_vertical_gradient_level(i) <= 0.0 )  THEN
+                gradient = pt_vertical_gradient(i) / 100.0
+                pt_vertical_gradient_level_ind(i) = k + 1
+                i = i + 1
+                IF ( i > 10 )  THEN
+                   IF ( myid == 0 )  THEN
+                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
+                              ' "pt_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= nzt )  THEN
+                   pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
+                ELSE
+                   pt_init(k)   = pt_surface - 0.5 * dzu(k+1) * gradient
+                   pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient
+                ENDIF
+             ELSE
+                pt_init(k) = pt_init(k+1)
+             ENDIF
+          ENDDO
+
+       ENDIF
 
 !
 !--    In case of no given temperature gradients, choose gradient of neutral
 !--    stratification
-       IF ( pt_vertical_gradient_level(1) == -1.0 )  THEN
+       IF ( pt_vertical_gradient_level(1) == -9999999.9 )  THEN
           pt_vertical_gradient_level(1) = 0.0
        ENDIF
 
 !
-!--    Store temperature gradient at the top boundary for possile Neumann
+!--    Store temperature gradient at the top boundary for possible Neumann
 !--    boundary condition
        bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
@@ -609,4 +662,41 @@
        ENDIF
 
+!
+!--    If required, compute initial salinity profile using the given salinity
+!--    gradients
+       IF ( ocean )  THEN
+
+          i = 1
+          gradient = 0.0
+
+          sa_vertical_gradient_level_ind(1) = nzt+1
+          DO  k = nzt, 0, -1
+             IF ( sa_vertical_gradient_level(i) > zu(k)  .AND. &
+                  sa_vertical_gradient_level(i) <= 0.0 )  THEN
+                gradient = sa_vertical_gradient(i) / 100.0
+                sa_vertical_gradient_level_ind(i) = k + 1
+                i = i + 1
+                IF ( i > 10 )  THEN
+                   IF ( myid == 0 )  THEN
+                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
+                              ' "sa_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= nzt )  THEN
+                   sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
+                ELSE
+                   sa_init(k)   = sa_surface - 0.5 * dzu(k+1) * gradient
+                   sa_init(k+1) = sa_surface + 0.5 * dzu(k+1) * gradient
+                ENDIF
+             ELSE
+                sa_init(k) = sa_init(k+1)
+             ENDIF
+          ENDDO
+
+       ENDIF
+
     ENDIF
 
@@ -1769,14 +1859,16 @@
 !-- y-value range of the coordinate system (PROFIL).
 !-- x-value range determined in plot_1d.
-    cross_uymin = 0.0
-    IF ( z_max_do1d == -1.0 )  THEN
-       cross_uymax = zu(nzt+1)
-    ELSEIF ( z_max_do1d < zu(nzb+1)  .OR.  z_max_do1d > zu(nzt+1) )  THEN
-       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do1d=',  &
-                                 z_max_do1d,' must be >= ', zu(nzb+1),  &
-                                 ' or <= ', zu(nzt+1)
-       CALL local_stop
-    ELSE
-       cross_uymax = z_max_do1d
+    IF ( .NOT. ocean )  THEN
+       cross_uymin = 0.0
+       IF ( z_max_do1d == -1.0 )  THEN
+          cross_uymax = zu(nzt+1)
+       ELSEIF ( z_max_do1d < zu(nzb+1)  .OR.  z_max_do1d > zu(nzt+1) )  THEN
+          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do1d=',  &
+                                    z_max_do1d,' must be >= ', zu(nzb+1),  &
+                                    ' or <= ', zu(nzt+1)
+          CALL local_stop
+       ELSE
+          cross_uymax = z_max_do1d
+       ENDIF
     ENDIF
 
Index: /palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_e.f90	(revision 93)
+++ /palm/trunk/SOURCE/diffusion_e.f90	(revision 94)
@@ -4,5 +4,8 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Adjustment of mixing length calculation for the ocean version. zw added to
+! argument list.
+! This is also a bugfix, because the height above the topography is now
+! used instead of the height above level k=0.
 !
 ! Former revisions:
@@ -48,5 +51,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
-                            rif, tend, zu )
+                            rif, tend, zu, zw )
 
        USE control_parameters
@@ -60,5 +63,5 @@
        REAL            ::  dpt_dz, l_stable, phi_m
        REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
-                           l_grid(1:nzt), zu(0:nzt+1)
+                           l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
        REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
        REAL, DIMENSION(:,:), POINTER   ::  rif
@@ -98,7 +101,10 @@
 !--                Adjustment of the mixing length
                    IF ( wall_adjustment )  THEN
-                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
-                                     l_stable )
-                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
+                      l(k,j)  = MIN( wall_adjustment_factor *          &
+                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                     l_grid(k), l_stable )
+                      ll(k,j) = MIN( wall_adjustment_factor *          &
+                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                     l_grid(k) )
                    ELSE
                       l(k,j)  = MIN( l_grid(k), l_stable )
@@ -106,6 +112,10 @@
                    ENDIF
                    IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
-                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
+                      l(k,j)  = MIN( l(k,j),  kappa *                          &
+                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
+                                              / phi_m )
+                      ll(k,j) = MIN( ll(k,j), kappa *                          &
+                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
+                                              / phi_m )
                    ENDIF
 
@@ -180,7 +190,10 @@
 !--                Adjustment of the mixing length
                    IF ( wall_adjustment )  THEN
-                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
-                                     l_stable )
-                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
+                      l(k,j)  = MIN( wall_adjustment_factor *          &
+                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                     l_grid(k), l_stable )
+                      ll(k,j) = MIN( wall_adjustment_factor *          &
+                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                     l_grid(k) )
                    ELSE
                       l(k,j)  = MIN( l_grid(k), l_stable )
@@ -188,6 +201,10 @@
                    ENDIF
                    IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
-                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
+                      l(k,j)  = MIN( l(k,j),  kappa *                          &
+                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
+                                              / phi_m )
+                      ll(k,j) = MIN( ll(k,j), kappa *                          &
+                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
+                                              / phi_m )
                    ENDIF
 
@@ -253,5 +270,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
-                               theta, rif, tend, zu )
+                               theta, rif, tend, zu, zw )
 
        USE control_parameters
@@ -265,5 +282,5 @@
        REAL            ::  dpt_dz, l_stable, phi_m
        REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
-                           l_grid(1:nzt), zu(0:nzt+1)
+                           l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
        REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
        REAL, DIMENSION(:,:), POINTER   ::  rif
@@ -301,6 +318,9 @@
 !--       Adjustment of the mixing length
           IF ( wall_adjustment )  THEN
-             l(k)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
-             ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
+             l(k)  = MIN( wall_adjustment_factor *                     &
+                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
+                          l_stable )
+             ll(k) = MIN( wall_adjustment_factor *                     &
+                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
           ELSE
              l(k)  = MIN( l_grid(k), l_stable )
@@ -308,6 +328,8 @@
           ENDIF
           IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-             l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
-             ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
+             l(k)  = MIN( l(k),  kappa * &
+                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
+             ll(k) = MIN( ll(k), kappa * &
+                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
           ENDIF
 
Index: /palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- /palm/trunk/SOURCE/diffusivities.f90	(revision 93)
+++ /palm/trunk/SOURCE/diffusivities.f90	(revision 94)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Adjustment of mixing length calculation for the ocean version.
+! This is also a bugfix, because the height above the topography is now
+! used instead of the height above level k=0.
 ! 
 ! Former revisions:
@@ -114,6 +116,8 @@
              ENDIF
              IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
-                ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
+                l(k)  = MIN( l(k),  kappa * &
+                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
+                ll(k) = MIN( ll(k), kappa * &
+                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
              ENDIF
 
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 93)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 94)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation of Deardorff velocity scale adjusted to be used with the ocean
+! version (HAS STILL TO BE COMPLETED!!!)
 !
 ! Former revisions:
@@ -784,26 +785,30 @@
 !
 !--    Determine the boundary layer height using two different schemes.
-!--    First scheme: Starting from the Earth's surface, look for the first 
-!--    relative minimum of the total heat flux. The corresponding height is
-!--    accepted as the boundary layer height, if it is less than 1.5 times the
-!--    height where the heat flux becomes negative for the first time.
+!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
+!--    first relative minimum (maximum) of the total heat flux.
+!--    The corresponding height is assumed as the boundary layer height, if it
+!--    is less than 1.5 times the height where the heat flux becomes negative
+!--    (positive) for the first time.
 !--    NOTE: This criterion is still capable of improving!
        z_i(1) = 0.0
        first = .TRUE.
-       DO  k = nzb, nzt-1
-          IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
-             first = .FALSE.
-             height = zw(k)
-          ENDIF
-          IF ( hom(k,1,18,sr) < 0.0  .AND. &
-               hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
-             IF ( zw(k) < 1.5 * height )  THEN
-                z_i(1) = zw(k)
-             ELSE
-                z_i(1) = height
+!       IF ( .NOT. ocean )  THEN
+          DO  k = nzb, nzt-1
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
+                first = .FALSE.
+                height = zw(k)
              ENDIF
-             EXIT
-          ENDIF
-       ENDDO
+             IF ( hom(k,1,18,sr) < 0.0  .AND. &
+                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
+                IF ( zw(k) < 1.5 * height )  THEN
+                   z_i(1) = zw(k)
+                ELSE
+                   z_i(1) = height
+                ENDIF
+                EXIT
+             ENDIF
+          ENDDO
+!       ELSE
+!       ENDIF
 
 !
@@ -829,5 +834,6 @@
        IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
           hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
-                                       hom(nzb,1,18,sr) * z_i(1) )**0.333333333
+                                       hom(nzb,1,18,sr) *      &
+                                       ABS( z_i(1) ) )**0.333333333
 !--       so far this only works if Prandtl layer is used
           hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 93)
+++ /palm/trunk/SOURCE/header.f90	(revision 94)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Output of model height, stretch level, Prandtl-layer height and initial
+! temperature profile adjusted to be used also with the ocean version. 
 !
 ! Former revisions:
@@ -287,8 +288,16 @@
 !
 !-- Computational grid
-    WRITE ( io, 250 )  dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(nzt+1)
-    IF ( dz_stretch_level_index < nzt+1 )  THEN
-       WRITE ( io, 252 )  dz_stretch_level, dz_stretch_level_index, &
-                          dz_stretch_factor, dz_max
+    IF ( .NOT. ocean )  THEN
+       WRITE ( io, 250 )  dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(nzt+1)
+       IF ( dz_stretch_level_index < nzt+1 )  THEN
+          WRITE ( io, 252 )  dz_stretch_level, dz_stretch_level_index, &
+                             dz_stretch_factor, dz_max
+       ENDIF
+    ELSE
+       WRITE ( io, 250 )  dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(0)
+       IF ( dz_stretch_level_index > 0 )  THEN
+          WRITE ( io, 252 )  dz_stretch_level, dz_stretch_level_index, &
+                             dz_stretch_factor, dz_max
+       ENDIF
     ENDIF
     WRITE ( io, 254 )  nx, ny, nzt+1, MIN( nnx, nx+1 ), MIN( nny, ny+1 ), &
@@ -440,5 +449,6 @@
 
     IF ( prandtl_layer )  THEN
-       WRITE ( io, 305 )  zu(1), roughness_length, kappa, rif_min, rif_max
+       WRITE ( io, 305 )  0.5 * (zu(1)-zu(0)), roughness_length, kappa, &
+                          rif_min, rif_max
        IF ( .NOT. constant_heatflux )  WRITE ( io, 308 )
        IF ( humidity  .AND.  .NOT. constant_waterflux )  THEN
@@ -950,15 +960,15 @@
     DO  WHILE ( pt_vertical_gradient_level_ind(i) /= -9999 )
 
-       WRITE (coor_chr,'(F6.2,4X)')  pt_init(pt_vertical_gradient_level_ind(i))
-       temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F6.2,4X)')  pt_vertical_gradient(i)
-       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(I6,4X)')  pt_vertical_gradient_level_ind(i)
-       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F6.1,4X)')  pt_vertical_gradient_level(i)
-       coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
+       WRITE (coor_chr,'(F7.2)')  pt_init(pt_vertical_gradient_level_ind(i))
+       temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F7.2)')  pt_vertical_gradient(i)
+       gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(I7)')  pt_vertical_gradient_level_ind(i)
+       slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F7.1)')  pt_vertical_gradient_level(i)
+       coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
 
        i = i + 1
@@ -1384,8 +1394,8 @@
 470 FORMAT (//' Actions during the simulation:'/ &
               ' -----------------------------'/)
-471 FORMAT ('    Disturbance impulse (u,v) every :  ',F6.2,' s'/             &
-            '    Disturbance amplitude           :    ',F4.2, ' m/s'/        &
-            '    Lower disturbance level         : ',F7.2,' m (GP ',I4,')'/  &
-            '    Upper disturbance level         : ',F7.2,' m (GP ',I4,')')
+471 FORMAT ('    Disturbance impulse (u,v) every :   ',F6.2,' s'/            &
+            '    Disturbance amplitude           :     ',F4.2, ' m/s'/       &
+            '    Lower disturbance level         : ',F8.2,' m (GP ',I4,')'/  &
+            '    Upper disturbance level         : ',F8.2,' m (GP ',I4,')')
 472 FORMAT ('    Disturbances continued during the run from i/j =',I4, &
                  ' to i/j =',I4)
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 93)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 94)
@@ -7,5 +7,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Initialization of salinity
 !
 ! Former revisions:
@@ -187,4 +187,12 @@
     ENDIF
 
+    IF ( ocean )  THEN
+       ALLOCATE( saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE( r(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
+                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
+                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
+                 sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+    ENDIF
+
 !
 !-- 3D-array for storing the dissipation, needed for calculating the sgs
@@ -285,5 +293,5 @@
           qsws   => qsws_1
           qswst  => qswst_1
-          q      => q_1;     q_p  => q_2;     tq_m => q_3;    q_m => q_3
+          q      => q_1;     q_p  => q_2;     tq_m  => q_3;    q_m => q_3
           IF ( humidity )        vpt  => vpt_1
           IF ( cloud_physics )   ql   => ql_1
@@ -292,4 +300,9 @@
              ql_c => ql_2
           ENDIF
+       ENDIF
+
+       IF ( ocean )  THEN
+          saswst => saswst_1
+          sa     => sa_1;    sa_p => sa_2;    tsa_m => sa_3
        ENDIF
 
@@ -457,4 +470,11 @@
           ENDIF
 
+          IF ( ocean )  THEN
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   sa(:,j,i) = sa_init
+                ENDDO
+             ENDDO
+          ENDIF
           
           IF ( constant_diffusion )  THEN
@@ -566,6 +586,6 @@
 !
 !--    Initialize fluxes at top surface
-!--    Currently, only the heatflux can be prescribed. The latent flux is
-!--    zero in this case!
+!--    Currently, only the heatflux and salinity flux can be prescribed.
+!--    The latent flux is zero in this case!
        IF ( use_top_fluxes )  THEN
 
@@ -579,4 +599,8 @@
                 qswst = 0.0
                 IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
+             ENDIF
+
+             IF ( ocean )  THEN
+                saswst = top_salinityflux
              ENDIF
          ENDIF
@@ -744,4 +768,9 @@
        ENDIF
 
+       IF ( ocean )  THEN
+          tsa_m = 0.0
+          sa_p  = sa
+       ENDIF
+
 !
 !--    Initialize old timelevels needed for radiation boundary conditions
@@ -783,4 +812,5 @@
        e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
        IF ( humidity  .OR.  passive_scalar )  q_p = q
+       IF ( ocean )  sa_p = sa
 
     ELSE
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 93)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 94)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Grid definition for ocean version
 !
 ! Former revisions:
@@ -70,31 +70,73 @@
        CALL local_stop
     ENDIF
-!
-!-- Since the w-level lies on the surface, the first u-level (staggered!) lies
-!-- below the surface (used for "mirror" boundary condition).
-!-- The first u-level above the surface corresponds to the top of the 
-!-- Prandtl-layer.
-    zu(0) = - dz * 0.5
-    zu(1) =   dz * 0.5
-
-    dz_stretch_level_index = nzt+1
-    dz_stretched = dz
-    DO  k = 2, nzt+1
-       IF ( dz_stretch_level <= zu(k-1)  .AND.  dz_stretched < dz_max )  THEN
-          dz_stretched = dz_stretched * dz_stretch_factor
-          dz_stretched = MIN( dz_stretched, dz_max )
-          IF ( dz_stretch_level_index == nzt+1 )  dz_stretch_level_index = k-1
-       ENDIF
-       zu(k) = zu(k-1) + dz_stretched
-    ENDDO
-
-!
-!-- Compute the u-levels. They are always staggered half-way between the 
-!-- corresponding w-levels. The top w-level is extrapolated linearly.
-    zw(0) = 0.0
-    DO  k = 1, nzt
-       zw(k) = ( zu(k) + zu(k+1) ) * 0.5
-    ENDDO
-    zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
+
+!
+!-- Define the vertical grid levels
+    IF ( .NOT. ocean )  THEN
+!
+!--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
+!--    Since the w-level lies on the surface, the first u-level (staggered!)
+!--    lies below the surface (used for "mirror" boundary condition).
+!--    The first u-level above the surface corresponds to the top of the 
+!--    Prandtl-layer.
+       zu(0) = - dz * 0.5
+       zu(1) =   dz * 0.5
+
+       dz_stretch_level_index = nzt+1
+       dz_stretched = dz
+       DO  k = 2, nzt+1
+          IF ( dz_stretch_level <= zu(k-1)  .AND.  dz_stretched < dz_max )  THEN
+             dz_stretched = dz_stretched * dz_stretch_factor
+             dz_stretched = MIN( dz_stretched, dz_max )
+             IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
+          ENDIF
+          zu(k) = zu(k-1) + dz_stretched
+       ENDDO
+
+!
+!--    Compute the w-levels. They are always staggered half-way between the 
+!--    corresponding u-levels. The top w-level is extrapolated linearly.
+       zw(0) = 0.0
+       DO  k = 1, nzt
+          zw(k) = ( zu(k) + zu(k+1) ) * 0.5
+       ENDDO
+       zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
+
+    ELSE
+!
+!--    Grid for ocean with solid surface at z=0 (k=0, w-grid). The free water
+!--    surface is at k=nzt (w-grid).
+!--    Since the w-level lies always on the surface, the first/last u-level
+!--    (staggered!) lies below the bottom surface / above the free surface.
+!--    It is used for "mirror" boundary condition.
+!--    The first u-level above the bottom surface corresponds to the top of the 
+!--    Prandtl-layer.
+       zu(nzt+1) =   dz * 0.5
+       zu(nzt)   = - dz * 0.5
+
+       dz_stretch_level_index = 0
+       dz_stretched = dz
+       DO  k = nzt-1, 0, -1
+          IF ( dz_stretch_level <= ABS( zu(k+1) )  .AND.  &
+               dz_stretched < dz_max )  THEN
+             dz_stretched = dz_stretched * dz_stretch_factor
+             dz_stretched = MIN( dz_stretched, dz_max )
+             IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
+          ENDIF
+          zu(k) = zu(k+1) - dz_stretched
+       ENDDO
+
+!
+!--    Compute the w-levels. They are always staggered half-way between the 
+!--    corresponding u-levels.
+!--    The top w-level (nzt+1) is not used but set for consistency, since
+!--    w and all scalar variables are defined up tp nzt+1.
+       zw(nzt+1) = dz
+       zw(nzt)   = 0.0
+       DO  k = 0, nzt
+          zw(k) = ( zu(k) + zu(k+1) ) * 0.5
+       ENDDO
+
+    ENDIF
 
 !
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 93)
+++ /palm/trunk/SOURCE/modules.f90	(revision 94)
@@ -5,5 +5,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +ocean, r, + salinity variables
+! defaults of .._vertical_gradient_levels changed from -1.0 to -9999999.9
 !
 ! Former revisions:
@@ -94,5 +95,5 @@
     REAL, DIMENSION(:), ALLOCATABLE ::                                         &
           ddzu, dd2zu, dzu, ddzw, dzw, km_damp_x, km_damp_y, l_grid, pt_init,  &
-          q_init, rdf, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,               &
+          q_init, rdf, sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,      &
           v_nzb_p1_for_vfc, zu, zw
 
@@ -101,9 +102,9 @@
 
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
-          qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, shf_1, shf_2,        &
-          tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
+          qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, saswst_1, shf_1,     &
+          shf_2, tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
 
     REAL, DIMENSION(:,:), POINTER ::                                           &
-          qsws, qsws_m, qswst, qswst_m, rif, rif_m, shf, shf_m, tswst,         &
+          qsws, qsws_m, qswst, qswst_m, rif, rif_m, saswst, shf, shf_m, tswst, &
           tswst_m, usws, usws_m, vsws, vsws_m
 
@@ -117,11 +118,11 @@
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
           e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1,     &
-          q_2, q_3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2,    &
-          w_1, w_2, w_3
+          q_2, q_3, ql_1, ql_2, r, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1, v_2,  &
+          v_3, vpt_1, vpt_2, w_1, w_2, w_3
 
     REAL, DIMENSION(:,:,:), POINTER ::                                         &
           e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql,    &
-          ql_c, te_m, tpt_m, tq_m, tu_m, tv_m, tw_m, u, u_m, u_p, v, v_m, v_p, &
-          vpt, vpt_m, w, w_m, w_p
+          ql_c, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, u_m,  &
+          u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
 
     REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
@@ -288,4 +289,5 @@
                 pt_vertical_gradient_level_ind(10) = -9999, &
                 q_vertical_gradient_level_ind(10) = -9999, &
+                sa_vertical_gradient_level_ind(10) = -9999, &
                 section(100,3), section_xy(100) = -9999, &
                 section_xz(100) = -9999, section_yz(100) = -9999, &
@@ -312,5 +314,5 @@
                 mg_switch_to_pe0 = .FALSE., &
                 netcdf_output = .FALSE., netcdf_64bit = .FALSE., &
-                netcdf_64bit_3d = .TRUE., &
+                netcdf_64bit_3d = .TRUE., ocean = .FALSE., &
                 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
                 outflow_s = .FALSE., passive_scalar = .FALSE., &
@@ -362,5 +364,5 @@
              rayleigh_damping_height = -1.0, residual_limit = 1.0E-4, &
              restart_time = 9999999.9, rif_max = 1.0, rif_min = -5.0, &
-             roughness_length = 0.1, simulated_time = 0.0, &
+             roughness_length = 0.1, sa_surface = 35.0, simulated_time = 0.0, &
              simulated_time_at_begin, sin_alpha_surface, &
              skip_time_data_output = 0.0, skip_time_data_output_av = 9999999.9,&
@@ -378,5 +380,6 @@
              time_do_sla = 0.0, time_dvrp = 0.0, time_prel = 0.0, &
              time_restart = 9999999.9, time_run_control = 0.0, &
-             top_heatflux = 9999999.9, ug_surface = 0.0, u_gtrans = 0.0, &
+             top_heatflux = 9999999.9, top_salinityflux = 0.0, &
+             ug_surface = 0.0, u_gtrans = 0.0, &
              ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, &
              ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, &
@@ -386,9 +389,11 @@
     REAL ::  do2d_xy_last_time(0:1) = -1.0, do2d_xz_last_time(0:1) = -1.0, &
              do2d_yz_last_time(0:1) = -1.0, pt_vertical_gradient(10) = 0.0, &
-             pt_vertical_gradient_level(10) = -1.0, &
+             pt_vertical_gradient_level(10) = -9999999.9, &
              q_vertical_gradient(10) = 0.0, &
              q_vertical_gradient_level(10) = -1.0, &
              s_vertical_gradient(10) = 0.0, &
-             s_vertical_gradient_level(10) = -1.0, threshold(20) = 0.0, &
+             s_vertical_gradient_level(10) = -1.0,
+             sa_vertical_gradient(10) = 0.0, &
+             sa_vertical_gradient_level(10) = -1.0, threshold(20) = 0.0, &
              tsc(10) = (/ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), &
              ug_vertical_gradient(10) = 0.0, &
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 93)
+++ /palm/trunk/SOURCE/parin.f90	(revision 94)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level,
+! top_salinity_flux in inipar,
+! sa_init is allocated
 !
 ! Former revisions:
@@ -61,40 +63,34 @@
     NAMELIST /inipar/  adjust_mixing_length, alpha_surface, bc_e_b, bc_lr, &
                        bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b, &
-                       bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t, &
-                       building_height, building_length_x, building_length_y, &
-                       building_wall_left, building_wall_south, &
-                       cloud_droplets, cloud_physics, conserve_volume_flow, &
-                       cut_spline_overshoot, damp_level_1d, dissipation_1d, &
-                       dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, &
-                       dz_stretch_factor, dz_stretch_level, e_min, &
-                       end_time_1d, fft_method, galilei_transformation, &
-                       grid_matching, humidity, inflow_disturbance_begin, &
-                       inflow_disturbance_end, initializing_actions, &
-                       km_constant, km_damp_max, long_filter_factor, &
-                       loop_optimization, mixing_length_1d, &
-                       momentum_advec, netcdf_precision, npex, npey, nsor_ini, &
-                       nx, ny, nz, omega, outflow_damping_width, &
-                       overshoot_limit_e, &
-                       overshoot_limit_pt, overshoot_limit_u, &
-                       overshoot_limit_v, overshoot_limit_w, passive_scalar, &
-                       phi, prandtl_layer, precipitation, pt_reference, &
-                       pt_surface, pt_surface_initial_change, &
-                       pt_vertical_gradient, pt_vertical_gradient_level, &
-                       q_surface, q_surface_initial_change, &
-                       q_vertical_gradient, q_vertical_gradient_level, &
-                       radiation, random_generator, random_heatflux, rif_max, &
-                       rif_min, roughness_length, scalar_advec, &
-                       statistic_regions, surface_heatflux, surface_pressure, &
-                       surface_scalarflux, surface_waterflux, s_surface, &
-                       s_surface_initial_change, s_vertical_gradient, &
-                       s_vertical_gradient_level, top_heatflux, &
-                       timestep_scheme, topography, ug_surface, &
-                       ug_vertical_gradient, ug_vertical_gradient_level, &
-                       ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, &
-                       ups_limit_w, use_surface_fluxes, use_top_fluxes, &
-                       use_ug_for_galilei_tr, use_upstream_for_tke, &
-                       vg_surface, vg_vertical_gradient, &
-                       vg_vertical_gradient_level, wall_adjustment, &
-                       wall_heatflux
+             bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t, building_height, &
+             building_length_x, building_length_y, building_wall_left, &
+             building_wall_south, cloud_droplets, cloud_physics, &
+             conserve_volume_flow, cut_spline_overshoot, damp_level_1d, &
+             dissipation_1d, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, &
+             dz_max, dz_stretch_factor, dz_stretch_level, e_min, end_time_1d, &
+             fft_method, galilei_transformation, grid_matching, humidity, &
+             inflow_disturbance_begin, inflow_disturbance_end, &
+             initializing_actions, km_constant, km_damp_max, &
+             long_filter_factor, loop_optimization, mixing_length_1d, &
+             momentum_advec, netcdf_precision, npex, npey, nsor_ini, nx, ny, &
+             nz, ocean, omega, outflow_damping_width, overshoot_limit_e, &
+             overshoot_limit_pt, overshoot_limit_u, overshoot_limit_v, &
+             overshoot_limit_w, passive_scalar, phi, prandtl_layer, &
+             precipitation, pt_reference, pt_surface, &
+             pt_surface_initial_change, pt_vertical_gradient, &
+             pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
+             q_vertical_gradient, q_vertical_gradient_level, radiation, &
+             random_generator, random_heatflux, rif_max, rif_min, &
+             roughness_length, sa_surface, sa_vertical_gradient, &
+             sa_vertical_gradient_level, scalar_advec, statistic_regions, &
+             surface_heatflux, surface_pressure, surface_scalarflux, &
+             surface_waterflux, s_surface, s_surface_initial_change, &
+             s_vertical_gradient, s_vertical_gradient_level, top_heatflux, &
+             top_salinity_flux, timestep_scheme, topography, ug_surface, &
+             ug_vertical_gradient, ug_vertical_gradient_level, ups_limit_e, &
+             ups_limit_pt, ups_limit_u, ups_limit_v, ups_limit_w, &
+             use_surface_fluxes, use_top_fluxes, use_ug_for_galilei_tr, &
+             use_upstream_for_tke, vg_surface, vg_vertical_gradient, &
+             vg_vertical_gradient_level, wall_adjustment, wall_heatflux
 
 
@@ -232,6 +228,6 @@
 
        ALLOCATE( ug(0:nz+1), vg(0:nz+1), &
-                 pt_init(0:nz+1), q_init(0:nz+1), u_init(0:nz+1), &
-                 v_init(0:nz+1),                                  &
+                 pt_init(0:nz+1), q_init(0:nz+1), sa_init(0:nz+1), &
+                 u_init(0:nz+1), v_init(0:nz+1),                   &
                  hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions) )
        hom = 0.0
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 93)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 94)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! new argument zw in calls of diffusion_e
 !
 ! Former revisions:
@@ -557,8 +557,8 @@
                 IF ( .NOT. humidity )  THEN
                    CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                     l_grid, pt, rif, tend, zu )
+                                     l_grid, pt, rif, tend, zu, zw )
                 ELSE
                    CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                     l_grid, vpt, rif, tend, zu )
+                                     l_grid, vpt, rif, tend, zu, zw )
                 ENDIF
              ELSE
@@ -582,16 +582,18 @@
                    IF ( .NOT. humidity )  THEN
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
-                                        km_m, l_grid, pt_m, rif_m, tend, zu )
+                                        km_m, l_grid, pt_m, rif_m, tend, zu, &
+                                        zw )
                    ELSE
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
-                                        km_m, l_grid, vpt_m, rif_m, tend, zu )
+                                        km_m, l_grid, vpt_m, rif_m, tend, zu, &
+                                        zw )
                    ENDIF
                 ELSE
                    IF ( .NOT. humidity )  THEN
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, pt, rif, tend, zu )
+                                        l_grid, pt, rif, tend, zu, zw )
                    ELSE
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, vpt, rif, tend, zu )
+                                        l_grid, vpt, rif, tend, zu, zw )
                    ENDIF
                 ENDIF
@@ -961,16 +963,18 @@
                    IF ( .NOT. humidity )  THEN
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
-                                        km_m, l_grid, pt_m, rif_m, tend, zu )
+                                        km_m, l_grid, pt_m, rif_m, tend, zu, &
+                                        zw )
                    ELSE
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
-                                        km_m, l_grid, vpt_m, rif_m, tend, zu )
+                                        km_m, l_grid, vpt_m, rif_m, tend, zu, &
+                                        zw )
                    ENDIF
                 ELSE
                    IF ( .NOT. humidity )  THEN
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, pt, rif, tend, zu )
+                                        l_grid, pt, rif, tend, zu, zw )
                    ELSE
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, vpt, rif, tend, zu )
+                                        l_grid, vpt, rif, tend, zu, zw )
                    ENDIF
                 ENDIF
@@ -1500,8 +1504,8 @@
           IF ( .NOT. humidity )  THEN
              CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
-                               rif, tend, zu )
+                               rif, tend, zu, zw )
           ELSE
              CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
-                               rif, tend, zu )
+                               rif, tend, zu, zw )
           ENDIF
        ELSE
@@ -1523,16 +1527,16 @@
              IF ( .NOT. humidity )  THEN
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
-                                  pt_m, rif_m, tend, zu )
+                                  pt_m, rif_m, tend, zu, zw )
              ELSE
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
-                                  vpt_m, rif_m, tend, zu )
+                                  vpt_m, rif_m, tend, zu, zw )
              ENDIF
           ELSE
              IF ( .NOT. humidity )  THEN
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
-                                  rif, tend, zu )
+                                  rif, tend, zu, zw )
              ELSE
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
-                                  rif, tend, zu )
+                                  rif, tend, zu, zw )
              ENDIF
           ENDIF
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 93)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 94)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +sa, saswst
 !
 ! Former revisions:
@@ -282,4 +282,8 @@
              ALLOCATE( s_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
              READ ( 13 )  s_av
+          CASE ( 'sa' )
+             READ ( 13 )  sa
+          CASE ( 'saswst' )
+             READ ( 13 )  saswst
           CASE ( 'shf' )
              READ ( 13 )  shf
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 93)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 94)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +ocean, sa_init, sa_surface, sa_vertical_gradient,
+! sa_vertical_gradient_level, top_salinity_flux
 !
 ! Former revisions:
@@ -107,5 +108,5 @@
     READ ( 13 )  statistic_regions
     ALLOCATE( ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1), v_init(0:nz+1), & 
-              pt_init(0:nz+1), q_init(0:nz+1), &
+              pt_init(0:nz+1), q_init(0:nz+1), sa_init(0:nz+1), &
               hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions) )
 
@@ -256,4 +257,6 @@
           CASE ( 'ny' )
              READ ( 13 )  ny
+          CASE ( 'ocean' )
+             READ ( 13 )  ocean
           CASE ( 'old_dt' )
              READ ( 13 )  old_dt
@@ -320,4 +323,12 @@
           CASE ( 'runnr' )
              READ ( 13 )  runnr
+          CASE ( 'sa_init' )
+             READ ( 13 )  sa_init
+          CASE ( 'sa_surface' )
+             READ ( 13 )  sa_surface
+          CASE ( 'sa_vertical_gradient' )
+             READ ( 13 )  sa_vertical_gradient
+          CASE ( 'sa_vertical_gradient_level' )
+             READ ( 13 )  sa_vertical_gradient_level
           CASE ( 'scalar_advec' )
              READ ( 13 )  scalar_advec
@@ -378,4 +389,6 @@
           CASE ( 'top_heatflux' )
              READ ( 13 )  top_heatflux
+          CASE ( 'top_salinityflux' )
+             READ ( 13 )  top_salinityflux
           CASE ( 'tsc' )
              READ ( 13 )  tsc
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 93)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 94)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +sa, saswst
 !
 ! Former revisions:
@@ -131,4 +131,8 @@
           WRITE ( 14 )  'qswst_m             ';  WRITE ( 14 ) qswst_m
        ENDIF
+    ENDIF
+    IF ( ocean )  THEN
+       WRITE ( 14 )  'sa                  ';  WRITE ( 14 )  sa
+       WRITE ( 14 )  'saswst              ';  WRITE ( 14 )  saswst
     ENDIF
     IF ( ALLOCATED( ql_c_av ) )  THEN
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 93)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 94)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +ocean, sa_init, sa_surface, sa_vertical_gradient,
+! sa_vertical_gradient_level, top_salinity_flux
 !
 ! Former revisions:
@@ -202,4 +203,6 @@
     WRITE ( 14 )  'ny                            '
     WRITE ( 14 )  ny
+    WRITE ( 14 )  'ocean                         '
+    WRITE ( 14 )  ocean
     WRITE ( 14 )  'old_dt                        '
     WRITE ( 14 )  old_dt
@@ -266,4 +269,12 @@
     WRITE ( 14 )  'runnr                         '
     WRITE ( 14 )  runnr
+    WRITE ( 14 )  'sa_init                       '
+    WRITE ( 14 )  sa_init
+    WRITE ( 14 )  'sa_surface                    '
+    WRITE ( 14 )  sa_surface
+    WRITE ( 14 )  'sa_vertical_gradient          '
+    WRITE ( 14 )  sa_vertical_gradient
+    WRITE ( 14 )  'sa_vertical_gradient_level    '
+    WRITE ( 14 )  sa_vertical_gradient_level
     WRITE ( 14 )  'scalar_advec                  '
     WRITE ( 14 )  scalar_advec
@@ -324,4 +335,6 @@
     WRITE ( 14 )  'top_heatflux                  '
     WRITE ( 14 )  top_heatflux
+    WRITE ( 14 )  'top_salinityflux              '
+    WRITE ( 14 )  top_salinityflux
     WRITE ( 14 )  'tsc                           '
     WRITE ( 14 )  tsc
