Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 95)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 96)
@@ -4,5 +4,5 @@
 + inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux
 
-boundary_conds, check_parameters, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, parin, prognostic_equations, read_var_list, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary
+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:
@@ -11,7 +11,7 @@
 Changed:
 -------
-hydro_press renamed hyp
+hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile
 
-advec_particles
+advec_particles, buoyancy, calc_liquid_water_content, init_cloud_physics, modules, prognostic_equations
 
 
Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 95)
+++ /palm/trunk/SOURCE/Makefile	(revision 96)
@@ -209,5 +209,6 @@
 	buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \
 	diffusion_e.o diffusion_s.o diffusion_u.o diffusion_v.o diffusion_w.o \
-	impact_of_latent_heat.o production_e.o user_interface.o
+	eqn_state_seawater.o impact_of_latent_heat.o production_e.o \
+	user_interface.o
 random_gauss.o: random_function.o
 read_3d_binary.o: modules.o random_function.o
Index: /palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/average_3d_data.f90	(revision 95)
+++ /palm/trunk/SOURCE/average_3d_data.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Averaging of density and salinity
 !
 ! Former revisions:
@@ -166,4 +166,13 @@
              ENDDO
              
+          CASE ( 'rho' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d )
+                   ENDDO
+                ENDDO
+             ENDDO
+             
           CASE ( 's' )
              DO  i = nxl-1, nxr+1
@@ -171,4 +180,13 @@
                    DO  k = nzb, nzt+1
                       s_av(k,j,i) = s_av(k,j,i) / REAL( average_count_3d )
+                   ENDDO
+                ENDDO
+             ENDDO
+             
+          CASE ( 'sa' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d )
                    ENDDO
                 ENDDO
Index: /palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- /palm/trunk/SOURCE/buoyancy.f90	(revision 95)
+++ /palm/trunk/SOURCE/buoyancy.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! calc_mean_pt_profile renamed calc_mean_profile
 !
 ! Former revisions:
@@ -29,5 +29,5 @@
 
     PRIVATE
-    PUBLIC buoyancy, calc_mean_pt_profile
+    PUBLIC buoyancy, calc_mean_profile
 
     INTERFACE buoyancy
@@ -36,7 +36,7 @@
     END INTERFACE buoyancy
 
-    INTERFACE calc_mean_pt_profile
-       MODULE PROCEDURE calc_mean_pt_profile
-    END INTERFACE calc_mean_pt_profile
+    INTERFACE calc_mean_profile
+       MODULE PROCEDURE calc_mean_profile
+    END INTERFACE calc_mean_profile
 
  CONTAINS
@@ -207,5 +207,5 @@
 
 
-    SUBROUTINE calc_mean_pt_profile( theta, pr )
+    SUBROUTINE calc_mean_profile( var, pr )
 
 !------------------------------------------------------------------------------!
@@ -224,19 +224,17 @@
 
        INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
-       REAL, DIMENSION(:,:,:), POINTER  ::  theta
-
-!
-!--    Computation of the horizontally averaged temperature profile, unless
+       REAL, DIMENSION(:,:,:), POINTER  ::  var
+
+!
+!--    Computation of the horizontally averaged profile of variable var, unless
 !--    already done by the relevant call from flow_statistics. The calculation
 !--    is done only for the first respective intermediate timestep in order to
 !--    spare communication time and to produce identical model results with jobs
 !--    which are calling flow_statistics at different time intervals.
-!--    Although this calculation is not required for model runs with a slope,
-!--    it is nevertheless also computed.
        IF ( .NOT. flow_statistics_called  .AND. &
             intermediate_timestep_count == 1 )  THEN
 
 !
-!--       Horizontal average of the potential temperature
+!--       Horizontal average of variable var
           tn           =   0  ! Default thread number in case of one thread
           !$OMP PARALLEL PRIVATE( i, j, k, tn )
@@ -247,5 +245,5 @@
              DO  j =  nys, nyn
                 DO  k = nzb_s_outer(j,i), nzt+1
-                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + theta(k,j,i)
+                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
                 ENDDO
              ENDDO
