Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 94)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 95)
@@ -2,14 +2,16 @@
 ---
 ocean version including prognostic equation for salinity and equation of state for seawater
-+ inipar-parameters ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinity_flux
++ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux
 
-check_parameters, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, parin, prognostic_equations, read_var_list, write_var_list, write_3d_binary
+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
 
+New:
+eqn_state_seawater, init_ocean
 
 Changed:
 -------
+hydro_press renamed hyp
 
-
-
+advec_particles
 
 
@@ -36,2 +38,6 @@
 
 default disturbance level
+
+Bott-Chlond scheme to be extended for salinity
+
+Upstream-spline scheme to be extended for salinity?????
Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 94)
+++ /palm/trunk/SOURCE/Makefile	(revision 95)
@@ -4,9 +4,12 @@
 # Actual revisions:
 # -----------------
-# +local_flush
+# +eqn_state_seawater, init_ocean
 #
 # Former revisions:
 # -----------------
 # $Id$
+#
+# 82 2007-04-16 15:40:52Z raasch
+# +local_flush
 #
 # 58 2007-03-09 14:27:38Z raasch
@@ -45,9 +48,9 @@
 	data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 \
 	diffusion_v.f90 diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
-	disturb_heatflux.f90 \exchange_horiz.f90 exchange_horiz_2d.f90 \
+	disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 exchange_horiz_2d.f90 \
 	fft_xy.f90 flow_statistics.f90 global_min_max.f90 header.f90 \
 	impact_of_latent_heat.f90 init_1d_model.f90 init_3d_model.f90 \
 	init_advec.f90 init_cloud_physics.f90 init_dvrp.f90 init_grid.f90 \
-	init_particles.f90 init_pegrid.f90 init_pt_anomaly.f90 \
+	init_ocean.f90 init_particles.f90 init_pegrid.f90 init_pt_anomaly.f90 \
 	init_rankine.f90 init_slope.f90 interaction_droplets_ptq.f90 \
 	local_flush.f90 local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
@@ -76,8 +79,8 @@
 	data_output_3d.o diffusion_e.o diffusion_s.o diffusion_u.o \
 	diffusion_v.o diffusion_w.o diffusivities.o disturb_field.o \
-	disturb_heatflux.o exchange_horiz.o exchange_horiz_2d.o fft_xy.o \
+	disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o exchange_horiz_2d.o fft_xy.o \
 	flow_statistics.o global_min_max.o header.o impact_of_latent_heat.o \
 	init_1d_model.o init_3d_model.o init_advec.o init_cloud_physics.o \
-	init_dvrp.o init_grid.o init_particles.o init_pegrid.o \
+	init_dvrp.o init_grid.o init_ocean.o init_particles.o init_pegrid.o \
 	init_pt_anomaly.o init_rankine.o init_slope.o \
 	interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \
@@ -164,4 +167,5 @@
 disturb_field.o: modules.o random_function.o
 disturb_heatflux.o: modules.o
+eqn_state_seawater.o: modules.o
 exchange_horiz.o: modules.o
 exchange_horiz_2d.o: modules.o
@@ -177,4 +181,5 @@
 init_dvrp.o: modules.o
 init_grid.o: modules.o
+init_ocean.o: modules.o eqn_state_seawater.o
 init_particles.o: modules.o random_function.o
 init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
Index: /palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- /palm/trunk/SOURCE/advec_particles.f90	(revision 94)
+++ /palm/trunk/SOURCE/advec_particles.f90	(revision 95)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! hydro_press renamed hyp
 ! TEST: PRINT statements on unit 9 (commented out)
 !
@@ -238,6 +239,5 @@
 !
 !--       Calculate real temperature and saturation vapor pressure
-          p_int = hydro_press(k) + ( particles(n)%z - zu(k) ) / dz * &
-                                   ( hydro_press(k+1) - hydro_press(k) )
+          p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) )
           t_int = pt_int * ( p_int / 100000.0 )**0.286
 
Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 94)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Boundary conditions for salinity added
 !
 ! Former revisions:
