Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 96)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 97)
@@ -1,17 +1,8 @@
 New:
 ---
-ocean version including prognostic equation for salinity and equation of state for seawater
-+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux
 
-average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, parin, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary
-
-New:
-eqn_state_seawater, init_ocean
 
 Changed:
 -------
-hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile
-
-advec_particles, buoyancy, calc_liquid_water_content, init_cloud_physics, modules, prognostic_equations
 
 
@@ -19,25 +10,3 @@
 ------
 
-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
-
-Bott-Chlond scheme to be extended for salinity
-
-Upstream-spline scheme to be extended for salinity?????
Index: palm/trunk/SOURCE/advec_s_bc.f90
===================================================================
--- palm/trunk/SOURCE/advec_s_bc.f90	(revision 96)
+++ palm/trunk/SOURCE/advec_s_bc.f90	(revision 97)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Advection of salinity included
+! Bugfix: Error in boundary condition for TKE removed
 !
 ! Former revisions:
@@ -755,29 +756,37 @@
        ENDIF
 
+    ELSEIF ( sk_char == 'sa' )  THEN
+
+!
+!--    Salinity boundary condition at the bottom boundary.
+!--    So far, always Neumann (i.e. here zero gradient) is used
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+          ENDDO
+       ENDDO
+
+!
+!--    Salinity boundary condition at the top boundary.
+!--    Dirichlet or Neumann (zero gradient)
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+             sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+          ENDDO
+       ENDDO
+
     ELSEIF ( sk_char == 'q' )  THEN
 
 !
-!--    Specific humidity boundary condition at the bottom boundary
-       IF ( ibc_q_b == 0 )  THEN
-!
-!--       Dirichlet (fixed surface humidity)
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
-                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
-             ENDDO
-          ENDDO
-
-       ELSE
-!
-!--       Neumann (i.e. here zero gradient)
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
-                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
-             ENDDO
-          ENDDO
-
-       ENDIF
+!--    Specific humidity boundary condition at the bottom boundary.
+!--    Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+          ENDDO
+       ENDDO
 
 !
@@ -809,8 +818,12 @@
 !
 !--    TKE boundary condition at bottom and top boundary (generally Neumann)
-       sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
-       sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
-       sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
-       sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+             sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
+             sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
+          ENDDO
+       ENDDO
 
     ELSE
Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 96)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 97)
@@ -4,4 +4,7 @@
 ! Actual revisions:
 ! -----------------
+! Routine reneralized to be used with temperature AND density:
+! argument theta renamed var, new argument var_reference,
+! use_pt_reference renamed use_reference,
 ! calc_mean_pt_profile renamed calc_mean_profile
 !
@@ -46,5 +49,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy( theta, wind_component, pr )
+    SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
 
        USE arrays_3d
@@ -57,5 +60,6 @@
 
        INTEGER ::  i, j, k, pr, wind_component
-       REAL, DIMENSION(:,:,:), POINTER  ::  theta
+       REAL    ::  var_reference
+       REAL, DIMENSION(:,:,:), POINTER  ::  var
 
 
@@ -63,11 +67,12 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_pt_reference )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                 &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
+          IF ( use_reference )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
+                                                            (                  &
+                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
+                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
                                                             )
                    ENDDO
@@ -78,7 +83,8 @@
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                  &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
+                                                            (                  &
+                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
+                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
                                                             )
                    ENDDO
@@ -136,5 +142,5 @@
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy_ij( i, j, theta, wind_component, pr )
+    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
 
        USE arrays_3d
@@ -147,5 +153,6 @@
 
        INTEGER ::  i, j, k, pr, wind_component
-       REAL, DIMENSION(:,:,:), POINTER  ::  theta
+       REAL    ::  var_reference
+       REAL, DIMENSION(:,:,:), POINTER  ::  var
 
 
@@ -153,17 +160,17 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_pt_reference )  THEN
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
-                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                      &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
-                                                       )
-             ENDDO
-          ELSE
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
-                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                       &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                       )
+          IF ( use_reference )  THEN
+             DO  k = nzb_s_inner(j,i)+1, nzt-1
+                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (   &
+                         ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
+                         ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
+                                                                          )
+             ENDDO
+          ELSE
+             DO  k = nzb_s_inner(j,i)+1, nzt-1
+                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (    &
+                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
+                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                                                                          )
              ENDDO
           ENDIF
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 96)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 97)
@@ -4,7 +4,9 @@
 ! Actual revisions:
 ! -----------------
-! Initial salinity profile is calculated.
+! Initial salinity profile is calculated, salinity boundary conditions are
+! checked,
 ! z_max_do1d is checked only in case of ocean = .f.,
-! +initial temperature profile for the ocean version,
+! +initial temperature and geostrophic velocity profiles for the ocean version,
+! use_pt_reference renamed use_reference
 !
 ! Former revisions:
@@ -163,4 +165,7 @@
           WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
        ENDIF
+       IF ( momentum_advec == 'ups-scheme' )  THEN
+          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
+       ENDIF
        IF ( action /= ' ' )  THEN
           IF ( myid == 0 )  THEN
@@ -447,30 +452,64 @@
        i = 1
        gradient = 0.0
-       ug_vertical_gradient_level_ind(1) = 0
-       ug(0) = ug_surface
-       DO  k = 1, nzt+1
-          IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
-               ug_vertical_gradient_level(i) >= 0.0 )  THEN
-             gradient = ug_vertical_gradient(i) / 100.0
-             ug_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', &
-                           ' "ug_vertical_gradient_level_ind" exceeded'
-                ENDIF
-                CALL local_stop
-             ENDIF
-          ENDIF
-          IF ( gradient /= 0.0 )  THEN
-             IF ( k /= 1 )  THEN
-                ug(k) = ug(k-1) + dzu(k) * gradient
+
+       IF ( .NOT. ocean )  THEN
+
+          ug_vertical_gradient_level_ind(1) = 0
+          ug(0) = ug_surface
+          DO  k = 1, nzt+1
+             IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
+                  ug_vertical_gradient_level(i) >= 0.0 )  THEN
+                gradient = ug_vertical_gradient(i) / 100.0
+                ug_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', &
+                              ' "ug_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= 1 )  THEN
+                   ug(k) = ug(k-1) + dzu(k) * gradient
+                ELSE
+                   ug(k) = ug_surface + 0.5 * dzu(k) * gradient
+                ENDIF
              ELSE