@@ -272,5 +270,5 @@
        ENDIF
 
-    END SUBROUTINE calc_mean_pt_profile
+    END SUBROUTINE calc_mean_profile
 
  END MODULE buoyancy_mod
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 95)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 96)
@@ -1485,4 +1485,24 @@
              hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
 
+          CASE ( 'sa', '#sa' )
+             IF ( .NOT. ocean )  THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for ocean = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 23
+                dopr_unit(i)  = 'psu'
+                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 26
+                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
           CASE ( 'u*2' )
              dopr_index(i) = 30
@@ -1846,4 +1866,51 @@
              dopr_unit(i)  = 'K m/s'
              hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'rho' )
+             dopr_index(i) = 64
+             dopr_unit(i)  = 'kg/m3'
+             hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'w"sa"' )
+             IF ( .NOT. ocean ) THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for ocean = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 65
+                dopr_unit(i)  = 'psu m/s'
+                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'w*sa*' )
+             IF ( .NOT. ocean ) THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for ocean = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 66
+                dopr_unit(i)  = 'psu m/s'
+                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'wsa' )
+             IF ( .NOT. ocean ) THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for ocean = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 67
+                dopr_unit(i)  = 'psu m/s'
+                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
 
 
@@ -2068,4 +2135,14 @@
              unit = 'kg/kg'
 
+          CASE ( 'rho' )
+             IF ( .NOT. ocean )  THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
+                                '" requires ocean = .TRUE.'
+                ENDIF
+                CALL local_stop
+             ENDIF
+             unit = 'kg/m3'
+
           CASE ( 's' )
              IF ( .NOT. passive_scalar )  THEN
@@ -2077,4 +2154,14 @@
              ENDIF
              unit = 'conc'
+
+          CASE ( 'sa' )
+             IF ( .NOT. ocean )  THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
+                                '" requires ocean = .TRUE.'
+                ENDIF
+                CALL local_stop
+             ENDIF
+             unit = 'psu'
 
           CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 95)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Output of density and salinity
 !
 ! Former revisions:
@@ -370,4 +370,11 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rho
+                ELSE
+                   to_be_resorted => rho_av
+                ENDIF
+
              CASE ( 's_xy', 's_xz', 's_yz' )
                 IF ( av == 0 )  THEN
@@ -375,4 +382,11 @@
                 ELSE
                    to_be_resorted => q_av
+                ENDIF
+
+             CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => sa
+                ELSE
+                   to_be_resorted => sa_av
                 ENDIF
 
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 95)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Output of density and salinity
 !
 ! Former revisions:
@@ -251,4 +251,11 @@
              ENDIF
 