@@ -128,6 +128,28 @@
 
 !
+!--    Boundary conditions for salinity
+       IF ( ocean )  THEN
+!
+!--       Bottom boundary: Neumann condition because salinity flux is always
+!--       given
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
+             ENDDO
+          ENDDO
+
+!
+!--       Top boundary: Dirichlet or Neumann
+          IF ( ibc_sa_t == 0 )  THEN
+             sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
+          ELSEIF ( ibc_sa_t == 1 )  THEN
+             sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
+          ENDIF
+
+       ENDIF
+
+!
 !--    Boundary conditions for total water content or scalar,
-!--    bottom and surface boundary (see also temperature)
+!--    bottom and top boundary (see also temperature)
        IF ( humidity  .OR.  passive_scalar )  THEN
 !
Index: /palm/trunk/SOURCE/calc_liquid_water_content.f90
===================================================================
--- /palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 94)
+++ /palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! hydro_press renamed hyp
 !
 ! Former revisions:
@@ -59,6 +59,5 @@
 !
 !--          Compute approximation of saturation humidity
-             q_s = 0.622 * e_s / &
-                   ( hydro_press(k) - 0.378 * e_s )
+             q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
 
 !
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 94)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 95)
@@ -1005,4 +1005,38 @@
 
 !
+!-- Boundary conditions for salinity
+    IF ( ocean )  THEN
+       IF ( bc_sa_t == 'dirichlet' )  THEN
+          ibc_sa_t = 0
+       ELSEIF ( bc_sa_t == 'neumann' )  THEN
+          ibc_sa_t = 1
+       ELSE
+          IF ( myid == 0 )  THEN
+             PRINT*, '+++ check_parameters:'
+             PRINT*, '    unknown boundary condition: bc_sa_t = ', bc_sa_t
+          ENDIF
+          CALL local_stop
+       ENDIF
+
+       IF ( top_salinityflux == 9999999.9 )  constant_top_salinityflux = .FALSE.
+
+!
+!--    A fixed salinity at the top implies Dirichlet boundary condition for
+!--    salinity. In this case specification of a constant salinity flux is
+!--    forbidden.
+       IF ( ibc_sa_t == 0  .AND.   constant_top_salinityflux  .AND. &
+            top_salinityflux /= 0.0 )  THEN
+          IF ( myid == 0 )  THEN
+             PRINT*, '+++ check_parameters:'
+             PRINT*, '    boundary_condition: bc_sa_t = ', bc_sa_t
+             PRINT*, '    is not allowed with constant_top_salinityflux = ', &
+                          '.TRUE.'
+          ENDIF
+          CALL local_stop
+       ENDIF
+
+    ENDIF
+
+!
 !-- In case of humidity or passive scalar, set boundary conditions for total 
 !-- water content / scalar
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 94)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 95)
@@ -7,5 +7,5 @@
 ! Actual revisions:
 ! -----------------
-! Initialization of salinity
+! Initialization of salinity, call of init_ocean
 !
 ! Former revisions:
@@ -188,6 +188,7 @@
 
     IF ( ocean )  THEN