-                ug(k) = ug_surface + 0.5 * dzu(k) * gradient
-             ENDIF
-          ELSE
-             ug(k) = ug(k-1)
-          ENDIF
-       ENDDO
+                ug(k) = ug(k-1)
+             ENDIF
+          ENDDO
+
+       ELSE
+
+          ug_vertical_gradient_level_ind(1) = nzt+1
+          DO  k = nzt, 0, -1
+             IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
+                  ug_vertical_gradient_level(i) <= 0.0 )  THEN
+                gradient = ug_vertical_gradient(i) / 100.0
+                ug_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', &
+                              ' "ug_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= nzt )  THEN
+                   ug(k) = ug(k+1) - dzu(k+1) * gradient
+                ELSE
+                   ug(k)   = ug_surface - 0.5 * dzu(k+1) * gradient
+                   ug(k+1) = ug_surface + 0.5 * dzu(k+1) * gradient
+                ENDIF
+             ELSE
+                ug(k) = ug(k+1)
+             ENDIF
+          ENDDO
+
+       ENDIF
 
        u_init = ug
@@ -478,5 +517,5 @@
 !
 !--    In case of no given gradients for ug, choose a vanishing gradient
-       IF ( ug_vertical_gradient_level(1) == -1.0 )  THEN
+       IF ( ug_vertical_gradient_level(1) == -9999999.9 )  THEN
           ug_vertical_gradient_level(1) = 0.0
        ENDIF  
@@ -488,30 +527,64 @@
        i = 1
        gradient = 0.0
-       vg_vertical_gradient_level_ind(1) = 0
-       vg(0) = vg_surface
-       DO  k = 1, nzt+1
-          IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
-               vg_vertical_gradient_level(i) >= 0.0 )  THEN
-             gradient = vg_vertical_gradient(i) / 100.0
-             vg_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', &
-                           ' "vg_vertical_gradient_level_ind" exceeded'
-                ENDIF
-                CALL local_stop
-             ENDIF
-          ENDIF
-          IF ( gradient /= 0.0 )  THEN
-             IF ( k /= 1 )  THEN
-                vg(k) = vg(k-1) + dzu(k) * gradient
+
+       IF ( .NOT. ocean )  THEN
+
+          vg_vertical_gradient_level_ind(1) = 0
+          vg(0) = vg_surface
+          DO  k = 1, nzt+1
+             IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
+                  vg_vertical_gradient_level(i) >= 0.0 )  THEN
+                gradient = vg_vertical_gradient(i) / 100.0
+                vg_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', &
+                              ' "vg_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= 1 )  THEN
+                   vg(k) = vg(k-1) + dzu(k) * gradient
+                ELSE
+                   vg(k) = vg_surface + 0.5 * dzu(k) * gradient
+                ENDIF
              ELSE
-                vg(k) = vg_surface + 0.5 * dzu(k) * gradient
-             ENDIF
-          ELSE
-             vg(k) = vg(k-1)
-          ENDIF
-       ENDDO
+                vg(k) = vg(k-1)
+             ENDIF
+          ENDDO
+
+       ELSE
+
+          vg_vertical_gradient_level_ind(1) = 0
+          DO  k = nzt, 0, -1
+             IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
+                  vg_vertical_gradient_level(i) <= 0.0 )  THEN
+                gradient = vg_vertical_gradient(i) / 100.0
+                vg_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', &
+                              ' "vg_vertical_gradient_level_ind" exceeded'
+                   ENDIF
+                   CALL local_stop
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0 )  THEN
+                IF ( k /= nzt )  THEN
+                   vg(k) = vg(k+1) - dzu(k+1) * gradient
+                ELSE
+                   vg(k)   = vg_surface - 0.5 * dzu(k+1) * gradient
+                   vg(k+1) = vg_surface + 0.5 * dzu(k+1) * gradient
+                ENDIF
+             ELSE
+                vg(k) = vg(k+1)
+             ENDIF
+          ENDDO
+
+       ENDIF
 
        v_init = vg
@@ -519,5 +592,5 @@
 !
 !--    In case of no given gradients for vg, choose a vanishing gradient
-       IF ( vg_vertical_gradient_level(1) == -1.0 )  THEN
+       IF ( vg_vertical_gradient_level(1) == -9999999.9 )  THEN
           vg_vertical_gradient_level(1) = 0.0
        ENDIF
@@ -707,6 +780,19 @@
 
 !
-!-- Reference temperature to be used in buoyancy terms
-    IF ( pt_reference /= 9999999.9 )  use_pt_reference = .TRUE.
+!-- Ocean runs always use reference values in the buoyancy term. Therefore
+!-- set the reference temperature equal to the surface temperature.
+    IF ( ocean  .AND.  pt_reference == 9999999.9 )  pt_reference = pt_surface
+
+!
+!-- Reference value has to be used in buoyancy terms
+    IF ( pt_reference /= 9999999.9 )  use_reference = .TRUE.
+
+!
+!-- Sign of buoyancy/stability terms
+    IF ( ocean )  atmos_ocean_sign = -1.0
+
+!
+!-- Ocean version is using flux boundary conditions at the top
+    IF ( ocean )  use_top_fluxes = .TRUE.
 
 !
@@ -1020,4 +1106,12 @@
 
        IF ( top_salinityflux == 9999999.9 )  constant_top_salinityflux = .FALSE.
+       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9 )  THEN
+          IF ( myid == 0 )  THEN
+             PRINT*, '+++ check_parameters:'
+             PRINT*, '    boundary_condition: bc_sa_t = ', bc_sa_t
+             PRINT*, '    requires to set top_salinityflux '
+          ENDIF
+          CALL local_stop
+       ENDIF
 
 !
@@ -2474,7 +2568,12 @@
 !
 !-- Determine upper and lower hight level indices for random perturbations
-    IF ( disturbance_level_b == -1.0 )  THEN
-       disturbance_level_b     = zu(nzb+3)
-       disturbance_level_ind_b = nzb + 3
+    IF ( disturbance_level_b == -9999999.9 )  THEN
+       IF ( ocean ) THEN
+          disturbance_level_b     = zu((nzt*2)/3)
+          disturbance_level_ind_b = ( nzt * 2 ) / 3
+       ELSE
+          disturbance_level_b     = zu(nzb+3)
+          disturbance_level_ind_b = nzb + 3
+       ENDIF
     ELSEIF ( disturbance_level_b < zu(3) )  THEN
        IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
@@ -2496,7 +2595,12 @@
     ENDIF
 