+          CASE ( 'rho' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rho
+             ELSE
+                to_be_resorted => rho_av
+             ENDIF
+             
           CASE ( 's' )
              IF ( av == 0 )  THEN
@@ -256,4 +263,11 @@
              ELSE
                 to_be_resorted => q_av
+             ENDIF
+             
+          CASE ( 'sa' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => sa
+             ELSE
+                to_be_resorted => sa_av
              ENDIF
              
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 96)
@@ -0,0 +1,198 @@
+ MODULE eqn_state_seawater_mod
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id:$
+!
+! Initial revision
+!
+!
+! Description:
+! ------------
+! Equation of state for seawater as a function of potential temperature,
+! salinity, and pressure.
+! For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
+!------------------------------------------------------------------------------!
+
+    IMPLICIT NONE
+
+    PRIVATE
+    PUBLIC eqn_state_seawater, eqn_state_seawater_func
+
+    REAL, DIMENSION(12), PARAMETER ::  nom =                             &
+                      (/ 9.9984085444849347D2,   7.3471625860981584D0,   &
+                        -5.3211231792841769D-2,  3.6492439109814549D-4,  &
+                         2.5880571023991390D0,  -6.7168282786692354D-3,  &
+                         1.9203202055760151D-3,  1.1798263740430364D-2,  &
+                         9.8920219266399117D-8,  4.6996642771754730D-6,  &
+                        -2.5862187075154352D-8, -3.2921414007960662D-12 /)
+
+    REAL, DIMENSION(13), PARAMETER ::  den =                             &
+                      (/ 1.0,                    7.2815210113327091D-3,  &
+                        -4.4787265461983921D-5,  3.3851002965802430D-7,  &
+                         1.3651202389758572D-10, 1.7632126669040377D-3,  &
+                        -8.8066583251206474D-6, -1.8832689434804897D-10, &
+                         5.7463776745432097D-6,  1.4716275472242334D-9,  &
+                         6.7103246285651894D-6, -2.4461698007024582D-17, &
+                        -9.1534417604289062D-18 /)
+
+    INTERFACE eqn_state_seawater
+       MODULE PROCEDURE eqn_state_seawater
+       MODULE PROCEDURE eqn_state_seawater_ij
+    END INTERFACE eqn_state_seawater
+ 
+    INTERFACE eqn_state_seawater_func
+       MODULE PROCEDURE eqn_state_seawater_func
+    END INTERFACE eqn_state_seawater_func
+ 
+ CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE eqn_state_seawater
+
+       USE arrays_3d
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+
+       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_u_inner(j,i)+1, nzt
+!
+!--             Pressure is needed in dbar
+                p1 = hyp(k) * 1E-4
+                p2 = p1 * p1
+                p3 = p2 * p1
+
+!
+!--             Temperature needed in degree Celsius
+                pt1 = pt_p(k,j,i) - 273.15
+                pt2 = pt1 * pt1
+                pt3 = pt1 * pt2
+                pt4 = pt2 * pt2
+
+                sa1  = sa_p(k,j,i)
+                sa15 = sa1 * SQRT( sa1 )
+                sa2  = sa1 * sa1
+
+                rho(k,j,i) =                                                   &
+                        ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
+                          nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
+                          nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
+                          nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
+                        ) /                                                    &
+                        ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
+                          den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
+                          den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
+                          den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
+                          den(13)*p3*pt1                                       &
+                        )
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE eqn_state_seawater
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE eqn_state_seawater_ij( i, j )
+
+       USE arrays_3d
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+
+       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+
+       DO  k = nzb_u_inner(j,i)+1, nzt
+!
+!--       Pressure is needed in dbar
+          p1 = hyp(k) * 1E-4
+          p2 = p1 * p1
+          p3 = p2 * p1
+
+!
+!--       Temperature needed in degree Celsius
+          pt1 = pt_p(k,j,i) - 273.15
+          pt2 = pt1 * pt1
+          pt3 = pt1 * pt2
+          pt4 = pt2 * pt2
+
+          sa1  = sa_p(k,j,i)
+          sa15 = sa1 * SQRT( sa1 )
+          sa2  = sa1 * sa1
+
+          rho(k,j,i) = ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
+                         nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
+                         nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
+                         nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
+                       ) /                                                    &
+                       ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
+                         den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
+                         den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
+                         den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
+                         den(13)*p3*pt1                                       &
+                       )
+       ENDDO
+
+    END SUBROUTINE eqn_state_seawater_ij
+
+
+!------------------------------------------------------------------------------!
+! Equation of state as a function
+!------------------------------------------------------------------------------!
+    REAL FUNCTION eqn_state_seawater_func( p, pt, sa )
+
+       IMPLICIT NONE
+
+       REAL ::  p, p1, p2, p3, pt, pt1, pt2, pt3, pt4, sa, sa15, sa2
+
+!
+!--    Pressure is needed in dbar
+       p1 = p  * 1E-4
+       p2 = p1 * p1
+       p3 = p2 * p1
+
+!
+!--    Temperature needed in degree Celsius
+       pt1 = pt - 273.15
+       pt2 = pt1 * pt1
+       pt3 = pt1 * pt2
+       pt4 = pt2 * pt2
+
+       sa15 = sa * SQRT( sa )
+       sa2  = sa * sa
+
+
+       eqn_state_seawater_func =                                               &
+         ( nom(1)        + nom(2)*pt1       + nom(3)*pt2    + nom(4)*pt3     + &
+           nom(5)*sa     + nom(6)*sa*pt1    + nom(7)*sa2    + nom(8)*p1      + &
+           nom(9)*p1*pt2 + nom(10)*p1*sa    + nom(11)*p2    + nom(12)*p2*pt2   &
+         ) /                                                                   &
+         ( den(1)        + den(2)*pt1       + den(3)*pt2    + den(4)*pt3     + &
+           den(5)*pt4    + den(6)*sa        + den(7)*sa*pt1 + den(8)*sa*pt3  + &
+           den(9)*sa15   + den(10)*sa15*pt2 + den(11)*p1    + den(12)*p2*pt3 + &
+           den(13)*p3*pt1                                                      &
+         )
+
+
+    END FUNCTION eqn_state_seawater_func
+
+ END MODULE eqn_state_seawater_mod
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 95)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 96)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Calculation of Deardorff velocity scale adjusted to be used with the ocean
+! 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!!!)
 !