-       ALLOCATE( saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
-       ALLOCATE( r(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
+       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), &
@@ -303,4 +304,5 @@
 
        IF ( ocean )  THEN
+          saswsb => saswsb_1
           saswst => saswst_1
           sa     => sa_1;    sa_p => sa_2;    tsa_m => sa_3
@@ -602,4 +604,5 @@
 
              IF ( ocean )  THEN
+                saswsb = bottom_salinityflux
                 saswst = top_salinityflux
              ENDIF
@@ -821,4 +824,8 @@
        CALL local_stop
     ENDIF
+
+!
+!-- If required, initialize quantities needed for the ocean model
+    IF ( ocean )  CALL init_ocean
 
 !
Index: /palm/trunk/SOURCE/init_cloud_physics.f90
===================================================================
--- /palm/trunk/SOURCE/init_cloud_physics.f90	(revision 94)
+++ /palm/trunk/SOURCE/init_cloud_physics.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! hydro_press renamed hyp
 !
 ! Former revisions:
@@ -35,5 +35,5 @@
     REAL    ::  t_surface
 
-    ALLOCATE( hydro_press(nzb:nzt+1), pt_d_t(nzb:nzt+1), t_d_pt(nzb:nzt+1) )
+    ALLOCATE( hyp(nzb:nzt+1), pt_d_t(nzb:nzt+1), t_d_pt(nzb:nzt+1) )
 
 !
@@ -54,11 +54,11 @@
 !-- pt / t : ratio of potential and actual temperature (pt_d_t)
 !-- t / pt : ratio of actual and potential temperature (t_d_pt)
-!-- p_0(z) : vertical profile of the hydrostatic pressure (hydro_press)
+!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
     t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286
     DO  k = nzb, nzt+1
-       hydro_press(k) = surface_pressure * 100.0 * &
-                        ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
-       pt_d_t(k)      = ( 100000.0 / hydro_press(k) )**0.286
-       t_d_pt(k)      = 1.0 / pt_d_t(k)       
+       hyp(k)    = surface_pressure * 100.0 * &
+                   ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
+       pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
+       t_d_pt(k) = 1.0 / pt_d_t(k)       
     ENDDO
 
Index: /palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- /palm/trunk/SOURCE/init_ocean.f90	(revision 95)
+++ /palm/trunk/SOURCE/init_ocean.f90	(revision 95)
@@ -0,0 +1,55 @@
+ SUBROUTINE init_ocean
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! hydro_press renamed hyp
+!
+! Former revisions:
+! ------------------
+! $Id:$
+!
+! Initial revision (raasch 01/06/07)
+!
+! Description:
+! ------------
+! Initialization of quantities needed for the ocean version
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE control_parameters
+    USE eqn_state_seawater_mod
+    USE grid_variables
+    USE indices
+
+    IMPLICIT NONE
+
+    INTEGER ::  k
+
+    REAL    ::  sa_l, pt_l, rho_l
+
+    ALLOCATE( hyp(nzb:nzt+1) )
+
+!
+!-- Set water density near the ocean surface
+    rho_surface = 1027.62
+
+!
+!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
+    hyp(nzt+1) = surface_pressure * 100.0
+
+    hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
+
+    DO  k = nzt-1, 0, -1
+
+       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
+       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
+
+       rho_l = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
+
+       hyp(k) = hyp(k+1) + rho_l * g * dzu(k+1)
+
+    ENDDO
+
+
+ END SUBROUTINE init_ocean
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 94)
+++ /palm/trunk/SOURCE/modules.f90	(revision 95)
@@ -7,4 +7,5 @@
 ! +ocean, r, + salinity variables
 ! defaults of .._vertical_gradient_levels changed from -1.0 to -9999999.9
+! hydro_press renamed hyp
 !
 ! Former revisions:
@@ -94,7 +95,7 @@
 
     REAL, DIMENSION(:), ALLOCATABLE ::                                         &
-          ddzu, dd2zu, dzu, ddzw, dzw, km_damp_x, km_damp_y, l_grid, pt_init,  &
-          q_init, rdf, sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,      &
-          v_nzb_p1_for_vfc, zu, zw
+          ddzu, dd2zu, dzu, ddzw, dzw, hyp, km_damp_x, km_damp_y, l_grid,      &
+          pt_init, q_init, rdf, sa_init, ug, u_init, u_nzb_p1_for_vfc, vg,     &
+          v_init, v_nzb_p1_for_vfc, zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
@@ -102,10 +103,10 @@
 
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
-          qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, saswst_1, shf_1,     &
-          shf_2, tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
+          qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, saswsb_1, saswst_1,  &
+          shf_1, shf_2, tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
 
     REAL, DIMENSION(:,:), POINTER ::                                           &
-          qsws, qsws_m, qswst, qswst_m, rif, rif_m, saswst, shf, shf_m, tswst, &
-          tswst_m, usws, usws_m, vsws, vsws_m
+          qsws, qsws_m, qswst, qswst_m, rif, rif_m, saswsb, saswst, shf,       &
+          shf_m, tswst, tswst_m, usws, usws_m, vsws, vsws_m
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
@@ -118,6 +119,6 @@
     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, r, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1, v_2,  &
-          v_3, vpt_1, vpt_2, w_1, w_2, w_3
+          q_2, q_3, ql_1, ql_2, rho, 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 ::                                         &
@@ -167,8 +168,7 @@
               mass_of_solute, molecular_weight_of_solute,                      &
               prec_time_const = 0.001, ql_crit = 0.0005, rho_l = 1.0E3,        &
-              r_d = 287.0, r_v = 461.51, rho_surface,                          &
-              thermal_conductivity_l = 2.43E-2
-
-    REAL, DIMENSION(:), ALLOCATABLE :: hydro_press, pt_d_t, t_d_pt
+              r_d = 287.0, r_v = 461.51, thermal_conductivity_l = 2.43E-2
+
+    REAL, DIMENSION(:), ALLOCATABLE   ::  pt_d_t, t_d_pt
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::  precipitation_amount, &
@@ -241,4 +241,5 @@
                             bc_q_b = 'dirichlet', bc_q_t = 'neumann', &
                             bc_s_b = 'dirichlet', bc_s_t = 'neumann', &
+                            bc_sa_t = 'neumann', &
                             bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
                             dissipation_1d = 'as_in_3d_model', &
@@ -272,5 +273,5 @@
                 dvrp_filecount = 0, dz_stretch_level_index, gamma_mg, &
                 grid_level, ibc_e_b, ibc_p_b, ibc_p_t, ibc_pt_b, ibc_pt_t, &
-                ibc_q_b, ibc_q_t, ibc_uv_b, ibc_uv_t, &
+                ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b, ibc_uv_t, &
                 inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
                 intermediate_timestep_count, intermediate_timestep_count_max, &
@@ -302,4 +303,5 @@
                 conserve_volume_flow = .FALSE., constant_diffusion = .FALSE., &
                 constant_heatflux = .TRUE., constant_top_heatflux = .TRUE., &
+                constant_top_salinityflux = .TRUE., &
                 constant_waterflux = .TRUE., create_disturbances = .TRUE., &
                 cut_spline_overshoot = .TRUE., &
@@ -333,4 +335,5 @@
              averaging_interval = 0.0, averaging_interval_pr = 9999999.9, &
              averaging_interval_sp = 9999999.9, bc_pt_t_val, bc_q_t_val, &
+             bottom_salinityflux = 0.0, &
              building_height = 50.0, building_length_x = 50.0, &
              building_length_y = 50.0, building_wall_left = 9999999.9, &
@@ -363,7 +366,7 @@
              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, 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, &
+             restart_time = 9999999.9, rho_ref, 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, &
              skip_time_data_output = 0.0, skip_time_data_output_av = 9999999.9,&
              skip_time_dopr = 9999999.9, skip_time_dosp = 9999999.9, &
@@ -380,5 +383,5 @@
              time_do_sla = 0.0, time_dvrp = 0.0, time_prel = 0.0, &
              time_restart = 9999999.9, time_run_control = 0.0, &
-             top_heatflux = 9999999.9, top_salinityflux = 0.0, &
+             top_heatflux = 9999999.9, top_salinityflux = 9999999.9, &
              ug_surface = 0.0, u_gtrans = 0.0, &
              ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, &
@@ -393,7 +396,7 @@
              q_vertical_gradient_level(10) = -1.0, &
              s_vertical_gradient(10) = 0.0, &
-             s_vertical_gradient_level(10) = -1.0,
+             s_vertical_gradient_level(10) = -1.0, &
              sa_vertical_gradient(10) = 0.0, &
-             sa_vertical_gradient_level(10) = -1.0, threshold(20) = 0.0, &
+             sa_vertical_gradient_level(10) = -9999999.9, threshold(20) = 0.0, &
              tsc(10) = (/ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), &
              ug_vertical_gradient(10) = 0.0, &
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 94)
+++ /palm/trunk/SOURCE/parin.f90	(revision 95)
@@ -4,6 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! +ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level,
-! top_salinity_flux in inipar,
+! +bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient,
+! sa_vertical_gradient_level, top_salinityflux in inipar,
 ! sa_init is allocated
 !