-    IF ( disturbance_level_t == -1.0 )  THEN
-       disturbance_level_t     = zu(nzt/3)
-       disturbance_level_ind_t = nzt / 3
+    IF ( disturbance_level_t == -9999999.9 )  THEN
+       IF ( ocean )  THEN
+          disturbance_level_t     = zu(nzt-3)
+          disturbance_level_ind_t = nzt - 3
+       ELSE
+          disturbance_level_t     = zu(nzt/3)
+          disturbance_level_ind_t = nzt / 3
+       ENDIF
     ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
        IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 96)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 97)
@@ -8,4 +8,6 @@
 ! This is also a bugfix, because the height above the topography is now
 ! used instead of the height above level k=0.
+! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
+! use_pt_reference renamed use_reference
 !
 ! Former revisions:
@@ -50,6 +52,6 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
-                            rif, tend, zu, zw )
+    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, var, &
+                            var_reference, rif, tend, zu, zw )
 
        USE control_parameters
@@ -61,10 +63,10 @@
 
        INTEGER ::  i, j, k
-       REAL            ::  dpt_dz, l_stable, phi_m
+       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
        REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1: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
-       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
+       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
        REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
  
@@ -73,5 +75,5 @@
 !--    This if clause must be outside the k-loop because otherwise
 !--    runtime errors occur with -C hopt on NEC
-       IF ( use_pt_reference )  THEN
+       IF ( use_reference )  THEN
 
           DO  i = nxl, nxr
@@ -91,8 +93,9 @@
 !
 !--                Calculate the mixing length (for dissipation)
-                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
-                   IF ( dpt_dz > 0.0 ) THEN
+                   dvar_dz = atmos_ocean_sign * &
+                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                   IF ( dvar_dz > 0.0 ) THEN
                       l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                        SQRT( g / pt_reference * dpt_dz ) + 1E-5
+                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
                    ELSE
                       l_stable = l_grid(k)
@@ -180,8 +183,9 @@
 !
 !--                Calculate the mixing length (for dissipation)
-                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
-                   IF ( dpt_dz > 0.0 ) THEN
+                   dvar_dz = atmos_ocean_sign * &
+                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                   IF ( dvar_dz > 0.0 ) THEN
                       l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                                        SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
                    ELSE
                       l_stable = l_grid(k)
@@ -270,5 +274,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
-                               theta, rif, tend, zu, zw )
+                               var, var_reference, rif, tend, zu, zw )
 
        USE control_parameters
@@ -280,10 +284,10 @@
 
        INTEGER         ::  i, j, k
-       REAL            ::  dpt_dz, l_stable, phi_m
+       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
        REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1: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
-       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
+       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
        REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
 
@@ -303,12 +307,13 @@
 !--    Calculate the mixing length (for dissipation)
        DO  k = nzb_s_inner(j,i)+1, nzt
-          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
-          IF ( dpt_dz > 0.0 ) THEN
-             IF ( use_pt_reference )  THEN
+          dvar_dz = atmos_ocean_sign * &
+                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+          IF ( dvar_dz > 0.0 ) THEN
+             IF ( use_reference )  THEN
                 l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                  SQRT( g / pt_reference * dpt_dz ) + 1E-5
+                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
              ELSE
                 l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
              ENDIF
           ELSE
Index: palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- palm/trunk/SOURCE/diffusivities.f90	(revision 96)
+++ palm/trunk/SOURCE/diffusivities.f90	(revision 97)
@@ -1,3 +1,3 @@
- SUBROUTINE diffusivities( theta )
+ SUBROUTINE diffusivities( var, var_reference )
 
 !------------------------------------------------------------------------------!
@@ -7,4 +7,6 @@
 ! This is also a bugfix, because the height above the topography is now
 ! used instead of the height above level k=0.
+! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
+! use_pt_reference renamed use_reference
 ! 
 ! Former revisions:
@@ -41,7 +43,9 @@
     INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
 
-    REAL    ::  dpt_dz, l_stable, phi_m = 1.0
+    REAL    ::  dvar_dz, l_stable, var_reference
 
-    REAL    ::  theta(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+    REAL, SAVE ::  phi_m = 1.0
+
+    REAL    ::  var(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
 
     REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
@@ -58,5 +62,5 @@
 !
 !-- Compute the turbulent diffusion coefficient for momentum
-    !$OMP PARALLEL PRIVATE (dpt_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
+    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
 !$  tn = omp_get_thread_num()
 
@@ -94,12 +98,13 @@
 !--       Determine the mixing length
           DO  k = nzb_s_inner(j,i)+1, nzt
-             dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
-             IF ( dpt_dz > 0.0 ) THEN
-                IF ( use_pt_reference )  THEN
+             dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
+                       ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+             IF ( dvar_dz > 0.0 ) THEN
+                IF ( use_reference )  THEN
                    l_stable = 0.76 * sqrt_e(k) / &
-                                     SQRT( g / pt_reference * dpt_dz ) + 1E-5
+                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
                 ELSE
                    l_stable = 0.76 * sqrt_e(k) / &
-                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
                 ENDIF
              ELSE
Index: palm/trunk/SOURCE/eqn_state_seawater.f90
===================================================================
--- palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 96)
+++ palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 97)
@@ -8,5 +8,5 @@
 ! Former revisions:
 ! -----------------
-! $Id:$
+! $Id$
 !
 ! Initial revision
@@ -18,4 +18,6 @@
 ! salinity, and pressure.
 ! For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
+! eqn_state_seawater calculates the potential density referred at hyp(0).
+! eqn_state_seawater_func calculates density.
 !------------------------------------------------------------------------------!
 
@@ -70,8 +72,8 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_u_inner(j,i)+1, nzt
+             DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--             Pressure is needed in dbar
-                p1 = hyp(k) * 1E-4
+                p1 = hyp(0) * 1E-4
                 p2 = p1 * p1
                 p3 = p2 * p1
@@ -102,4 +104,8 @@
 
              ENDDO
+!
+!--          Neumann conditions are assumed at bottom and top boundary
+             rho(nzt+1,j,i)            = rho(nzt,j,i)
+             rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
           ENDDO
        ENDDO
@@ -122,8 +128,8 @@
        REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
-       DO  k = nzb_u_inner(j,i)+1, nzt
+       DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--       Pressure is needed in dbar
-          p1 = hyp(k) * 1E-4
+          p1 = hyp(0) * 1E-4
           p2 = p1 * p1
           p3 = p2 * p1
@@ -152,4 +158,8 @@
                        )
        ENDDO