@@ -123,4 +124,18 @@
           ENDDO
        ENDDO
+
+!
+!--    Horizontally averaged profile of salinity
+       IF ( ocean )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                DO  k = nzb_s_outer(j,i), nzt+1
+                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
+                                      sa(k,j,i) * rmask(j,i,sr)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
 
 !
@@ -177,4 +192,7 @@
              sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
              sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
+             IF ( ocean )  THEN
+                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
+             ENDIF
              IF ( humidity )  THEN
                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
@@ -200,4 +218,8 @@
        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
                            MPI_SUM, comm2d, ierr )
+       IF ( ocean )  THEN
+          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+       ENDIF
        IF ( humidity ) THEN
           CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
@@ -221,4 +243,5 @@
        sums(:,2) = sums_l(:,2,0)
        sums(:,4) = sums_l(:,4,0)
+       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
        IF ( humidity ) THEN
           sums(:,44) = sums_l(:,44,0)
@@ -241,4 +264,11 @@
        hom(:,1,2,sr) = sums(:,2)             ! v
        hom(:,1,4,sr) = sums(:,4)             ! pt
+
+!
+!--    Salinity
+       IF ( ocean )  THEN
+          sums(:,23) = sums(:,23) / ngp_2dh_outer(:,sr)
+          hom(:,1,23,sr) = sums(:,23)             ! sa
+       ENDIF
 
 !
@@ -377,4 +407,13 @@
 
 !
+!--             Salinity flux w"sa"
+                IF ( ocean )  THEN
+                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
+                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
+                                               * ddzu(k+1) * rmask(j,i,sr)
+                ENDIF
+
+!
 !--             Buoyancy flux, water flux (humidity flux) w"q"
                 IF ( humidity ) THEN
@@ -420,4 +459,8 @@
                 sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                     0.0 * rmask(j,i,sr)           ! v"pt"
+                IF ( ocean )  THEN
+                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
+                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
+                ENDIF
                 IF ( humidity )  THEN
                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
@@ -450,4 +493,8 @@
                 sums_l(nzt,61,tn) = sums_l(nzt,61,tn) + &
                                     0.0 * rmask(j,i,sr)           ! v"pt"
+                IF ( ocean )  THEN
+                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
+                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
+                ENDIF
                 IF ( humidity )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
@@ -508,4 +555,16 @@
                 sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
                                                     rmask(j,i,sr)
+
+!
+!--             Salinity flux and density (density does not belong to here,
+!--             but so far there is no suitable place to calculate)
+                IF ( ocean )  THEN
+                   pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
+                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
+                   sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
+                                                       rmask(j,i,sr)
+                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
+                                                       rmask(j,i,sr)
+                ENDIF
 
 !
@@ -739,5 +798,6 @@
        hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
        hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
-                                       ! profiles 23-29 left empty for initial 
+                                       ! profile 24 is initial profile (sa)
+                                       ! profiles 25-29 left empty for initial 
                                        ! profiles
        hom(:,1,30,sr) = sums(:,30)     ! u*2
@@ -771,4 +831,8 @@
        hom(:,1,62,sr) = sums(:,62)     ! v*pt*
        hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