@@ -63,5 +63,6 @@
     NAMELIST /inipar/  adjust_mixing_length, alpha_surface, bc_e_b, bc_lr, &
                        bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b, &
-             bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t, building_height, &
+             bc_q_t,bc_s_b, bc_s_t, bc_sa_t, bc_uv_b, bc_uv_t, &
+             bottom_salinityflux, building_height, &
              building_length_x, building_length_y, building_wall_left, &
              building_wall_south, cloud_droplets, cloud_physics, &
@@ -87,5 +88,5 @@
              surface_waterflux, s_surface, s_surface_initial_change, &
              s_vertical_gradient, s_vertical_gradient_level, top_heatflux, &
-             top_salinity_flux, timestep_scheme, topography, ug_surface, &
+             top_salinityflux, timestep_scheme, topography, ug_surface, &
              ug_vertical_gradient, ug_vertical_gradient_level, ups_limit_e, &
              ups_limit_pt, ups_limit_u, ups_limit_v, ups_limit_w, &
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 94)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 95)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! prognostic equation for salinity
 ! new argument zw in calls of diffusion_e
 !
@@ -414,4 +415,93 @@
 
 !
+!-- If required, compute prognostic equation for salinity
+    IF ( ocean )  THEN
+
+       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
+
+!
+!--    sa-tendency terms with communication
+       sat = tsc(1)
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
+!--          switched on. Thus:
+             sat = 1.0
+             sbt = 1.0
+          ENDIF
+          tend = 0.0
+          CALL advec_s_bc( sa, 'sa' )
+       ELSE
+          IF ( tsc(2) /= 2.0 )  THEN
+             IF ( scalar_advec == 'ups-scheme' )  THEN
+                tend = 0.0
+                CALL advec_s_ups( sa, 'sa' )
+             ENDIF
+          ENDIF
+       ENDIF
+
+!
+!--    sa terms with no communication
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+!
+!--          Tendency-terms
+             IF ( scalar_advec == 'bc-scheme' )  THEN
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
+                                  tend )
+             ELSE
+                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) THEN
+                   tend(:,j,i) = 0.0
+                   CALL advec_s_pw( i, j, sa )
+                ELSE
+                   IF ( scalar_advec /= 'ups-scheme' )  THEN
+                      tend(:,j,i) = 0.0
+                      CALL advec_s_up( i, j, sa )
+                   ENDIF
+                ENDIF
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
+                                  tend )
+             ENDIF
+       
+             CALL user_actions( i, j, 'sa-tendency' )
+
+!
+!--          Prognostic equation for salinity
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                sa_p(k,j,i) = sat * sa(k,j,i) +                                &
+                              dt_3d * (                                        &
+                                     sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
+                                      ) -                                      &
+                              tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
+                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
+             ENDDO
+
+!
+!--          Calculate tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tsa_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count < &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                      5.3125 * tsa_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
+
+    ENDIF
+
+!
 !-- If required, compute prognostic equation for total water content / scalar
     IF ( humidity  .OR.  passive_scalar )  THEN