+!
+!--    Neumann conditions are assumed at bottom and top boundary
+       rho(nzt+1,j,i)            = rho(nzt,j,i)
+       rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
 
     END SUBROUTINE eqn_state_seawater_ij
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 96)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 97)
@@ -5,6 +5,6 @@
 ! -----------------
 ! Statistics for ocean version (salinity, density) added,
-! calculation of Deardorff velocity scale adjusted to be used with the ocean
-! version (HAS STILL TO BE COMPLETED!!!)
+! calculation of z_i and Deardorff velocity scale adjusted to be used with
+! the ocean version
 !
 ! Former revisions:
@@ -558,5 +558,5 @@
 !
 !--             Salinity flux and density (density does not belong to here,
-!--             but so far there is no suitable place to calculate)
+!--             but so far there is no other suitable place to calculate)
                 IF ( ocean )  THEN
                    pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
@@ -611,4 +611,8 @@
 
 !
+!--    Density at top follows Neumann condition
+       IF ( ocean )  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
+
+!
 !--    Divergence of vertical flux of resolved scale energy and pressure
 !--    fluctuations. First calculate the products, then the divergence.
@@ -709,10 +713,10 @@
 !
 !--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
-          sums(nzb,58) = 0.0
-          sums(nzb,59) = 0.0
-          sums(nzb,60) = 0.0
-          sums(nzb,61) = 0.0
-          sums(nzb,62) = 0.0
-          sums(nzb,63) = 0.0
+          sums_l(nzb,58,tn) = 0.0
+          sums_l(nzb,59,tn) = 0.0
+          sums_l(nzb,60,tn) = 0.0
+          sums_l(nzb,61,tn) = 0.0
+          sums_l(nzb,62,tn) = 0.0
+          sums_l(nzb,63,tn) = 0.0
 
        ENDIF
@@ -857,5 +861,21 @@
        z_i(1) = 0.0
        first = .TRUE.
-!       IF ( .NOT. ocean )  THEN
+       IF ( ocean )  THEN
+          DO  k = nzt, nzb+1, -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
+                ENDIF
+                EXIT
+             ENDIF
+          ENDDO
+       ELSE
           DO  k = nzb, nzt-1
              IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
@@ -873,19 +893,28 @@
              ENDIF
           ENDDO
-!       ELSE
-!       ENDIF
-
-!
-!--    Second scheme: Starting from the top model boundary, look for the first 
-!--    characteristic kink in the temperature profile, where the originally 
-!--    stable stratification notably weakens.
+       ENDIF
+
+!
+!--    Second scheme: Starting from the top/bottom model boundary, look for
+!--    the first characteristic kink in the temperature profile, where the
+!--    originally stable stratification notably weakens.
        z_i(2) = 0.0
-       DO  k = nzt-1, nzb+1, -1
-          IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
-               2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
-             z_i(2) = zu(k)
-             EXIT
-          ENDIF
-       ENDDO
+       IF ( ocean )  THEN
+          DO  k = nzb+1, nzt-1
+             IF ( ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) > &
+                  2.0 * ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) )  THEN
+                z_i(2) = zu(k)
+                EXIT
+             ENDIF
+          ENDDO
+       ELSE
+          DO  k = nzt-1, nzb+1, -1
+             IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
+                  2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
+                z_i(2) = zu(k)
+                EXIT
+             ENDIF
+          ENDDO
+       ENDIF
 
        hom(nzb+6,1,pr_palm,sr) = z_i(1)
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 96)
+++ palm/trunk/SOURCE/header.f90	(revision 97)
@@ -4,6 +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. 
+! Adjustments for the ocean version.
+! use_pt_reference renamed use_reference
 !
 ! Former revisions:
@@ -103,4 +103,9 @@
           PRINT*,'+++ header:  unknown action(s): ',initializing_actions
        ENDIF
+    ENDIF
+    IF ( ocean )  THEN
+       run_classification = 'ocean - ' // run_classification
+    ELSE
+       run_classification = 'atmosphere - ' // run_classification
     ENDIF
 
@@ -391,35 +396,44 @@
        roben = 'e(nzt+1) = e(nzt) = e(nzt-1)'
 
-       WRITE ( io, 301 )  runten, roben       
-
-    ENDIF
-
-    IF ( humidity  .OR.  passive_scalar )  THEN
-       IF ( humidity )  THEN
-          IF ( ibc_q_b == 0 )  THEN
-             runten = 'q(0)     = q_surface'
-          ELSE
-             runten = 'q(0)     = q(1)'
-          ENDIF
-          IF ( ibc_q_t == 0 )  THEN
-             roben =  'q(nzt)   = q_top'
-          ELSE
-             roben =  'q(nzt)   = q(nzt-1) + dq/dz'
-          ENDIF
-       ELSE
-          IF ( ibc_q_b == 0 )  THEN
-             runten = 's(0)     = s_surface'
-          ELSE
-             runten = 's(0)     = s(1)'
-          ENDIF
-          IF ( ibc_q_t == 0 )  THEN
-             roben =  's(nzt)   = s_top'
-          ELSE
-             roben =  's(nzt)   = s(nzt-1) + ds/dz'
-          ENDIF
-       ENDIF
-
-       WRITE ( io, 302 ) runten, roben
-
+       WRITE ( io, 301 )  'e', runten, roben       
+
+    ENDIF
+
+    IF ( ocean )  THEN
+       runten = 'sa(0)    = sa(1)'
+       IF ( ibc_sa_t == 0 )  THEN
+          roben =  'sa(nzt+1) = sa_surface'
+       ELSE
+          roben =  'sa(nzt+1) = sa(nzt)'
+       ENDIF
+       WRITE ( io, 301 ) 'sa', runten, roben
+    ENDIF
+
+    IF ( humidity )  THEN
+       IF ( ibc_q_b == 0 )  THEN
+          runten = 'q(0)     = q_surface'
+       ELSE
+          runten = 'q(0)     = q(1)'
+       ENDIF
+       IF ( ibc_q_t == 0 )  THEN
+          roben =  'q(nzt)   = q_top'
+       ELSE
+          roben =  'q(nzt)   = q(nzt-1) + dq/dz'
+       ENDIF
+       WRITE ( io, 301 ) 'q', runten, roben
+    ENDIF
+
+    IF ( passive_scalar )  THEN
+       IF ( ibc_q_b == 0 )  THEN
+          runten = 's(0)     = s_surface'
+       ELSE
+          runten = 's(0)     = s(1)'
+       ENDIF
+       IF ( ibc_q_t == 0 )  THEN
+          roben =  's(nzt)   = s_top'
+       ELSE
+          roben =  's(nzt)   = s(nzt-1) + ds/dz'
+       ENDIF
+       WRITE ( io, 301 ) 's', runten, roben
     ENDIF
 