+       hom(:,1,64,sr) = sums(:,64)     ! rho
+       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
+       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
+       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 95)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 96)
@@ -190,8 +190,10 @@
        ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
                  saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
-       ALLOCATE( rho(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), &
+       ALLOCATE( rho_1(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) )
+       rho => rho_1  ! routine calc_mean_profile requires density to be a
+                     ! pointer
     ENDIF
 
@@ -826,17 +828,19 @@
 
 !
-!-- If required, initialize quantities needed for the ocean model
-    IF ( ocean )  CALL init_ocean
-
-!
 !-- If required, initialize dvrp-software
     IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
 
-!
-!-- If required, initialize quantities for handling cloud physics
-!-- This routine must be called before init_particles, because
-!-- otherwise, array pt_d_t, needed in data_output_dvrp (called by
-!-- init_particles) is not defined.
-    CALL init_cloud_physics
+    IF ( ocean )  THEN
+!
+!--    Initialize quantities needed for the ocean model
+       CALL init_ocean
+    ELSE
+!
+!--    Initialize quantities for handling cloud physics
+!--    This routine must be called before init_particles, because
+!--    otherwise, array pt_d_t, needed in data_output_dvrp (called by
+!--    init_particles) is not defined.
+       CALL init_cloud_physics
+    ENDIF
 
 !
Index: /palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- /palm/trunk/SOURCE/init_ocean.f90	(revision 95)
+++ /palm/trunk/SOURCE/init_ocean.f90	(revision 96)
@@ -8,5 +8,5 @@
 ! Former revisions:
 ! ------------------
-! $Id:$
+! $Id$
 !
 ! Initial revision (raasch 01/06/07)
@@ -37,7 +37,9 @@
 !
 !-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
+!-- and the reference density (used later in buoyancy term)
     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)
 
     DO  k = nzt-1, 0, -1
@@ -46,10 +48,14 @@
        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)
+       hyp(k)  = hyp(k+1) + rho_l * g * dzu(k+1)
+       rho_ref = rho_ref + rho_l * dzu(k+1)
 
     ENDDO
 
+    rho_ref = rho_ref / ( zw(nzt) - zu(nzb) )
+    print*, '*** rho_ref = ', rho_ref
+
 
  END SUBROUTINE init_ocean
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 95)
+++ /palm/trunk/SOURCE/modules.f90	(revision 96)
@@ -119,11 +119,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, rho, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1,     &
+          q_2, q_3, ql_1, ql_2, rho_1, 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, 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
+          ql_c, rho, 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
@@ -149,5 +149,5 @@
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
           e_av, p_av, pc_av, pr_av, pt_av, q_av, ql_av, ql_c_av, ql_v_av, &
-          ql_vp_av, qv_av, s_av, u_av, v_av, vpt_av, w_av
+          ql_vp_av, qv_av, rho_av, s_av, sa_av, u_av, v_av, vpt_av, w_av
 
  END MODULE averaging
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 95)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 96)
@@ -4,6 +4,8 @@
 ! Actual revisions:
 ! -----------------
-! prognostic equation for salinity
-! new argument zw in calls of diffusion_e
+! 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
+! calc_mean_profile
 !
 ! Former revisions:
@@ -41,4 +43,5 @@
     USE control_parameters
     USE cpulog
+    USE eqn_state_seawater_mod
     USE grid_variables
     USE indices
@@ -107,6 +110,7 @@
 !-- Calculate those variables needed in the tendency terms which need
 !-- global communication
-    CALL calc_mean_pt_profile( pt, 4 )
-    IF ( humidity )  CALL calc_mean_pt_profile( vpt, 44 )
+    CALL calc_mean_profile( pt, 4 )
+    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
 
 !
@@ -495,4 +499,8 @@
                 ENDIF
              ENDIF
+
+!
+!--          Calculate density by the equation of state for seawater
+             CALL eqn_state_seawater( i, j )
 
           ENDDO
@@ -758,6 +766,7 @@
 !-- Calculate those variables needed in the tendency terms which need
 !-- global communication