@@ -647,6 +737,7 @@
 ! be called for the standard Piascek-Williams advection scheme.
 !
-! The call of this subroutine is embedded in two DO loops over i and j, thus
-! communication between CPUs is not allowed in this subroutine.
+! Here the calls of most subroutines are embedded in two DO loops over i and j,
+! so communication between CPUs is not allowed (does not make sense) within
+! these loops.
 !
 ! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
@@ -883,4 +974,51 @@
                    ENDDO
                 ENDIF
+             ENDIF
+
+!
+!--          If required, compute prognostic equation for salinity
+             IF ( ocean )  THEN
+
+!
+!--             Tendency-terms for salinity
+                tend(:,j,i) = 0.0
+                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
+                THEN
+                   CALL advec_s_pw( i, j, sa )
+                ELSE
+                   CALL advec_s_up( i, j, sa )
+                ENDIF
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
+                                  tend )
+       
+                CALL user_actions( i, j, 'sa-tendency' )
+
+!
+!--             Prognostic equation for salinity
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   sa_p(k,j,i) = tsc(1) * sa(k,j,i) +                          &
+                                 dt_3d * (                                     &
+                                  tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
+                                         ) -                                   &
+                                tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
+                   IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
+                ENDDO
+
+!
+!--             Calculate tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tsa_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                         5.3125 * tsa_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
              ENDIF
 