@@ -442,4 +456,7 @@
        IF ( constant_top_heatflux )  THEN
           WRITE ( io, 306 )  top_heatflux
+       ENDIF
+       IF ( ocean  .AND.  constant_top_salinityflux )  THEN
+          WRITE ( io, 309 )  top_salinityflux
        ENDIF
        IF ( humidity  .OR.  passive_scalar )  THEN
@@ -887,5 +904,11 @@
 !-- Other quantities
     WRITE ( io, 411 )  g
-    IF ( use_pt_reference )  WRITE ( io, 412 )  pt_reference
+    IF ( use_reference )  THEN
+       IF ( ocean )  THEN
+          WRITE ( io, 412 )  prho_reference
+       ELSE
+          WRITE ( io, 413 )  pt_reference
+       ENDIF
+    ENDIF
 
 !
@@ -1011,4 +1034,34 @@
                              TRIM( gradients ), TRIM( slices )
        ENDIF
+    ENDIF
+
+!
+!-- Initial salinity profile
+!-- Building output strings, starting with surface salinity
+    IF ( ocean )  THEN
+       WRITE ( temperatures, '(F6.2)' )  sa_surface
+       gradients = '------'
+       slices = '     0'
+       coordinates = '   0.0'
+       i = 1
+       DO  WHILE ( sa_vertical_gradient_level_ind(i) /= -9999 )
+
+          WRITE (coor_chr,'(F7.2)')  sa_init(sa_vertical_gradient_level_ind(i))
+          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F7.2)')  sa_vertical_gradient(i)
+          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I7)')  sa_vertical_gradient_level_ind(i)
+          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F7.1)')  sa_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          i = i + 1
+       ENDDO
+
+       WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
+                          TRIM( gradients ), TRIM( slices )
     ENDIF
 
@@ -1135,7 +1188,7 @@
 
  99 FORMAT (1X,78('-'))
-100 FORMAT (/1X,'*************************',11X,28('-')/        &
+100 FORMAT (/1X,'*************************',11X,42('-')/        &
             1X,'* ',A,' *',11X,A/                               &
-            1X,'*************************',11X,28('-')//        &
+            1X,'*************************',11X,42('-')//        &
             ' Date:            ',A8,11X,'Run:       ',A20/      &
             ' Time:            ',A8,11X,'Run-No.:   ',I2.2/     &
@@ -1256,8 +1309,5 @@
              ' B. bound.: ',A/ &
              ' T. bound.: ',A)
-301 FORMAT (/'                     e'// &
-             ' B. bound.: ',A/ &
-             ' T. bound.: ',A)
-302 FORMAT (/'                     q'// &
+301 FORMAT (/'                     ',A// &
              ' B. bound.: ',A/ &
              ' T. bound.: ',A)
@@ -1268,7 +1318,8 @@
              '       zp = ',F6.2,' m   z0 = ',F6.4,' m   kappa = ',F4.2/ &
              '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
-306 FORMAT ('       Predefined constant heatflux:   ',F6.3,' K m/s')
+306 FORMAT ('       Predefined constant heatflux:   ',F9.6,' K m/s')
 307 FORMAT ('       Heatflux has a random normal distribution')
 308 FORMAT ('       Predefined surface temperature')
+309 FORMAT ('       Predefined constant salinityflux:   ',F9.6,' psu m/s')
 310 FORMAT (//'    1D-Model:'// &
              '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
@@ -1352,5 +1403,6 @@
             '                            f*    = ',F9.6,' 1/s')
 411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
-412 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
+412 FORMAT (/'    Reference density in buoyancy terms: ',F8.3,' kg/m**3')
+413 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
 415 FORMAT (/'    Cloud physics parameters:'/ &
              '    ------------------------'/)
@@ -1382,7 +1434,12 @@
 424 FORMAT (/'    Characteristic levels of the geo. wind component vg:'// &
             '       Height:      ',A,'  m'/ &
-            '       vg:          ',A,'  m/S'/ &
+            '       vg:          ',A,'  m/s'/ &
             '       Gradient:    ',A,'  1/100s'/ &
             '       Gridpoint:   ',A)
+425 FORMAT (/'    Characteristic levels of the initial salinity profile:'// &
+            '       Height:     ',A,'  m'/ &
+            '       Salinity:   ',A,'  psu'/ &
+            '       Gradient:   ',A,'  psu/100m'/ &
+            '       Gridpoint:  ',A)
 450 FORMAT (//' LES / Turbulence quantities:'/ &
               ' ---------------------------'/)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 96)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 97)
@@ -529,4 +529,9 @@
        hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
 
+       IF ( ocean )  THEN
+!
+!--       Store initial salinity profile
+          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
+       ENDIF
 
        IF ( humidity )  THEN
Index: palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- palm/trunk/SOURCE/init_ocean.f90	(revision 96)
+++ palm/trunk/SOURCE/init_ocean.f90	(revision 97)
@@ -40,6 +40,6 @@
     hyp(nzt+1) = surface_pressure * 100.0
 
-    hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
-    rho_ref    = rho_surface * 0.5 * dzu(nzt+1)
+    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
+    rho_reference = rho_surface * 0.5 * dzu(nzt+1)
 
     DO  k = nzt-1, 0, -1
@@ -48,13 +48,27 @@
        pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
 
-       rho_l   = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
+       rho_l = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
 
-       hyp(k)  = hyp(k+1) + rho_l * g * dzu(k+1)
-       rho_ref = rho_ref + rho_l * dzu(k+1)
+       hyp(k)        = hyp(k+1) + rho_l * g * dzu(k+1)
+       rho_reference = rho_reference + rho_l * dzu(k+1)
 
     ENDDO
 
-    rho_ref = rho_ref / ( zw(nzt) - zu(nzb) )
-    print*, '*** rho_ref = ', rho_ref
+    rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
+
+!
+!-- Calculate the reference potential density
+    prho_reference = 0.0
+    DO  k = 0, nzt
+
+       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
+       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
+
+       prho_reference = prho_reference + dzu(k+1) * &
+                        eqn_state_seawater_func( hyp(0), pt_l, sa_l )
+
+    ENDDO
+
+    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
 
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 96)
+++ palm/trunk/SOURCE/modules.f90	(revision 97)
@@ -5,7 +5,7 @@
 ! Actual revisions:
 ! -----------------