-    CALL calc_mean_pt_profile( pt, 4 )
-    IF ( humidity )  CALL calc_mean_pt_profile( vpt, 44 )
+    CALL calc_mean_profile( pt, 4 )
+    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
     IF ( .NOT. constant_diffusion )  CALL production_e_init
 
@@ -1020,4 +1029,8 @@
                    ENDIF
                 ENDIF
+
+!
+!--             Calculate density by the equation of state for seawater
+                CALL eqn_state_seawater( i, j )
 
              ENDIF
@@ -1178,6 +1191,7 @@
 !-- Calculate those variables needed in the tendency terms which need
 !-- global communication
-    CALL calc_mean_pt_profile( pt, 4 )
-    IF ( humidity )  CALL calc_mean_pt_profile( vpt, 44 )
+    CALL calc_mean_profile( pt, 4 )
+    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
 
 !
@@ -1589,4 +1603,10 @@
 
        CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
+
+!
+!--    Calculate density by the equation of state for seawater
+       CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
+       CALL eqn_state_seawater
+       CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
 
     ENDIF
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 95)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +sa, saswsb, saswst
+! +rho_av, sa, sa_av, saswsb, saswst
 !
 ! Former revisions:
@@ -273,4 +273,7 @@
              READ ( 13 )  random_iv
              READ ( 13 )  random_iy
+          CASE ( 'rho_av' )
+             ALLOCATE( rho_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+             READ ( 13 )  rho_av
           CASE ( 'rif' )
              READ ( 13 )  rif
@@ -284,4 +287,7 @@
           CASE ( 'sa' )
              READ ( 13 )  sa
+          CASE ( 'sa_av' )
+             ALLOCATE( sa_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+             READ ( 13 )  sa_av
           CASE ( 'saswsb' )
              READ ( 13 )  saswsb
Index: /palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 95)
+++ /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +sum-up of density and salinity
 !
 ! Former revisions:
@@ -132,4 +132,10 @@
                 qv_av = 0.0
 
+             CASE ( 'rho' )
+                IF ( .NOT. ALLOCATED( rho_av ) )  THEN
+                   ALLOCATE( rho_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                rho_av = 0.0
+
              CASE ( 's' )
                 IF ( .NOT. ALLOCATED( s_av ) )  THEN
@@ -137,4 +143,10 @@
                 ENDIF
                 s_av = 0.0
+
+             CASE ( 'sa' )
+                IF ( .NOT. ALLOCATED( sa_av ) )  THEN
+                   ALLOCATE( sa_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+                ENDIF
+                sa_av = 0.0
 
              CASE ( 't*' )
@@ -337,4 +349,13 @@
              ENDDO
 
+          CASE ( 'rho' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      rho_av(k,j,i) = rho_av(k,j,i) + rho(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+             
           CASE ( 's' )
              DO  i = nxl-1, nxr+1
@@ -342,4 +363,13 @@
                    DO  k = nzb, nzt+1
                       s_av(k,j,i) = s_av(k,j,i) + q(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+             
+          CASE ( 'sa' )
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   DO  k = nzb, nzt+1
+                      sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
                    ENDDO
                 ENDDO
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 95)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 96)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +sa, saswsb, saswst
+! +rho_av, sa, sa_av, saswsb, saswst
 !
 ! Former revisions:
@@ -133,5 +133,11 @@
     ENDIF
     IF ( ocean )  THEN
+       IF ( ALLOCATED( rho_av ) )  THEN
+          WRITE ( 14 )  'rho_av              ';  WRITE ( 14 )  rho_av
+       ENDIF
        WRITE ( 14 )  'sa                  ';  WRITE ( 14 )  sa
+       IF ( ALLOCATED( sa_av ) )  THEN
+          WRITE ( 14 )  'sa_av               ';  WRITE ( 14 )  sa_av
+       ENDIF
        WRITE ( 14 )  'saswsb              ';  WRITE ( 14 )  saswsb
        WRITE ( 14 )  'saswst              ';  WRITE ( 14 )  saswst