@@ -1363,4 +1501,96 @@
 
 !
+!-- If required, compute prognostic equation for salinity
+    IF ( ocean )  THEN
+
+       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
+
+!
+!--    sa-tendency terms with communication
+       sat = tsc(1)
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
+!--          switched on. Thus:
+             sat = 1.0
+             sbt = 1.0
+          ENDIF
+          tend = 0.0
+          CALL advec_s_bc( sa, 'sa' )
+       ELSE
+          IF ( tsc(2) /= 2.0 )  THEN
+             IF ( scalar_advec == 'ups-scheme' )  THEN
+                tend = 0.0
+                CALL advec_s_ups( sa, 'sa' )
+             ENDIF
+          ENDIF
+       ENDIF
+
+!
+!--    Scalar/q-tendency terms with no communication
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+          CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, tend )
+       ELSE
+          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
+             tend = 0.0
+             CALL advec_s_pw( sa )
+          ELSE
+             IF ( scalar_advec /= 'ups-scheme' )  THEN
+                tend = 0.0
+                CALL advec_s_up( sa )
+             ENDIF
+          ENDIF
+          CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, tend )
+       ENDIF
+       
+       CALL user_actions( 'sa-tendency' )
+
+!
+!--    Prognostic equation for salinity
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                sa_p(k,j,i) = sat * sa(k,j,i) +                                &
+                              dt_3d * (                                        &
+                                     sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
+                                      ) -                                      &
+                              tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
+                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Calculate tendencies for the next Runge-Kutta step
+       IF ( timestep_scheme(1:5) == 'runge' )  THEN
+          IF ( intermediate_timestep_count == 1 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tsa_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSEIF ( intermediate_timestep_count < &
+                   intermediate_timestep_count_max )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                      5.3125 * tsa_m(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
+
+       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
+
+    ENDIF
+
+!
 !-- If required, compute prognostic equation for total water content / scalar
     IF ( humidity  .OR.  passive_scalar )  THEN
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 94)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +sa, saswst
+! +sa, saswsb, saswst
 !
 ! Former revisions:
@@ -284,4 +284,6 @@
           CASE ( 'sa' )
              READ ( 13 )  sa
+          CASE ( 'saswsb' )
+             READ ( 13 )  saswsb
           CASE ( 'saswst' )
              READ ( 13 )  saswst
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 94)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 95)
@@ -4,6 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! +ocean, sa_init, sa_surface, sa_vertical_gradient,
-! sa_vertical_gradient_level, top_salinity_flux
+! +bc_sa_t, ocean, sa_init, sa_surface, sa_vertical_gradient,
+! sa_vertical_gradient_level, bottom/top_salinity_flux
 !
 ! Former revisions:
@@ -161,8 +161,12 @@
           CASE ( 'bc_s_t' )
              READ ( 13 )  bc_s_t
+          CASE ( 'bc_sa_t' )
+             READ ( 13 )  bc_sa_t
           CASE ( 'bc_uv_b' )
              READ ( 13 )  bc_uv_b
           CASE ( 'bc_uv_t' )
              READ ( 13 )  bc_uv_t
+          CASE ( 'bottom_salinityflux' )
+             READ ( 13 )  bottom_salinityflux
           CASE ( 'building_height' )
              READ ( 13 )  building_height
Index: /palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- /palm/trunk/SOURCE/swap_timelevel.f90	(revision 94)
+++ /palm/trunk/SOURCE/swap_timelevel.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Swaping of salinity
 !
 ! Former revisions:
@@ -106,8 +106,11 @@
              pt => pt_1;  pt_p => pt_2
              IF ( .NOT. constant_diffusion )  THEN
-                e => e_1;  e_p => e_2
-             ENDIF
-             IF ( humidity  .OR.  passive_scalar )  THEN
-                q => q_1;  q_p => q_2
+                e => e_1;    e_p => e_2
+             ENDIF
+             IF ( ocean )  THEN
+                sa => sa_1;  sa_p => sa_2
+             ENDIF
+             IF ( humidity  .OR.  passive_scalar )  THEN
+                q => q_1;    q_p => q_2
              ENDIF
 
@@ -152,8 +155,11 @@
              pt => pt_2;  pt_p => pt_1
              IF ( .NOT. constant_diffusion )  THEN
-                e => e_2;  e_p => e_1
-             ENDIF
-             IF ( humidity  .OR.  passive_scalar )  THEN
-                q => q_2;  q_p => q_1
+                e => e_2;    e_p => e_1
+             ENDIF
+             IF ( ocean )  THEN
+                sa => sa_2;  sa_p => sa_1
+             ENDIF
+             IF ( humidity  .OR.  passive_scalar )  THEN
+                q => q_2;    q_p => q_1
              ENDIF
 
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 94)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Ghostpoint exchange for salinity and density
 !
 ! Former revisions:
@@ -150,4 +150,8 @@
           CALL exchange_horiz( pt_p )
           IF ( .NOT. constant_diffusion       )  CALL exchange_horiz( e_p )
+          IF ( ocean )  THEN
+             CALL exchange_horiz( sa_p )
+             CALL exchange_horiz( rho )
+          ENDIF
           IF ( humidity  .OR.  passive_scalar )  CALL exchange_horiz( q_p )
           IF ( cloud_droplets )  THEN
Index: /palm/trunk/SOURCE/user_interface.f90
===================================================================
--- /palm/trunk/SOURCE/user_interface.f90	(revision 94)
+++ /palm/trunk/SOURCE/user_interface.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! user action for salinity added
 !
 ! Former revisions:
@@ -393,4 +393,7 @@
 
 
+          CASE ( 'sa-tendency' )
+
+
           CASE ( 'e-tendency' )
 
@@ -443,4 +446,7 @@
 
           CASE ( 'pt-tendency' )
+
+
+          CASE ( 'sa-tendency' )
 
 
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 94)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 95)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +sa, saswst
+! +sa, saswsb, saswst
 !
 ! Former revisions:
@@ -134,4 +134,5 @@
     IF ( ocean )  THEN
        WRITE ( 14 )  'sa                  ';  WRITE ( 14 )  sa
+       WRITE ( 14 )  'saswsb              ';  WRITE ( 14 )  saswsb
        WRITE ( 14 )  'saswst              ';  WRITE ( 14 )  saswst
     ENDIF
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 94)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 95)
@@ -4,6 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! +ocean, sa_init, sa_surface, sa_vertical_gradient,
-! sa_vertical_gradient_level, top_salinity_flux
+! +bc_sa_t, ocean, sa_init, sa_surface, sa_vertical_gradient,
+! sa_vertical_gradient_level, bottom/top_salinity_flux
 !
 ! Former revisions:
@@ -107,8 +107,12 @@
     WRITE ( 14 )  'bc_s_t                        '
     WRITE ( 14 )  bc_s_t
+    WRITE ( 14 )  'bc_sa_t                       '
+    WRITE ( 14 )  bc_sa_t
     WRITE ( 14 )  'bc_uv_b                       '
     WRITE ( 14 )  bc_uv_b
     WRITE ( 14 )  'bc_uv_t                       '
     WRITE ( 14 )  bc_uv_t
+    WRITE ( 14 )  'bottom_salinityflux           '
+    WRITE ( 14 )  bottom_salinityflux
     WRITE ( 14 )  'building_height               '
     WRITE ( 14 )  building_height