-! +ocean, r, + salinity variables
+! +atmos_ocean_sign, ocean, r, + salinity variables
 ! defaults of .._vertical_gradient_levels changed from -1.0 to -9999999.9
-! hydro_press renamed hyp
+! hydro_press renamed hyp, use_pt_reference renamed use_reference
 !
 ! Former revisions:
@@ -324,5 +324,5 @@
                 sloping_surface = .FALSE., stop_dt = .FALSE., &
                 terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,&
-                use_pt_reference = .FALSE., use_surface_fluxes = .FALSE., &
+                use_reference = .FALSE., use_surface_fluxes = .FALSE., &
                 use_top_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., &
                 use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE.
@@ -333,4 +333,5 @@
     REAL ::  advected_distance_x = 0.0, advected_distance_y = 0.0, &
              alpha_surface = 0.0, asselin_filter_factor = 0.1, &
+             atmos_ocean_sign = 1.0, &
              averaging_interval = 0.0, averaging_interval_pr = 9999999.9, &
              averaging_interval_sp = 9999999.9, bc_pt_t_val, bc_q_t_val, &
@@ -340,6 +341,8 @@
              building_wall_south = 9999999.9, cfl_factor = -1.0, &
              cos_alpha_surface, disturbance_amplitude = 0.25, &
-             disturbance_energy_limit = 0.01, disturbance_level_b = -1.0, &
-             disturbance_level_t = -1.0, dt = -1.0, dt_averaging_input = 0.0, &
+             disturbance_energy_limit = 0.01, &
+             disturbance_level_b = -9999999.9, &
+             disturbance_level_t = -9999999.9, &
+             dt = -1.0, dt_averaging_input = 0.0, &
              dt_averaging_input_pr = 9999999.9, dt_data_output = 9999999.9, &
              dt_data_output_av = 9999999.9, dt_disturb = 9999999.9, &
@@ -360,11 +363,11 @@
              overshoot_limit_w = 0.0, particle_maximum_age = 9999999.9, &
              phi = 55.0, prandtl_number = 1.0, &
-             precipitation_amount_interval = 9999999.9, &
-             pt_reference = 9999999.9, &
-             pt_slope_offset = 0.0, pt_surface = 300.0, &
-             pt_surface_initial_change = 0.0, q_surface = 0.0, &
-             q_surface_initial_change = 0.0, rayleigh_damping_factor = -1.0, &
-             rayleigh_damping_height = -1.0, residual_limit = 1.0E-4, &
-             restart_time = 9999999.9, rho_ref, rho_surface, rif_max = 1.0, &
+             precipitation_amount_interval = 9999999.9, prho_reference, &
+             pt_reference = 9999999.9, pt_slope_offset = 0.0, &
+             pt_surface = 300.0, pt_surface_initial_change = 0.0, &
+             q_surface = 0.0, q_surface_initial_change = 0.0, &
+             rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, &
+             residual_limit = 1.0E-4, restart_time = 9999999.9, rho_reference, &
+             rho_surface, rif_max = 1.0, &
              rif_min = -5.0, roughness_length = 0.1, sa_surface = 35.0, &
              simulated_time = 0.0, simulated_time_at_begin, sin_alpha_surface, &
@@ -401,7 +404,7 @@
              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, &
-             ug_vertical_gradient_level(10) = -1.0, &
+             ug_vertical_gradient_level(10) = -9999999.9, &
              vg_vertical_gradient(10) = 0.0, &
-             vg_vertical_gradient_level(10) = -1.0, &
+             vg_vertical_gradient_level(10) = -9999999.9, &
              volume_flow(1:2) = 0.0, volume_flow_area(1:2) = 0.0, &
              volume_flow_initial(1:2) = 0.0, wall_heatflux(0:4) = 0.0
Index: palm/trunk/SOURCE/netcdf.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 96)
+++ palm/trunk/SOURCE/netcdf.f90	(revision 97)
@@ -7,5 +7,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Grids defined for rho and sa
 !
 ! Former revisions:
@@ -345,5 +345,5 @@
 !--             Most variables are defined on the scalar grid
                 CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
-                       'ql_vp', 'qv', 's', 'vpt' )
+                       'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -870,5 +870,5 @@
                       CASE ( 'e_xy', 'p_xy', 'pc_xy', 'pr_xy', 'pt_xy', 'q_xy',&
                              'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
-                             'qv_xy', 's_xy', 'vpt_xy' )
+                             'qv_xy', 'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
 
                          grid_x = 'x'
@@ -1414,5 +1414,5 @@
                    CASE ( 'e_xz', 'p_xz', 'pc_xz', 'pr_xz', 'pt_xz', 'q_xz',  &
                           'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qv_xz', &
-                          's_xz', 'vpt_xz' )
+                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
 
                       grid_x = 'x'
@@ -1925,5 +1925,5 @@
                    CASE ( 'e_yz', 'p_yz', 'pc_yz', 'pr_yz', 'pt_yz', 'q_yz',  &
                           'ql_yz', 'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qv_yz', &
-                          's_yz', 'vpt_yz' )
+                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
 
                       grid_x = 'x'
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 96)
+++ palm/trunk/SOURCE/palm.f90	(revision 97)
@@ -63,5 +63,5 @@
     INTEGER           ::  i, run_description_header_i(80)
 
-    version = 'PALM 3.2b'
+    version = 'PALM 3.3'
 
 #if defined( __parallel )
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 96)
+++ palm/trunk/SOURCE/production_e.f90	(revision 97)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! energy production by density flux (in ocean) added
+! use_pt_reference renamed use_reference
 !
 ! Former revisions:
@@ -364,41 +365,76 @@
           IF ( .NOT. humidity )  THEN
 
-             IF ( use_pt_reference )  THEN
-
-                DO  j = nys, nyn
-                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * &
+             IF ( use_reference )  THEN
+
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux in
+!--                the bottom and top surface layer
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i), nzt
+                         tend(k,j,i) = tend(k,j,i) +                    &
+                                       kh(k,j,i) * g / prho_reference * &
+                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   DO  j = nys, nyn
+                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                         tend(k,j,i) = tend(k,j,i) +                  &
+                                       kh(k,j,i) * g / pt_reference * &
                                        ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                      ENDDO
+
+                      IF ( use_surface_fluxes )  THEN
+                         k = nzb_diff_s_inner(j,i)-1
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
+                      ENDIF
+
+                      IF ( use_top_fluxes )  THEN
+                         k = nzt
+                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
+                                                     tswst(j,i)
+                      ENDIF
                    ENDDO
 
-                   IF ( use_surface_fluxes )  THEN
-                      k = nzb_diff_s_inner(j,i)-1
-                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
-                   ENDIF
-
-                   IF ( use_top_fluxes )  THEN
-                      k = nzt
-                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
-                   ENDIF
-                ENDDO
+                ENDIF
 
              ELSE
 
-                DO  j = nys, nyn
-                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux in
+!--                the bottom and top surface layer
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i), nzt
+                         tend(k,j,i) = tend(k,j,i) -                &
+                                       kh(k,j,i) * g / rho(k,j,i) * &
+                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   DO  j = nys, nyn
+                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                         tend(k,j,i) = tend(k,j,i) -               &
+                                       kh(k,j,i) * g / pt(k,j,i) * &
                                        ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                      ENDDO
+
+                      IF ( use_surface_fluxes )  THEN
+                         k = nzb_diff_s_inner(j,i)-1
+                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
+                      ENDIF
+
+                      IF ( use_top_fluxes )  THEN
+                         k = nzt
+                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+                      ENDIF
                    ENDDO
 
-                   IF ( use_surface_fluxes )  THEN
-                      k = nzb_diff_s_inner(j,i)-1
-                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-                   ENDIF
-
-                   IF ( use_top_fluxes )  THEN
-                      k = nzt
-                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-                   ENDIF
-                ENDDO
+                ENDIF
 
              ENDIF
@@ -762,39 +798,65 @@
        IF ( .NOT. humidity )  THEN
 
-          IF ( use_pt_reference )  THEN
-
-             DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
+          IF ( use_reference )  THEN
+
+             IF ( ocean )  THEN
+!
+!--             So far in the ocean no special treatment of density flux in the
+!--             bottom and top surface layer
+                DO  k = nzb_s_inner(j,i), nzt
+                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
+                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
+                ENDDO
+
+             ELSE
+
+                DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
                                        ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-             ENDDO
-
-             IF ( use_surface_fluxes )  THEN
-                k = nzb_diff_s_inner(j,i)-1
-                tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
+                ENDDO
+
+                IF ( use_surface_fluxes )  THEN
+                   k = nzb_diff_s_inner(j,i)-1
+                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
+                ENDIF
+
+                IF ( use_top_fluxes )  THEN
+                   k = nzt
+                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
+                ENDIF
+
              ENDIF
 
-             IF ( use_top_fluxes )  THEN
-                k = nzt
-                tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
+          ELSE
+
+             IF ( ocean )  THEN
+!
+!--             So far in the ocean no special treatment of density flux in the
+!--             bottom and top surface layer
+                DO  k = nzb_s_inner(j,i), nzt
+                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
+                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
+                ENDDO
+
+             ELSE
+
+                DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
+                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                ENDDO
+
+                IF ( use_surface_fluxes )  THEN
+                   k = nzb_diff_s_inner(j,i)-1
+                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
+                ENDIF
+
+                IF ( use_top_fluxes )  THEN
+                   k = nzt
+                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+                ENDIF
+
              ENDIF
 
-          ELSE
-
-             DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
-                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-             ENDDO
-
-             IF ( use_surface_fluxes )  THEN
-                k = nzb_diff_s_inner(j,i)-1
-                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-             ENDIF
-
-             IF ( use_top_fluxes )  THEN
-                k = nzt
-                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-             ENDIF
-
-         ENDIF
+          ENDIF
 
        ELSE
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 96)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 97)
@@ -5,6 +5,9 @@
 ! -----------------
 ! prognostic equation for salinity, density is calculated from equation of
-! state for seawater, +eqn_state_seawater_mod
-! new argument zw in calls of diffusion_e, calc_mean_pt_profile renamed
+! state for seawater and is used for calculation of buoyancy,
+! +eqn_state_seawater_mod
+! diffusion_e is called with argument rho in case of ocean runs,
+! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
+! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
 ! calc_mean_profile
 !
@@ -148,5 +151,5 @@
           ENDIF
           CALL coriolis( i, j, 1 )
-          IF ( sloping_surface )  CALL buoyancy( i, j, pt, 1, 4 )
+          IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
           CALL user_actions( i, j, 'u-tendency' )
 
@@ -283,8 +286,12 @@
           ENDIF
           CALL coriolis( i, j, 3 )
-          IF ( .NOT. humidity )  THEN
-             CALL buoyancy( i, j, pt, 3, 4 )
+          IF ( ocean )  THEN
+             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
           ELSE
-             CALL buoyancy( i, j, vpt, 3, 44 )
+             IF ( .NOT. humidity )  THEN
+                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
+             ELSE
+                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+             ENDIF
           ENDIF
           CALL user_actions( i, j, 'w-tendency' )
@@ -654,9 +661,17 @@
                   .NOT. use_upstream_for_tke )  THEN
                 IF ( .NOT. humidity )  THEN
-                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                     l_grid, pt, rif, tend, zu, zw )
+                   IF ( ocean )  THEN
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,   &
+                                        l_grid, rho, prho_reference, rif, tend, &
+                                        zu, zw )
+                   ELSE
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
+                                        l_grid, pt, pt_reference, rif, tend,  &
+                                        zu, zw )
+                   ENDIF
                 ELSE
-                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                     l_grid, vpt, rif, tend, zu, zw )
+                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,     &
+                                     l_grid, vpt, pt_reference, rif, tend, zu, &
+                                     zw )
                 ENDIF
              ELSE
@@ -680,18 +695,26 @@
                    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, &
-                                        zw )
+                                        km_m, l_grid, pt_m, pt_reference,   &
+                                        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, &
-                                        zw )
+                                        km_m, l_grid, vpt_m, pt_reference,  &
+                                        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, zw )
+                      IF ( ocean )  THEN
+                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
+                                           km, l_grid, rho, prho_reference,   &
+                                           rif, tend, zu, zw )
+                      ELSE
+                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
+                                           km, l_grid, pt, pt_reference, rif, &
+                                           tend, zu, zw )
+                      ENDIF
                    ELSE
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, vpt, rif, tend, zu, zw )
+                                        l_grid, vpt, pt_reference, rif, tend, &
+                                        zu, zw )
                    ENDIF
                 ENDIF
@@ -797,5 +820,6 @@
              ENDIF
              CALL coriolis( i, j, 1 )
-             IF ( sloping_surface )  CALL buoyancy( i, j, pt, 1, 4 )
+             IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, &
+                                                    4 )
              CALL user_actions( i, j, 'u-tendency' )
 
@@ -896,8 +920,12 @@
              ENDIF
              CALL coriolis( i, j, 3 )
-             IF ( .NOT. humidity )  THEN
-                CALL buoyancy( i, j, pt, 3, 4 )
+             IF ( ocean )  THEN
+                CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
              ELSE
-                CALL buoyancy( i, j, vpt, 3, 44 )
+                IF ( .NOT. humidity )  THEN
+                   CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
+                ELSE
+                   CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+                ENDIF
              ENDIF
              CALL user_actions( i, j, 'w-tendency' )
@@ -1114,18 +1142,26 @@
                    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, &
-                                        zw )
+                                        km_m, l_grid, pt_m, pt_reference,   &
+                                        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, &
-                                        zw )
+                                        km_m, l_grid, vpt_m, pt_reference,  &
+                                        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, zw )
+                      IF ( ocean )  THEN
+                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
+                                           km, l_grid, rho, prho_reference,  &
+                                           rif, tend, zu, zw )
+                      ELSE
+                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
+                                           km, l_grid, pt, pt_reference, rif, &
+                                           tend, zu, zw )
+                      ENDIF
                    ELSE
                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, vpt, rif, tend, zu, zw )
+                                        l_grid, vpt, pt_reference, rif, tend, &
+                                        zu, zw )
                    ENDIF
                 ENDIF
@@ -1224,5 +1260,5 @@
     ENDIF
     CALL coriolis( 1 )
-    IF ( sloping_surface )  CALL buoyancy( pt, 1, 4 )
+    IF ( sloping_surface )  CALL buoyancy( pt, pt_reference, 1, 4 )
     CALL user_actions( 'u-tendency' )
 
@@ -1367,8 +1403,12 @@
     ENDIF
     CALL coriolis( 3 )
-    IF ( .NOT. humidity )  THEN
-       CALL buoyancy( pt, 3, 4 )
+    IF ( ocean )  THEN
+       CALL buoyancy( rho, prho_reference, 3, 64 )
     ELSE
-       CALL buoyancy( vpt, 3, 44 )
+       IF ( .NOT. humidity )  THEN
+          CALL buoyancy( pt, pt_reference, 3, 4 )
+       ELSE
+          CALL buoyancy( vpt, pt_reference, 3, 44 )
+       ENDIF
     ENDIF
     CALL user_actions( 'w-tendency' )
@@ -1753,9 +1793,14 @@
        THEN
           IF ( .NOT. humidity )  THEN
-             CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
-                               rif, tend, zu, zw )
+             IF ( ocean )  THEN
+                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
+                                  prho_reference, rif, tend, zu, zw )
+             ELSE
+                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
+                                  pt_reference, rif, tend, zu, zw )
+             ENDIF
           ELSE
              CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
-                               rif, tend, zu, zw )
+                               pt_reference, rif, tend, zu, zw )
           ENDIF
        ELSE
@@ -1777,16 +1822,21 @@
              IF ( .NOT. humidity )  THEN
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
-                                  pt_m, rif_m, tend, zu, zw )
+                                  pt_m, pt_reference, 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, zw )
+                                  vpt_m, pt_reference, 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, zw )
+                IF ( ocean )  THEN
+                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
+                                     rho, prho_reference, rif, tend, zu, zw )
+                ELSE
+                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
+                                     pt, pt_reference, rif, tend, zu, zw )
+                ENDIF
              ELSE
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
-                                  rif, tend, zu, zw )
+                                  pt_reference, rif, tend, zu, zw )
              ENDIF
           ENDIF
Index: palm/trunk/SOURCE/run_control.f90
===================================================================
--- palm/trunk/SOURCE/run_control.f90	(revision 96)
+++ palm/trunk/SOURCE/run_control.f90	(revision 97)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Timestep format changed
 !
 ! Former revisions:
@@ -113,15 +113,15 @@
 100 FORMAT (///'Run-control output:'/ &
               &'------------------'// &
-           &'RUN  ITER. HH:MM:SS.SS   DT(E)     UMAX     VMAX     WMAX     U*', &
-           &'    W*   THETA*   Z_I     ENERG.   DISTENERG    DIVOLD     DIVNE', &
-           &'W     UMAX(KJI)    VMAX(KJI)    WMAX(KJI)   ADVECX   ADVECY   MG', &
-           &'CYC'/ &
-           &'----------------------------------------------------------------', &
-           &'----------------------------------------------------------------', &
-           &'----------------------------------------------------------------', &
-           &'--')
-101 FORMAT (I3,1X,I6,1X,A8,F3.2,1X,F7.4,A1,A1,F8.4,A1,F8.4,A1,F8.4,2X,F5.3,2X, &
+          &'RUN  ITER. HH:MM:SS.SS    DT(E)     UMAX     VMAX     WMAX     U', &
+          &'*    W*   THETA*    Z_I     ENERG.   DISTENERG    DIVOLD     DIV', &
+          &'NEW     UMAX(KJI)    VMAX(KJI)    WMAX(KJI)   ADVECX   ADVECY   ', &
+          &'MGCYC'/                                                            &
+          &'----------------------------------------------------------------', &
+          &'----------------------------------------------------------------', &
+          &'----------------------------------------------------------------', &
+          &'-----')
+101 FORMAT (I3,1X,I6,1X,A8,F3.2,1X,F8.4,A1,A1,F8.4,A1,F8.4,A1,F8.4,2X,F5.3,2X, &
             F4.2, &
-            2X,F6.3,2X,F5.0,1X,4(E10.3,1X),3(3(I4),1X),F8.3,1X,F8.3,5X,I3)
+            2X,F6.3,2X,F6.0,1X,4(E10.3,1X),3(3(I4),1X),F8.3,1X,F8.3,5X,I3)
 
  END SUBROUTINE run_control
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 96)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 97)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! Ghostpoint exchange for salinity and density
+! diffusivities is called with argument rho in case of ocean runs,
+! new argument pt_/prho_reference in calls of diffusivities,
+! ghostpoint exchange for salinity and density
 !
 ! Former revisions:
@@ -249,7 +251,11 @@
              CALL cpu_log( log_point(17), 'diffusivities', 'start' )
              IF ( .NOT. humidity ) THEN
-                CALL diffusivities( pt )
+                IF ( ocean )  THEN
+                   CALL diffusivities( rho, prho_reference )
+                ELSE
+                   CALL diffusivities( pt, pt_reference )
+                ENDIF
              ELSE
-                CALL diffusivities( vpt )
+                CALL diffusivities( vpt, pt_reference )
              ENDIF
              CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
